7

I have some points in EPSG:2284 and I'm trying to convert them to EPSG:4326 using pyproj, however I get weird results. When I go to the spatial reference website I see that there is an endpoint that does the conversion. When I use that instead I get the correct result.

Here is some code illustrating the different results:

from pyproj import Proj, transform
import requests

if __name__ == '__main__':

    # points I want to convert from EPSG:2284 -> EPSG:4326
    targets = [
        (3669486.63386,11293013.10592,709.08004),
        (3669559.13811,11292972.72831,711.60055),
        (3669639.51308,11292851.48264,712.22258),
        (3669800.62304,11292949.38118,714.75766),
    ]

    for northing, easting, up in targets:
        print "Converting", easting, northing

        # Transform using pyproj (gives wrong answer)
        WGS84 = Proj(init='EPSG:4326')
        inp = Proj(init='EPSG:2284')
        x, y = transform(inp, WGS84, easting, northing)
        print x, y, "(with proj4)"

        # Transform using spatialreference.org URL (gives correct answer)
        url = "http://spatialreference.org/projection/?json=%7B%22type%22%3A%22Feature%22%2C%20%22geometry%22%3A%7B%22type%22%3A%22Point%22%2C%20%22coordinates%22%3A%5B{lon}%2C%20{lat}%5D%7D%2C%22properties%22%3A%7B%7D%7D&inref=EPSG:{inp}&outref=epsg:{outp}"
        url = url.format(lon=easting, lat=northing, inp='2284', outp='4326')
        x, y = requests.get(url).json()['coordinates']
        print x, y, '(using spatialreference.org URL)'

Does anyone know why they are giving different result and what is wrong with the Pyproj code?

til_b
  • 5,074
  • 1
  • 19
  • 36
nickponline
  • 825
  • 1
  • 12
  • 22

2 Answers2

9

Simply use (Converting elevations into correct units with pyproj?, Proj4 String for NAD83(2011) / Louisiana South (ftUS), ...)

preserve_units=True (as you say, pyproj assumes that your coordinates are in meters, therefore)

for northing, easting, up in targets:
    print "Converting", easting, northing
    # Transform using pyproj (gives wrong answer)
    WGS84 = Proj(init='EPSG:4326')
    inp = Proj(init='EPSG:2284',preserve_units=True)
    ...

Results

Converting 11293013.1059 3669486.63386
-79.1537547735 37.3989918968 (with proj4)
-79.153755 37.398992 (using spatialreference.org URL)
Converting 11292972.7283 3669559.13811
-79.1538955003 37.3991902547 (with proj4)
-79.153896 37.39919 (using spatialreference.org URL)
Converting 11292851.4826 3669639.51308
-79.1543148015 37.3994086895 (with proj4)
-79.154315 37.399409 (using spatialreference.org URL)
Converting 11292949.3812 3669800.62304
-79.1539816309 37.3998530254 (with proj4)
-79.153982 37.399853 (using spatialreference.org URL)

Supplementary control with GDAL/OSR

from osgeo import osr
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
inp= osr.SpatialReference()
inp.ImportFromEPSG(2284)
transformation = osr.CoordinateTransformation(inp,wgs84)
for northing, easting, up in targets:
    print transformation.TransformPoint(easting,northing)

(-79.15375477347949, 37.3989918968496, 0.0)
(-79.15389550031416, 37.39919025468435, 0.0)
(-79.15431480154861, 37.39940868945351, 0.0)
(-79.15398163085317, 37.39985302535012, 0.0)
gene
  • 54,868
  • 3
  • 110
  • 187
1

Pyproj expects meters not feet.

nickponline
  • 825
  • 1
  • 12
  • 22