I have a shapefile that has a .prj file, with the following wkt contents:
I am using the following python code to convert that to a pyproj projection:
prjText = open( prjPath, 'r').read() srs = osr.SpatialReference() if ( srs.ImportFromWkt( prjText ) ): return ( False, "error importing .prj information from " + prjPath ) inProjection = pyproj.Proj( srs.ExportToProj4() )
I then create another projection:
outProjection = pyproj.Proj( "+proj=latlong +datum=WGS84" )
and then go through each point
p from the shapefile and convert it:
print "before transform: " + str( p[ 0 ] ) + ", " + str( p[ 1 ] ) p = pyproj.transform( inProjection, outProjection, p[ 0 ], p[ 1 ] ) print " -> after transform: " + str( p[ 0 ] ) + ", " + str( p[ 1 ] )
when I print the srs.ExportToProj4() text, it gives:
+proj=tmerc +lat_0=40 +lon_0=-78.58333333333333 +k=0.9999375 +x_0=152400.3048006096 +y_0=0 +ellps=clrk66 +to_meter=0.3048006096012192 +no_defs
Furthermore, using another program that I have, I import the same shapefile without using the .prj file. To do this, I specify that the projection is "State Plane", and specifically that it is "New York West NAD27 (ftUS)". This other program then determines that that means the same as epsg:32017, and indeed when I look up epsg:32017 I see that its srs string matches the one above.
Thus, it seems like my inProjection variable (above) should be set up correctly. However, when I run the code above to transform each of the points, I get, for the first point:
before transform: 759232.003438, 1149854.52147 -> after transform: -70.1071453809, 50.036825617
That lon/lat, approximately (-70, 50), is not in the right place. By importing the file into either the above-mentioned program, or into ArcGIS Explorer Desktop, I see the point (correctly) in northern New York, just south of Lake Ontario. But (-70, 50) is up in central Quebec. So it's off by like 1000km to the north and east.
So the questions are:
(1) Does it seems likely something is going wrong in pyproj (I'm using version 1.8.9)? (since according to these other programs my srs string is correct, so everything is good up to the point I use pyproj)
(2) Based on the wkt and srs string shown above, what lon/lat should I expect the the shapefile point (759232.003438, 1149854.52147) to correspond to? Is that point the x/y in meters relative to the lon_0 and lat_0 of the srs string? If that were the case, I would indeed expect a point somewhere up in Quebec, and not in northern New York. However, since both of these other programs are correctly mapping the point to northern New York, that must not be the correct interpretation of the file coordinates with respect to the file's projection?
Thanks for any help.