Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
    
Some code expected? –  huckfinn Feb 10 at 22:13
    
GPS coordinates are normally WGS84, so they are probably global already. If the ground control points are in a local projection, with a different datum than the GPS (eg NAD83), the datum must be converted. PROJ4 do support datum shifts as far as I know. –  Oyvind Mar 14 at 9:21
add comment

3 Answers 3

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.

  1. Find the transformation matrix from known pairings of input and output coordinates (your GPS points and their respective locations in your locally-defined grid).
  2. Use this transformation matrix to georeference your point cloud.

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(fromCoordinates))
        toCoords = N.array(toCoordinates)
        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.array) # 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.

share|improve this answer
add comment

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 abel to transform vector data using GCP points. If you can share some of your data, we could work it out.

share|improve this answer
add comment

Grass transform function does exactly what you need, even though it is not python or proj based as requested:

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

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.