projection/pyproj puzzle, and understanding SRS format

Question

I have a shapefile that has a .prj file, with the following wkt contents:

PROJCS["NAD_1927_StatePlane_New_York_West_FIPS_3103",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-78.58333333333333],PARAMETER["Scale_Factor",0.9999375],PARAMETER["Latitude_Of_Origin",40.0],UNIT["Foot_US",0.3048006096012192]]

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

The puzzle:

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.

Answer 1

I've run your coordinates through gdaltransform:

$ gdaltransform -s_srs EPSG:32017 -t_srs EPSG:4326
759232.003438, 1149854.52147    
-77.6116223688997 43.1517747887723 0

And it appears to come up with the right answer. This means that proj4 (which GDAL and PyProj are based on) is doing the right thing.

Sometimes these sorts of errors can be caused by having the coordinates swapped, but trying that with your sample point doesn't yield the wrong location you're getting. So the only thing I can think of is that PyProj is having trouble with a US-foot based projection. Try multiplying your input coordinates by 0.3048006096012192 which is the +to_meter scale value in the Proj4 definition.