Problem doing a coordinate transform with GDAL python

Question

I'm trying to do a concordant transform from the format of the shapefile (3309) to google maps(4326).

This code works on one Centos 6.6 box with the SCL python 2.7 and a custom gdal 2.0 build. Another similarly configured box gives an error. I seem to have the two datums installed:

$ /opt/gdal-custom/bin/gdalsrsinfo epsg:3309
PROJ.4 : '+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD27 +units=m +no_defs '
OGC WKT :
PROJCS["NAD27 / California Albers",
 GEOGCS["NAD27",
 DATUM["North_American_Datum_1927",
 SPHEROID["Clarke 1866",6378206.4,294.9786982138982,
 AUTHORITY["EPSG","7008"]],
 AUTHORITY["EPSG","6267"]],
 PRIMEM["Greenwich",0,
 AUTHORITY["EPSG","8901"]],
 UNIT["degree",0.0174532925199433,
 AUTHORITY["EPSG","9122"]],
 AUTHORITY["EPSG","4267"]],
 PROJECTION["Albers_Conic_Equal_Area"],
 PARAMETER["standard_parallel_1",34],
 PARAMETER["standard_parallel_2",40.5],
 PARAMETER["latitude_of_center",0],
 PARAMETER["longitude_of_center",-120],
 PARAMETER["false_easting",0],
 PARAMETER["false_northing",-4000000],
 UNIT["metre",1,
 AUTHORITY["EPSG","9001"]],
 AXIS["X",EAST],
 AXIS["Y",NORTH],
 AUTHORITY["EPSG","3309"]]
$ /opt/gdal-custom/bin/gdalsrsinfo epsg:4326
PROJ.4 : '+proj=longlat +datum=WGS84 +no_defs '
OGC WKT :
GEOGCS["WGS 84",
 DATUM["WGS_1984",
 SPHEROID["WGS 84",6378137,298.257223563,
 AUTHORITY["EPSG","7030"]],
 AUTHORITY["EPSG","6326"]],
 PRIMEM["Greenwich",0,
 AUTHORITY["EPSG","8901"]],
 UNIT["degree",0.0174532925199433,
 AUTHORITY["EPSG","9122"]],
 AUTHORITY["EPSG","4326"]]

Snippet in question:

 def coordinateTransformation(self,EpsgFrom,EpsgTo):
 from_ogrSpatialReference = osr.SpatialReference()
 to_ogrSpatialReference = osr.SpatialReference()
 from_ogrSpatialReference.ImportFromEPSG(EpsgFrom)
 to_ogrSpatialReference.ImportFromEPSG(EpsgTo)
 retObj = osr.CoordinateTransformation(from_ogrSpatialReference,to_ogrSpatialReference)
 from pprint import pprint
 print("CoordinateTransformation: ")
 pprint(retObj)
 return retObj
. . . .
self.coordTransform = self.coordinateTransformation(4326,3309)
. . .
intersectgeometry = featuregeom.Intersection(poly)
from pprint import pprint
print ("Intersection: {0}".format(intersectgeometry.ExportToJson()))
pprint (self.coordTransform)
intersectgeometry.Transform(self.coordTransform)

Output:

<osgeo.osr.CoordinateTransformation; proxy of <Swig Object of type 'OSRCoordinateTransformationShadow *' at 0x7f9f64652120> >
Intersection: { "type": "LineString", "coordinates": [ [ 1286274.293049275642261, -89857.546843856573105, 0.0 ], [ 1284590.659783720737323, -90115.061295109800994, 0.0 ], [ 1279467.874491626396775, -92807.228362500201911, 0.0 ], [ 1273444.6429856531322, -94958.127677664160728, 0.0 ], [ 1270756.154574947431684, -98174.25142274517566, 0.0 ], [ 1264401.513631681911647, -104188.64924998767674, 0.0 ], [ 1255148.083878455683589, -111527.525126685388386, 0.0 ], [ 1247811.151507917791605, -119020.223048191517591, 0.0 ], [ 1243559.735155994538218, -124702.111474551726133, 0.0 ], [ 1235582.033464746084064, -130480.617263145744801, 0.0 ] ] }
<osgeo.osr.CoordinateTransformation; proxy of <Swig Object of type 'OSRCoordinateTransformationShadow *' at 0x7f9f64652120> >

Stack trace:

 File "/opt/rh/httpd24/root/var/www/MyFile.py", line 380, in MyFunction
 intersectgeometry.Transform(self.coordTransform)
 File "/opt/rh/httpd24/root/var/www/wsgi-virtualenv/lib/python2.7/site-packages/GDAL-2.0.0-py2.7-linux-x86_64.egg/osgeo/ogr.py", line 5236, in Transform
 return _ogr.Geometry_Transform(self, *args)
RuntimeError: OGR Error: General Error

What am I doing wrong?

Answer 1

self.coordTransform = self.coordinateTransformation(4326,3309)

It looks like you have your source and destination coordinate systems mixed up.

• Mike T: Correct. The error message from PROJ.4 that GDAL is not showing is "latitude or longitude exceeded limits". This is the error message from PostGIS or pyproj.
• Justin Dearing: @MikeT Is there any way I can do the transform outside of the range of that projection?
• Justin Dearing: So I actually want to do it this way here. I am transforming my input shape on openbox to the shapefile system. Later I'm going to do it the other way with the output to render that back to the map.
• Justin Dearing: You were actually right, that transform was invalid. Why did it work on another system? Can proj4 be compiled to not check bounds?
• Mintx: You're going from NAD27/California X/Ys to WGS84 lat/longs, so you'll need to know if your source or destination coordinate system requires lat/longs or X/Ys. Clearly -89857 is not a valid latitude, and no transform would be able to handle that value as latitude. You will want that bounds check in place if you want your results to make any sort of sense.