[GIS] Converting EPSG:2284 to EPSG:4326 with pyproj

coordinate systemepsgpyprojpython

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?

Best Answer

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)
Related Question