Projections and transformation

On his blog Richie Carmichael (First edition Second edition) showed how to use the ESRI Projection Engine from .NET. In this article I am going to do the same with Python.
My original problem was that I needed to convert some coordinates from ED50 Projected and ED50 geographic to ETRS89 projected for the UTM zone 31.
First of all we will need the pe.dll provided by ESRI which can be found in the BIN folder of the ArcGIS installation. Copy this dll in the same folder as your Python script.
To use this dll from Python you can use the ctypes package. This extension started has a separate project in Python 2.3 but it has been fully incorporated in Python 2.5. So if you're still using Python 2.4 that shipped with ArcGIS 9.2 You will have to download and install ctypes. This extensions creates the possibility to easily use C dll's from within Python.
Now lets start writing some code. First we import the ctypes package.
    import ctypes
The pe.dll is accessed in the following way.
    pe = ctypes.windll.pe
To transform from one projection to another we need to first create a projection object for the input spatial reference and the output spatial reference. Then we also need a geographic transformation object. We do this by calling the appropriate factory functions which can be found in the SDE C API documentation on EDN.
    input_proj = pe.pe_factory_projcs(input_proj_code)
    transformation = pe.pe_factory_geogtran(transformation_code)
    output_proj = pe.pe_factory_projcs(output_proj_code)
To lookup the code of your projection you will need to take a look at the following sites page 1 page 2 page 3 and page 4. The code of your geographic transformation can be found on the this sites page 1 page 2 and page 3.
To transform points with the projection engine we need to pass the points as a 2 dimensional array of points. In order to do this I created a small Python point class and methods to convert from a list of points to a 2 dimensional C array and back to a list of points. If you want more information about this I suggest you take a look at the ctypes documentation (tutorial reference) and also have a look at the needed datatypes for the pe functions.
    class Point:
        def __init__(self, x, y):
            self.x = x
            self.y = y
        def __str__(self):
            return 'x : ' + str(self.x) + ' y : ' + str(self.y)
    
    def _ctypes_to_points(point_array):
        points = []
        for point in point_array:
            point = Point(point[0], point[1])
            points.append(point)
        return points

    def _points_to_ctypes(points):
        array_type = ctypes.ARRAY(self.point_type, len(points))
        ctype_points = []
        for point in points:
            ctype_point = self.point_type(point.x, point.y)
            ctype_points.append(ctype_point)
        return array_type(*ctype_points)
Now we can create the function for doing the actual transformation from input projected reference system to the output projected reference system. As you can see we first need to convert the projected coordinates to geographic coordinates. Then transform these geographic coordinates with the give transformation. Finally we convert the geographic coordinates back to projected coordinates.
def transform(self, points):
    points_count = len(points)
    point_array = self._points_to_ctypes(points)
    pe.pe_proj_to_geog(self.input_proj, points_count, 
                       point_array)
    pe.pe_geog1_to_geog2(self.transformation, points_count, 
                         point_array, None)
    pe.pe_geog_to_proj(self.output_proj, points_count, 
                       point_array)
    return self._ctypes_to_points(point_array)
The complete code looks like this. I've put the projection related code in a class and at the bottom I added some test code.
import ctypes

class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y
    def __str__(self):
        return 'x : ' + str(self.x) + ' y : ' + str(self.y)

class Proj:
    def __init__(self, input_proj_code, transformation_code, 
        output_proj_code):
        self.pe = ctypes.windll.pe
        self.input_proj = self.pe.pe_factory_projcs(input_proj_code)
        self.output_proj = self.pe.pe_factory_projcs(output_proj_code)
        self.transformation = self.pe.pe_factory_geogtran(transformation_code)
        self.point_type = ctypes.c_double * 2

    def _points_to_ctypes(self, points):
        array_type = ctypes.ARRAY(self.point_type, len(points))
        ctype_points = []
        for point in points:
            ctype_point = self.point_type(point.x, point.y)
            ctype_points.append(ctype_point)
        return array_type(*ctype_points)

    def _ctypes_to_points(self, point_array):
        points = []
        for point in point_array:
            point = Point(point[0], point[1])
            points.append(point)
        return points

    def transform(self, points):
        points_count = len(points)
        point_array = self._points_to_ctypes(points)
        self.pe.pe_proj_to_geog(self.input_proj, points_count, 
                                point_array)
        self.pe.pe_geog1_to_geog2(self.transformation, 
                                  points_count, point_array, None)
        self.pe.pe_geog_to_proj(self.output_proj, points_count, 
                                point_array)
        return self._ctypes_to_points(point_array)

if __name__ == '__main__':
    ed50_proj_code = 23031
    ed50_etrs89_transformation = 1650
    etrs89_proj_code = 25831
    proj = Proj(ed50_proj_code, ed50_etrs89_transformation, etrs89_proj_code)
    points = [Point(585635, 5690770), Point(585654, 5690110)]
    points = proj.transform(points)
    for point in points:
        print point

I hope you liked this post. Have fun.

7 comments:

gene said...

look for Gdal/ogr and python bindings

Samuel Bosch said...

Thanks for the comment, I know ogr have used it before for other stuff like appending shapefiles and reading MapInfo files from C#. I used the ESRI pe.dll because I didn't want to install extra software at the client side and wanted to be sure that the results where the same as when you use the ESRI projection tool.

Scott said...

Great post! But I am having an issue you may be able to help me with. I'm trying to convert from GCS NAD27 to PCS UTM NAD83 17N but I can't get the transformation to work properly.

Basically, I removed the pe_proj_to_geog from the transform function. The script runs fine but the results don't seem to take the transformation into account. In fact it doesn't seem to matter if I do the transformation or not, I always get the same results.

Do you have any suggestions here? Any advice would be greatly appreciated. Thanks!

-Scott

Samuel said...

Did you try the projection in ArcMap or ArcCatalog ?
I'm not familiar with the American projection systems but I think you will first have to do a geographical transformation from NAD27 to NAD83 and the from NAD83 to the projected NAD83 17N.
Can you mail me some test data and your script ?

Succes!

Samuel

Scott said...

I've tried the projection in ArcCatalog and ArcObjects and get the same results. Unfortunately, I need this to work in Python. What's the best way to get the script to you? You can email me at swiegand@eqt.com. Thanks!

Anonymous said...

What I am trying to do is to load into memory a .dll module. It has been successfully loaded.

>>> mylib = CDLL("d:\\SHImaps\AS13910\ArealInterpolator.dll")
>>> print mylib


Manual use of this module will need input parameters (GIS map files) and the module will produce a new file.

I do not know what to do after loading it into the memory.

I am looking for excellent example codes and instructions.

I am trying to to programmatically run and use .dll file in Python.

Sincerely,

David
davidgshi@yahoo.co.uk

Samuel Bosch said...

Hi David,

I just took a brief look at the ArealInterpolator and sadly enough I couldn't find the documentation on the internet. But I found out that you it is a toolbar with a command for ArcMap. You can add it to ArcMap by clicking Tools -> Customize -> click Add from file (browse to ArealInterpolator.dll). Once you have done this you will see the Areal Interpolator in the list of toolbars. Then its up to you to find out what it exactly does.

I don't think it can be used with Python.

Kind regards,

Samuel