Использование библиотеки PROJ.4 для преобразования из локальной системы координат в глобальную систему координат с использованием наземных контрольных точек?

9

У меня есть облако точек, координаты которого относительно локальной системы координат. У меня также есть наземные контрольные точки со значениями GPS. Могу ли я преобразовать эти локальные координаты в глобальную систему координат, используя PROJ.4 или любую другую библиотеку?

Любой код на Python для вышеупомянутой проблемы был бы большим подспорьем.

user18953
источник
Какой-то код ожидается?
huckfinn
Координаты GPS обычно WGS84, поэтому они, вероятно, уже глобальные. Если наземные контрольные точки находятся в локальной проекции с данными, отличными от GPS (например, NAD83), данные должны быть преобразованы. Насколько мне известно, PROJ4 поддерживает базовые сдвиги.
Ойвинд
Вот аналогичный вопрос, но с гораздо большей детализацией: gis.stackexchange.com/questions/357910 .
trusktr

Ответы:

7

Похоже, вы пытаетесь провести аффинное преобразование между вашей локальной системой координат и системой координат с географической привязкой.

Аффинные преобразования лежат в основе всех систем координат и могут быть представлены матричным уравнением ниже.

|x_1 y_1 1| |a d|   |x'_1 y'_1|
|x_2 y_2 1| |b e| = |x'_2 y'_2|
|x_3 y_3 1| |c f|   |x'_3 y'_3|
input     transform.  output
coords    matrix      coords
(n x 3)   (3 x 2)     (n x 2)

Однако у вас есть двухэтапная проблема.

  1. Найдите матрицу преобразования из известных пар входных и выходных координат (ваши точки GPS и их соответствующие местоположения в вашей локальной сетке).
  2. Используйте эту матрицу преобразования для привязки вашего облака точек.

Proj.4 превосходит на # 2: перенос между системами координат с привязками с известными матрицами преобразования. Насколько мне известно, он не может быть использован для поиска матрицы преобразования из точечных данных. Тем не менее, вы можете сделать все это легко с помощью некоторой легкой линейной алгебры (инверсия матрицы наименьших квадратов) в Numpy. Я использовал версию этого класса для сокращения данных из нескольких полевых исследований:

import numpy as N 

def augment(a):
    """Add a final column of ones to input data"""
    arr = N.ones((a.shape[0],a.shape[1]+1))
    arr[:,:-1] = a
    return arr

class Affine(object):
    def __init__(self, array=None):
        self.trans_matrix = array

    def transform(self, points):
        """Transform locally projected data using transformation matrix"""
        return N.dot(augment(N.array(points)), self.trans_matrix)

    @classmethod
    def from_tiepoints(cls, fromCoords, toCoords):
        "Produce affine transform by ingesting local and georeferenced coordinates for tie points"""
        fromCoords = augment(N.array(fromCoords))
        toCoords = N.array(toCoords)
        trans_matrix, residuals, rank, sv = N.linalg.lstsq(fromCoords, toCoords)

        affine =  cls(trans_matrix) # Setup affine transform from transformation matrix
        sol = N.dot(fromCoords,affine.trans_matrix) # Compute model solution
        print "Pixel errors:"
        print (toCoords - sol)
        return affine

Может использоваться как таковой:

transform = Affine.from_tiepoints(gps_points_local,gps_points_geo)
projected_data = transform.transform(local_point_cloud)

projected_coordinatesтеперь в WGS84, UTM или в любой другой системе координат, записанной вашим GPS. Главной особенностью этого метода является то, что его можно использовать с любым количеством точек привязки (3 или более), и он получает точность при использовании большего количества точек привязки. По сути, вы находите наилучшее соответствие через все ваши точки привязки.

Давен Куинн
источник
Здравствуйте! Вы упоминаете, что Proj (Proj4) не может обрабатывать пользовательские части преобразования? Означает ли это, что технически нет простого ответа Proj на вопрос по адресу gis.stackexchange.com/questions/357910 ?
trusktr
0

Всегда проще определить локальную систему координат, как мы это делали здесь:

Стереографическая проекция эллипсоида WGS84 на плоскость [питон]

GDAL теперь может преобразовывать векторные данные, используя точки GCP.

Andrej
источник
Здравствуйте! Я новичок в этом деле. Является ли «стереографическая» проекция на самом деле тем, что мне нужно использовать для моего вопроса на gis.stackexchange.com/questions/357910 ?
trusktr
0

Я застрял в той же проблеме несколько недель назад, я нашел сценарий Python, который может помочь. Оригинальное решение здесь

import pyproj
import math
import numpy as np
from statistics import mean
import scipy.optimize as optimize

#This function converts the numbers into text
def text_2_CRS(params):
    # print(params)  # <-- you'll see that params is a NumPy array
    x_0, y_0, gamma, alpha, lat_0, lonc = params # <-- for readability you may wish to assign names to the component variables
    pm = '+proj=omerc +lat_0='+ str(lat_0) +' +lonc='+ str(lonc) +' +alpha=' + str(alpha) + ' +gamma=' + str(
        gamma) + ' +k=0.999585495 +x_0=' + str(x_0) + ' +y_0=' + str(y_0) + ' +ellps=GRS80 +units=m +no_defs'
    return pm

#Optimisation function
def convert(params):
    pm = text_2_CRS(params)
    trans_points = []
    #Put your control points in mine grid coordinates here
    points_local = [[5663.648, 7386.58],
                    [20265.326, 493.126],
                    [1000, -10000],
                    [-1000, -10000],
                    [1331.817, 2390.206],
                    [5794, -1033.6],
                    ]
    # Put your control points here mga here
    points_mga = [[567416.145863305, 7434410.3451835],
                  [579090.883705669, 7423265.25196681],
                  [557507.390559793, 7419390.6658927],
                  [555610.407664593, 7420021.64968145],
                  [561731.125709093, 7431037.98474379],
                  [564883.285081307, 7426382.75146683],
                  ]
    for i in range(len(points_local)):
        #note that EPSG:28350 is MGA94 Zone 50
        trans = pyproj.transform(pyproj.Proj(pm), pyproj.Proj("EPSG:28350"), points_local[i][0], points_local[i][1])
        trans_points.append(trans)
    error = []
    #this finds the difference between the control points
    for i in range(len(points_mga)):
        x1 = trans_points[i][0]
        y1 = trans_points[i][1]
        x2 = points_mga[i][0]
        y2 = points_mga[i][1]
        error.append(math.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2))

    print("Current Params are: ")
    with np.printoptions(precision=3, suppress=True):
        print(params)
    print("Current average error is: " + str(mean(error)) + " meters")
    print("String to use is: " + pm)
    print('')

    return mean(error)


#Add your inital guess
x_0 = 950
y_0 = -1200
gamma = -18.39841101
alpha=-0
lat_0 = -23.2583926082939
lonc = 117.589084840039


#define your control points
points_local = [[5663.648,7386.58],
          [20265.326,493.126],
          [1000,-10000],
          [-1000,-10000],
          [1331.817,2390.206],
          [5794,-1033.6],
          ]

points_mga = [[567416.145863305,7434410.3451835],
          [579090.883705669,7423265.25196681],
          [557507.390559793,7419390.6658927],
          [555610.407664593,7420021.64968145],
          [561731.125709093,7431037.98474379],
          [564883.285081307,7426382.75146683],
          ]


params = [x_0, y_0, gamma,alpha, lat_0, lonc]

error = convert(params)

print(error)

result = optimize.minimize(convert, params, method='Powell')
if result.success:
    fitted_params = result.x
    print(fitted_params)
else:
    raise ValueError(result.message)
Гарт Купер
источник