4

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:

enter image description here

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.

4
  • 2
    Keep in mind that, counter-intuitively, Latitude should be plotted on the Y-axis and Longitude on the X-axis. This is because these values ascend and descend on a perpendicular axis to the direction they traverse on the globe. Try swapping x for y and see if it produces your expected result? Aug 31, 2018 at 0:58
  • What is geodata.crs? It's not EPSG:4326. How are you importing/reading this? Aug 31, 2018 at 0:59
  • @Michael The question has been updated with geodata.crs. I am reading the CRS system of the geometry and them attempting to transform it.
    – kms
    Aug 31, 2018 at 1:05
  • @RobinHorner Yeah, good point. however, I don't think that is the case here - I tried swapping x for y.
    – kms
    Aug 31, 2018 at 1:08

2 Answers 2

3

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.

1
  • Yep - adding preserve_units=True worked for me. According to the documentation when not true it forces the units to be in meters.
    – kms
    Aug 31, 2018 at 15:19
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.