Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
add comment

1 Answer 1

up vote 5 down vote accepted

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.

share|improve this answer
1  
+1 The foot vs. meter difference would neatly explain the 1000 km discrepancy. (I compute a 975.9 km error in projected coordinates.) –  whuber Nov 29 '11 at 17:27
    
Yes, in question gis.stackexchange.com/questions/10209, @geographika said much the same thing. –  MerseyViking Nov 29 '11 at 17:36
    
Thanks. Yes, it was a foot vs. meter problem. And yes, I found gis.stackexchange.com/questions/10209 after asking my question and it is basically the same question. The answer was in code.google.com/p/pyproj/issues/detail?id=20 (i.e., pass preserve_units=True as a second parameter when calling the pyproj.Proj() constructor). I find it very strange that the default behavior of pyproj is to silently work in meters, even when the projections involved are in feet. From what I can tell, using preserve_units=True means "do it the right way" instead of "sometimes do it the wrong way." –  M Katz Nov 29 '11 at 18:40
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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