Getting inaccurate results converting from New York State projection to NAD83 with Python's pyproj

Question

I've been trying to do a standard coordinate translation between New York State Long Island (ESRI:102718) and latitude/longitude (EPSG:4269)...but I keep getting the wrong lat/lng results.

I don't think this has anything to do with the pyproj library in particular, but maybe I'm missing something in the transformation process.

Sample data comes from NYC.gov: http://www1.nyc.gov/site/planning/data-maps/open-data/dwn-selfac.page

Its documentation says:

The X coordinate of the XY coordinate pair which depicts the approximate location of the facility, program, or parkland property. (Also see Y COORDINATE). Facilities' X and Y coordinates must be used together for mapping selection. The XY coordinates are expressed in the New York Long Island State Plane coordinate system.

Here's a sample record:

```

Sample record:

"3861130100",1,"MN",104,"MANHATTAN COMMUNITY DISTRICT 4",731,1,"YAI - YOUNG ADULT INSTITUTE", "460 W 34 St","10001",3807,7,72,,98,,,,91,57," 103","0103",3,2,10,5500,1007310001.000, 984533.000,213998.000,"L" ```

So, the python code:

import pyproj NYSP1983 = pyproj.Proj(init="ESRI:102718") x, y = (984533, 213998) NYSP1983(x, y, inverse=True)

That inverse projection should result in the lat/lng coordinates. However, the result is nowhere near New York:

       (-65.75200831429146, 41.80318458470038)

According to Google geocoder, these are the proper coordinates:

    location: {
        lat: 40.7541045,
        lng: -73.9989178
    }

Here's another way to do transformations in pyproj, going from latlng to the specified X/Y projection, but still no dice:

import pyproj NYSP1983 = pyproj.Proj(init="ESRI:102718") lng, lat = (-73.9989178, 40.7541045) NYSP1983(lng, lat)

The result:

            (300091.3883012832, 65232.600609979076)

Doing an explicit transform between systems doesn't help either:

import pyproj NAD83 = pyproj.Proj(init="EPSG:4269") NYSP1983 = pyproj.Proj(init="ESRI:102718") x, y = (984533, 213998) pyproj.transform(NYSP1983, NAD83, x, y)

         (-65.75200831429146, 41.80318458470038)

What am I missing here?

Edit

Seems that I'm missing some parameters that need to be set up, such as False_Easting...but I don't know how in pyproj (or in any system) how to specify that. Namely, how do I use what's in this PROJCS string (example here: https://gist.github.com/veltman/6be84d216376814cdc67) to initialize the desired projection?

Answer 1

The problem was my reading of the documentation of the pyproj.Proj() constructor function. The preserve_units optional argument has to be set to True:

https://github.com/jswhit/pyproj/issues/40#issuecomment-165585816

NYSP1983 = pyproj.Proj(init="ESRI:102718", preserve_units=True)
x, y = (984533, 213998)
NYSP1983(x, y, inverse=True)

Output

   (-73.99897854592253, 40.75405155302713)

It turns out pyproj.Proj() will ignore the +to_meter param without the preserve_units=True being included as an argument. i.e., you have to do preserve_units=True in this situation:

projstr = '+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs'

myproj = pyproj.Proj(projstr, preserve_units=True)
myproj(x, y, inverse=True)