I have a point cloud whose coordinates are with respect to a local coordinate system. I also have ground control points with GPS values. Can I convert these local coordinates to a global coordinate system using PROJ.4 or any other library?
Any code in Python for the above stated problem would be a great help.
You seem to be looking to conduct an affine transformation between your local coordinate system and a georeferenced coordinate system.
Affine transforms underly all coordinate systems and can be represented by the matrix equation below.
|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)
However, you have a two-step problem.
Proj.4 excels at #2: transferring between georeferenced coordinate systems with known transformation matrices. It can't to my knowledge be used to find a transformation matrix from point data. However, you can do the entire thing easily by using some light linear algebra (a least-squares matrix inversion) in Numpy. I've used a version of this class for reducing data from several field studies:
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
It can be used as such:
transform = Affine.from_tiepoints(gps_points_local,gps_points_geo)
projected_data = transform.transform(local_point_cloud)
projected_coordinates is now in WGS84, UTM, or whatever coordinate system your recorded by your GPS. A major feature of this method is that it can be used with any number of tie points (3 or more) and gains accuracy the more tie points are used. You're essentially finding the best fit through all of your tie points.
I was stuck in the same problem a few weeks ago, I figured out a python script that can help. Original solution from here
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)
Grass transform function does exactly what you need, even though it is not python or proj based as requested:
It is always easier to identify the local coordinate system, as we did here:
Stereographic projection of WGS84 ellipsoid on a plane[python]
GDAL is now able to transform vector data using GCP points.
