2

I'm trying to calculate the area of a flood with Python, but the results seem way too small:

iceye_breaks = iceye_breaks.to_crs("+proj=cea +lon_0=0 +x_0=0 +y_0=0 +lat_ts=45 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

iceye_breaks['area'] = iceye_breaks['geometry'].area

iceye_breaks

These are the results:

enter image description here

The area such be much bigger than that. Could it be that geopandas doesn't calculate MultiPolygon areas correctly? Based on what I've read it should be able to calculate MultiPolygons.

I'd like the result in km2.

I've tried changing the projection and still get the same results. I can't use local epsg's since the same script needs to work worldwide, small errors aren't a worry since it's just supposed to be an estimation.

2
  • @BERA They have been projected to Cylindrical equal-area ("proj=cea")
    – UniStudent
    Apr 23 at 15:10
  • 1
    Is your original data in a geographic CRS (ie. lat/lons)? If so instead of converting to a projected CRS you can use methods outlined here gis.stackexchange.com/q/413349
    – Shawn
    Apr 23 at 15:41

1 Answer 1

3

You can create a function to find the utm zone for the centroid of the lat, long point then project to that zone and calculate area. I think it should work worldwide

import geopandas as gpd
import utm #pip install utm
from pyproj import CRS

df = gpd.read_file(r'/home/bera/Desktop/GIStest/4_buildings_4326.shp')

#  df.geometry.area
#0    1.499222e-06
#1    1.783689e-06
#2    7.173407e-07
#3    1.983493e-06


def findtheutm(aGeometry):
"""A function to find a coordinates UTM zone"""
    x, y, parallell, latband = utm.from_latlon(aGeometry.centroid.y, aGeometry.centroid.x)
    if latband in 'CDEFGHJKLM': #https://www.lantmateriet.se/contentassets/379fe00e09d74fa68550f4154350b047/utm-zoner.gif
        ns = 'S'
    else:
        ns = 'N'
    crs = "+proj=utm +zone={0} +{1}".format(parallell, ns) #https://gis.stackexchange.com/questions/365584/convert-utm-zone-into-epsg-code
    crs = CRS.from_string(crs)
    _, code = crs.to_authority()
    return int(code)

epsg = findtheutm(df.geometry.iloc[0]) #Input the dfs first geometry to the function and get utm epsg code back

df['correctArea'] = df.to_crs(epsg).area # /1e6 for km2

# df
#         id  ...   correctArea
# 0   585331  ...   9140.157424
# 1   585332  ...  10874.176357
# 2   585337  ...   4373.700636
# 3  2488763  ...  12093.625510
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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