I'm working with a geodataset of fields. I want to compute the surface of each parcel on the fly when they are displayed on the map (the rational for that is not relevant here).
When I was able to manipulate the dataset, I would have simply do the following:
# I know 3857 is not optimal everywhere but that will do for the example
gdf["surface"] = gdf.to_crs(3857).area
Now what I have is the following:
feat = gdf[gdf.id == change["new"]].squeeze()
# this is a shaely geometry in 4326 and before computing area,
# I would need to reproject it
surface = feat.geometry
from shapely.geometry import Polygon
from shapely.ops import transform
import pyproj
poly = Polygon([[20,63], [20.001,63], [20.001,63.001], [20, 63.001], [20,63]])
wgs84 = pyproj.CRS('EPSG:4326')
mercator = pyproj.CRS('EPSG:3857')
project = pyproj.Transformer.from_crs(wgs84, mercator, always_xy=True).transform
poly_3857 = transform(project, poly)
print(poly_3857.area)
27296.256818790367