3

I'm trying to convert Projected bounds to WGS84 using Pyproj but getting strange results.

The Projected Bounds are : 909126.0155, 110626.2880, 1610215.3590, 424498.0529

The corresponding WGS84 Bounds (the target) are: -74.2700, 40.4700, -71.7500, 41.3100

The data is from: http://spatialreference.org/ref/epsg/nad83-new-york-long-island-ftus/

>>> from pyproj import Proj, transform
>>> inProj  = Proj("+init=EPSG:2236")
>>> outProj = Proj("+init=EPSG:3857")
>>> x1,y1 =  (110626.2880, 909126.0155)
>>> print(transform(inProj,outProj,x1,y1))
(-9122788.926222347, 3833518.013375292)

The outProj is based on my understanding of what Google Maps uses.

Improve this question
6
  • Try "2263" instead.
    – mkennedy
    Apr 20 '18 at 22:44
  • Did for inProj and outProj. Not getting: -74.2700, 40.4700
    – zadrozny
    Apr 20 '18 at 22:48
  • If you want lat/lon values, use 4269 (NAD83) or 4326 (WGS84) for outProj.
    – mkennedy
    Apr 20 '18 at 22:52
  • Both outProj = Proj("+init=EPSG:4269") and outProj = Proj("+init=EPSG:4326") give me (-76.53059046753678, 48.30210692217332) That's much closer but still off.
    – zadrozny
    Apr 20 '18 at 22:59
  • EPSG:3857 is not a geographic CRS and doesn't use lon/lat coordinates, its horizontal units are metres (i.e. not degrees). You see lon/lat when you click or search in Google Maps as map coordinates are projected on the fly.
    – user2856
    Apr 20 '18 at 23:12
6

Look Converting EPSG:2284 to EPSG:4326 with pyproj (and many others...).

Pyproj expects degrees (lon, lat) or meters (x,y) as units but the unit of Projection: 2263 isUNIT["US survey foot"..., therefore you need to use preserve_units=True.

from pyproj import Proj, transform
inProj  = Proj("+init=EPSG:2263",preserve_units=True)
outProj = Proj("+init=EPSG:4326") # WGS84 in degrees and not EPSG:3857 in meters)
# swap x,y as mkennedy says
y1,x1 =  (110626.2880, 909126.0155)
print(transform(inProj,outProj,x1,y1))
(-74.2700000001129, 40.46999999990434)

Control with GDAL/OSR

from osgeo import osr
inp= osr.SpatialReference()
inp.ImportFromEPSG(2263)
out= osr.SpatialReference()
out.ImportFromEPSG(4326)
transformation = osr.CoordinateTransformation(inp,out)
print(transformation.TransformPoint(x1,y1))
(-74.27000000011289, 40.46999999990432, 0.0)

New : Now how do I reverse the result?

y1,x1 =  (110626.2880, 909126.0155) # and not x1,y1 =
lon,lat = transform(inProj,outProj,x1,y1)
print(lon, lat)
-74.2700000001129 40.46999999990434

Now reverse the result:

x, y = transform(outProj,inProj,lon,lat)
print(x,y)
909126.0155000015 110626.28799994898

And you get the original values if you use print(y,x)

print(y,x)
110626.28799994898 909126.0155000015
Improve this answer
3
  • Now how do I reverse the result? I have swapped inProj and outProj and experimented with all permutations of preserve_units and none gets me back to (110626.2880, 909126.0155). Sorry, I am new to GIS and lack a conceptual grasp.
    – zadrozny
    Apr 21 '18 at 17:17
  • If you haven't tried this yet, also try swapping the latitude and longitude values.
    – mkennedy
    Apr 21 '18 at 21:43
  • Did you understand the answer correctly (# swap x,y as mkennedy says) ? Look at New in my answer
    – gene
    Apr 22 '18 at 10:13

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.