1

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?

2

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

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

  • 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. – Mike T Aug 5 '15 at 3:02
  • @MikeT Is there any way I can do the transform outside of the range of that projection? – Justin Dearing Aug 5 '15 at 12:50
  • 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 Aug 5 '15 at 12:50
  • You were actually right, that transform was invalid. Why did it work on another system? Can proj4 be compiled to not check bounds? – Justin Dearing Aug 6 '15 at 2:21
  • 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. – Mintx Aug 12 '15 at 18:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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