Inverse X Y to Lat Lon, data not as expected

[Jestre]

Hello,

I am trying to convert a script I wrote in Perl years ago to Python, and while I believe I've followed the directions, the results do not match up with what I'd expect.

Code:

from pyproj import Proj
nc = Proj(init='epsg:3359') # I've tried most of the available NC options (site below)
y = [2110926.0]
x = [0738252.0]
lon, lat = nc(x, y, inverse=True)
print lon, lat

[-77.19861281130545] [52.46179322191929]

This address should return as -78.7122, 35.7684 or thereabouts.

Am I missing some step? I've also tried entering all of the parameters individually with the same results.

http://spatialreference.org/ref/?search=carolina

Thanks

[micahcochran]

Are the coordinates in feet or meters? Pyproj will assume that you want meters so you might try this nc = Proj(init='epsg:3359', preserve_units=True). BTW that didn't exactly produce your target coordinates.

Do you have a proj4 string, WKT, shapefile, or something else that defines the projection?

[Jestre]

Based on the setup I had in the Perl script (below), I am assuming they are in ft. These are data exported from a 911 CAD system, so I don't have access to a better definition than that.

my $proj = Geo::Proj4->new(
proj => "lcc",
units => "us-ft",
datum => "NAD83",
lat_0 => 33.75,
lon_0 => -79,
x_0 => 609601.22,
y_0 => 0,
lat_1 => 34.33,
lat_2 => 36.166667
)

Based on that config, the values match EPSG reference is 32019 (which I'm thinking is working in feet also based on the to_meters option.

http://spatialreference.org/ref/epsg/32019/proj4/

I will try the preserve units.

[micahcochran]

Hold on epsg 32019 uses the older NAD27 datum, and the above uses NAD83. Changing those can cause things willy nilly can make your data be off by hundreds of feet.

>>> from pyproj import Proj
>>> nc = Proj("+proj=lcc +lat_1=34.33333333333334 +lat_2=36.16666666666666 +lat_0=33.75 +lon_0=-79 +x_0=609601.2192024384 +y_0=0 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs ", preserve_units=True)
>>> nc(2110926, 738252, inverse=True)
(-78.62602489780987, 35.77787470941835)

That seems to work. I used the one that you defined in EPSG:32019 and combined it with the above. Or translate the one used in the code. I would encourage you to use a standard EPSG number or state plane projection if at all possible.

[jswhit]

Here's what I get with github master:

from pyproj import Proj
p = Proj(init='epsg:3359')
lon, lat = p(21109260., 7382520., inverse=True)
print lon, lat

68.1559778764 -40.6009157221

[Jestre]

Thank you all for your assistance. It appears that ESRI:102719 was what I was looking for. It differs minimally in the +x_0 setting from that posted above, but the answers are certainly within the limits allowed by this map.