geopandas reproject change the extent of the vector

Question

Geopandas reproject changes the extent of the vector. Below is an example, you can download the vector file from here

# First please download the shp file

import matplotlib.pyplot as plt
import geopandas
from shapely.geometry import Polygon

# Read the file using geopandas
domain = geopandas.read_file("data.shp")

#Create a polygone
polygon = Polygon([(-100,62),(-170,62),(-170,64),(-100,64)])
poly_gdf = geopandas.GeoDataFrame([1], geometry=[polygon], crs=domain.crs)

fig, ax1 = plt.subplots(1, 1, figsize=(12, 8))
domain.plot(ax=ax1)
poly_gdf.boundary.plot(ax = ax1,color="red")

Now let's reproject both and plot them again:

# Reproject both polygons to Albers Equal Area projection
domain_reproj = domain.to_crs('+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs')
poly_gdf_reproj = poly_gdf.to_crs('+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs')

# And plot them again
fig, ax1 = plt.subplots(1, 1, figsize=(12, 8))
domain_reproj.plot(ax=ax1)
poly_gdf_reproj.boundary.plot(ax = ax1, color = 'red')

As you can after reprojecting both polygons to the same projection system the overlapping area changes. The yellow circles help to spot the difference. How can I fix this?

Answer 1

Geopandas is just re-projecting each of the points in your polygons ie "(-100,62),(-170,62),(-170,64),(-100,64)" rather than projecting along the line.

You can check that poly_gdf.iloc[0].geometry.wkt only give you five points where as you would expect a wide curve made from lots of small line segments in the projected map.

If you want a polygon that will project nicely one option is to add more points along each edge of your polygon.

As Michael suggested, you can do this really cleanly with shapely's interpolate

original_polygon = Polygon([(-100,62),(-170,62),(-170,64),(-100,64),(-100,62)])

polygon = Polygon([original_polygon.boundary.interpolate(i, normalized=True) for i in np.linspace(0, 1, 1000)])

poly_gdf = geopandas.GeoDataFrame([1], geometry=[polygon], crs=domain.crs)