Używasz biblioteki PROJ.4 do transformacji z lokalnych współrzędnych układu współrzędnych na globalny układ współrzędnych z wykorzystaniem punktów kontroli na ziemi?

Mam chmurę punktów, której współrzędne odnoszą się do lokalnego układu współrzędnych. Mam także naziemne punkty kontrolne z wartościami GPS. Czy mogę przekonwertować te lokalne współrzędne na globalny układ współrzędnych za pomocą PROJ.4 lub innej biblioteki?

Każdy kod w Pythonie dla wyżej wymienionego problemu byłby bardzo pomocny.

Wygląda na to, że chcesz przeprowadzić afiniczną transformację między lokalnym układem współrzędnych a układem współrzędnych georeferencyjnych.

Afinit przekształca zasadniczo wszystkie układy współrzędnych i można je przedstawić za pomocą poniższego równania macierzowego.

|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)

Masz jednak problem dwuetapowy.

1. Znajdź macierz transformacji ze znanych par współrzędnych wejściowych i wyjściowych (punkty GPS i ich odpowiednie lokalizacje w zdefiniowanej lokalnie siatce).
2. Użyj tej macierzy transformacji do georeferencji chmury punktów.

Proj.4 wyróżnia się na nr 2: przenoszenie między układami współrzędnych georeferencyjnych ze znanymi matrycami transformacji. Według mojej wiedzy nie można go użyć do znalezienia macierzy transformacji z danych punktowych. Możesz jednak z łatwością zrobić to wszystko, używając lekkiej algebry liniowej (odwrócenie macierzy najmniejszych kwadratów) w Numpy. Użyłem wersji tej klasy do zmniejszenia danych z kilku badań terenowych:

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

Może być używany jako taki:

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

projected_coordinatesjest teraz w WGS84, UTM lub dowolnym innym układzie współrzędnych zarejestrowanym przez GPS. Główną cechą tej metody jest to, że można jej używać z dowolną liczbą punktów wiążących (3 lub więcej) i zwiększa dokładność, im więcej punktów wiążących jest używanych. Zasadniczo znajdujesz najlepsze dopasowanie we wszystkich punktach remisowych.

Funkcja transformacji trawy wykonuje dokładnie to, czego potrzebujesz, mimo że nie jest oparta na pythonie ani proj, zgodnie z żądaniem:

http://grass.osgeo.org/grass65/manuals/g.transform.html

Zawsze łatwiej jest zidentyfikować lokalny układ współrzędnych, tak jak tutaj:

Rzut stereograficzny elipsoidy WGS84 na płaszczyźnie [python]

GDAL może teraz przekształcać dane wektorowe przy użyciu punktów GCP.

Utknąłem w tym samym problemie kilka tygodni temu, wymyśliłem skrypt w języku Python, który może pomóc. Oryginalne rozwiązanie stąd

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)