Transforming EPSG:26918 to EPSG:4236 in GeoPandas

Question

I am trying to transform a shapefile of polygons projected on CRS: EPSG:26918 to EPSG:4326 so that I can perform a spatial join in GeoPandas over two shapefiles.

My previous code allowed me to simply set the CRS as:

gdf.crs = {'init': 'epsg:4326'}

However, with GeoPandas now using a new version of proj that reverses how it reads long,lat, my map is now completely reversed after that operation -- it is upside down and backwards.

I get the warning:

..\crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated.
'<authority>:<code>' is the preferred initialization method. 
When making the change, be mindful of axis order changes: 
https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 return _prepare_from_string(" ".join(pjargs))

Is there an easy way to handle this?

More info on gdf.crs before conversion:

<Projected CRS: EPSG:26918>
Name: NAD83 / UTM zone 18N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: North America - 78°W to 72°W and NAD83 by country
- bounds: (-78.0, 28.28, -72.0, 84.0)
Coordinate Operation:
- name: UTM zone 18N
- method: Transverse Mercator
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

CRS of shapefile I'd like to spatial join with:

<Geographic 2D CRS: +init=epsg:4326 +type=crs>
Name: WGS 84
Axis Info [ellipsoidal]:
- lon[east]: Longitude (degree)
- lat[north]: Latitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Edit: Here is the full operation:

Answer 1

Not sure what you were doing but the following do the reprojection without any issue (axis inversion). I've used the following sample https://labs.webgeodatavore.com/partage/nyc_census_blocks.zip (contains nyc_census_blocks.shp and related associated files that use EPSG 26918 like in your case)

import geopandas

# Read input shp
gdf_26918 = geopandas.read_file('nyc_census_blocks.shp')
# Reproject and assign to a new variable
gdf_4326 = gdf_26918.to_crs(epsg=4326)
# Export
gdf_4326.to_file("nyc_census_blocks_4326.shp")

An alternative could be to use pyproj CRS objet

import geopandas
from pyproj import CRS

# Read input shp
gdf_26918 = geopandas.read_file('nyc_census_blocks.shp')
# Explicitly prepare EPSG 4326 Pyproj CRS object
crs_4326 = CRS("epsg:4326")
# Reproject and assign to a new variable using Pyproj CRS object
gdf_4326 = gdf_26918.to_crs(crs=crs_4326)
# Export
gdf_4326.to_file("nyc_census_blocks_4326.shp")

Answer 2

With geopandas 0.6+ & pyproj 2.2+, the axis order shouldn't be an issue: https://github.com/geopandas/geopandas/pull/1122

So, you should be able to do:

gdf4326 = gdf.to_crs('EPSG:4326')

And it should work properly and the warning should go away.

Answer 3

First, you need to see if your layer's CRS is defined. You can check it by:

import geopandas as gpd    
gdf = gpd.read_file(<file path>)
print(gdf.crs)

If CRS is not defined, then you need to set the CRS:

gdf = gdf.set_crs(epsg=26918)

Then you can reproject it:

gdf = gdf.to_crs(epsg=4326)

Now you can perform your spatial join:

gdf2 = gpd.read_file(<your second layer>)    
result = gpd.sjoin(gdf, gdf2, how=<the join method>, op=<spatial relation>)

You need to check if the second layer has its CRS defined properly so you can get correct results without any warnings.