Coordinate reference system transformation

Question

I am looking to transform the CRS of the geometry from my initial shapefile data to standard lat-long coordinates - espg 4236.

Here's the geopandas dataframe:

I wrote some script below that transforms the geometry to the espg CRS.

geodata.crs 

{'proj': 'lcc',
'lat_1': 37.06666666666667,
'lat_2': 38.43333333333333,
'lat_0': 36.5,
'lon_0': -120.5,
'x_0': 2000000,
'y_0': 500000.0000000002,
'ellps': 'GRS80',
'units': 'us-ft',
'no_defs': True}

pointsEB is point datatype: 

POINT (6091769.359598021 2090180.343520682) 

# Function to transform CRS projection 
inProj = Proj(geodata.crs)
outProj = Proj(init='epsg:4326')

pointsEBt = []

for pt in pointsEB: 

    x,y = transform(inProj, outProj, pt.x, pt.y)
    pointsEBt.append(Point(x,y))

Transformed Point: 
pointsEBt[1].x, pointsEBt[1].y
(-68.42841348959915, 40.16978821064519) --> Located in Antartica. 

My question: There seems to be something going on with the conversion to lat-long CRS. I used this as reference and "+init=epsg:4326" seems to be the format. However, when I plot a transformed point on the map the location is not what I expect it to be, far from it.

Answer 1

Converting the coordinates in your geodata dataframe seems to work correctly if I use the geopandas to_crs method:

In [20]: import geopandas

In [21]: from shapely.geometry import Point


In [26]: geodata = geopandas.GeoDataFrame(
...          {'OBJECTID': [1, 2, 3], 
...           'geometry': [Point(6065941, 2104148), Point(6143913, 2209458), Point(5879149, 2203020)]},
...          crs={'proj': 'lcc', 'lat_1': 37.06666666666667, 'lat_2': 38.43333333333333, 'lat_0': 36.5, 'lon_0': -120.5, 'x_0': 2000000, 'y_0': 500000.0000000002, 'ellps': 'GRS80', 'units': 'us-ft', 'no_defs': True})

In [27]: geodata
Out[27]: 
   OBJECTID                 geometry
0         1  POINT (6065941 2104148)
1         2  POINT (6143913 2209458)
2         3  POINT (5879149 2203020)


In [29]: geodata.to_crs(epsg=4326)
Out[29]: 
   OBJECTID                                      geometry
0         1  POINT (-122.2149623348043 37.76112770415131)
1         2  POINT (-121.9509132153954 38.05391734121174)
2         3  POINT (-122.8695903725446 38.02140044156164)

Those coordinates are located somewhere near San Francisco.

---

The reason that your above script did not work was due to the unit not being respected (see also answer of @user30184). If I add preserve_units=True, it does give the correct result:

In [44]: import pyproj

In [45]: inProj = pyproj.Proj(geodata.crs, preserve_units=True)
    ...: outProj = pyproj.Proj(init='epsg:4326')
    ...: 

In [46]: pt = geodata.geometry[0]

In [47]: x,y = pyproj.transform(inProj, outProj, pt.x, pt.y)

In [48]: x, y
Out[48]: (-122.2149623348043, 37.76112770415131)

But still, I would personally recommend to use the geopandas to_crs method, which already does this correctly and is also easier to use.

Answer 2

I inserted your CRS and coordinatates directly into gdaltransform https://www.gdal.org/gdaltransform.html

gdaltransform -t_srs epsg:4326 -s_srs "+proj=lcc +lat_1=37.06666666666667
 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000 +y_0=500000.0000000002 +ellps=GRS80
 +units=us-ft +no_defs"
6091769.359598021 2090180.343520682
-122.124788765566 37.7240393183059 0

But if I remove the units=us-ftfrom the proj4 string the result is

6091769.359598021 2090180.343520682
-69.7789972302181 40.583399518226 0

The latter result seems to be the same that you get. BTW. as stated that is not at Antarctica because for Proj.4 the first coordinate is always longitude of easting and the second one latitude or northing.

Maybe in your sript the proj4 string in geodata.crs gets truncated too early. Check the syntax.