4

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.

zadrozny
  • 298
  • 1
  • 4
  • 14
  • 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
  • Oh, try swapping x1,y1 values. When I use x = 909126 and y = 110626, in Esri projection engine, I get -74.27 and 40.47 (more or less). Forgive me if you know this, but US coordinates are usually x = easting, y = northing. – mkennedy Apr 20 '18 at 23:19

1 Answers1

7

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
gene
  • 54,868
  • 3
  • 110
  • 187
  • 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