Utilisation de la bibliothèque PROJ.4 pour passer des coordonnées du système de coordonnées local au système de coordonnées global à l'aide de points de contrôle au sol?

J'ai un nuage de points dont les coordonnées sont par rapport à un système de coordonnées local. J'ai également des points de contrôle au sol avec des valeurs GPS. Puis-je convertir ces coordonnées locales en un système de coordonnées global à l'aide de PROJ.4 ou de toute autre bibliothèque?

Tout code en Python pour le problème mentionné ci-dessus serait d'une grande aide.

Vous semblez chercher à effectuer une transformation affine entre votre système de coordonnées local et un système de coordonnées géoréférencées.

Affine transforme sous tous les systèmes de coordonnées et peut être représenté par l'équation de la matrice ci-dessous.

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

Cependant, vous avez un problème en deux étapes.

1. Trouvez la matrice de transformation à partir d'appariements connus de coordonnées d'entrée et de sortie (vos points GPS et leurs emplacements respectifs dans votre grille définie localement).
2. Utilisez cette matrice de transformation pour géoréférencer votre nuage de points.

Le projet 4 excelle au n ° 2: le transfert entre des systèmes de coordonnées géoréférencées avec des matrices de transformation connues. À ma connaissance, il ne peut pas être utilisé pour trouver une matrice de transformation à partir de données ponctuelles. Cependant, vous pouvez faire tout cela facilement en utilisant une algèbre linéaire légère (une inversion de matrice des moindres carrés) dans Numpy. J'ai utilisé une version de cette classe pour réduire les données de plusieurs études sur le terrain:

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

Il peut être utilisé tel quel:

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

projected_coordinatesest maintenant en WGS84, UTM, ou tout autre système de coordonnées enregistré par votre GPS. Une caractéristique majeure de cette méthode est qu'elle peut être utilisée avec n'importe quel nombre de points de liaison (3 ou plus) et gagne en précision à mesure que plus de points de liaison sont utilisés. Vous trouvez essentiellement le meilleur ajustement à travers tous vos points de liaison.

La fonction Grass transform fait exactement ce dont vous avez besoin, même si elle n'est pas basée sur python ou proj comme demandé:

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

Il est toujours plus facile d'identifier le système de coordonnées local, comme nous l'avons fait ici:

Projection stéréographique d'un ellipsoïde WGS84 sur un plan [python]

GDAL est maintenant capable de transformer des données vectorielles en utilisant des points GCP.

J'étais coincé dans le même problème il y a quelques semaines, j'ai trouvé un script python qui peut aider. Solution originale d' ici

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)