2

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?

3

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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