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:
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 ] ] } 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?
أكثر...
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:
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 ] ] } 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?
أكثر...