From 673086927590c6b50cfd78c5df319e033a3835cf Mon Sep 17 00:00:00 2001
From: Cecilia Nievas <cnievas@gfz-potsdam.de>
Date: Mon, 8 Nov 2021 17:28:06 +0100
Subject: [PATCH] Fixed bug in get_area_polygon_on_Earth()

---
 GDE_TOOLS_world_grid.py | 52 +++++++++++++++++------------------------
 1 file changed, 21 insertions(+), 31 deletions(-)

diff --git a/GDE_TOOLS_world_grid.py b/GDE_TOOLS_world_grid.py
index 233d35d..86c0dac 100644
--- a/GDE_TOOLS_world_grid.py
+++ b/GDE_TOOLS_world_grid.py
@@ -349,48 +349,38 @@ def get_cells_in_bbox_as_pandas_with_WKT(sw_lon, sw_lat, ne_lon, ne_lat, zoomlev
     return out_df    
     
 
-def get_area_polygon_on_Earth(polygon_in,units='m2'):
-    """This function returns the area in m2 or km2 of a polygon on Earth whose geometry is
-    defined in epsg:4326 ("normal" latitude and longitude).
-
-    This function was written based on the code in:
-    https://gist.github.com/robinkraft/c6de2f988c9d3f01af3c
-    (general process)
-    and
-    https://stackoverflow.com/questions/51554602/
-    how-do-i-get-the-area-of-a-geojson-polygon-with-python
-    (improved projection, instead of epsg:3857 as suggested above)
-    and
-    https://stackoverflow.com/questions/4681737/
-    how-to-calculate-the-area-of-a-polygon-on-the-earths-surface-using-python
-    (adding parameters apart from lat1 and lat2 to partial(... pyproj.Proj), which are used to
-    define the Albers Equal Area projection)
+def get_area_polygon_on_Earth(input_polygon, units='m2'):
+    """This function returns the area (in m2 or km2) of a polygon on Earth calculated using the
+    Albers equal area projection. The polygon's geometry is assumed to be defined in ESPG:4326.
 
     Args:
-        polygon_in (Shapely polygon): Polygon defined in epsg:4326.
+        input_polygon (Shapely polygon): Polygon defined in epsg:4326.
         units (str): The units of the output, can be m2 or km2 (default: m2).
         
     Returns:
         projected_area: Area of the polygon in m2 or km2.
     """
+    
+    # Get the bounding box of the polygon
+    lon_w = input_polygon.bounds[0]
+    lat_s = input_polygon.bounds[1]
+    lon_e = input_polygon.bounds[2]
+    lat_n = input_polygon.bounds[3]
+
+    # Use the resulting coordinates to carry out the transformation
+    project_albers = pyproj.Proj(
+        "+proj=aea +lat_1={} +lat_2={} +lat_0={} +lon_0={}".format(
+            lat_s, lat_n, (lat_s + lat_n) / 2.0, (lon_w + lon_e) / 2.0
+        )
+    )
+    geometry = shapely.ops.transform(project_albers, input_polygon)
+
+    projected_area = geometry.area  # in m2
 
-    bbox = polygon_in.bounds # sw_lon, sw_lat, ne_lon, ne_lat
-    # Create projection from epsg:4326' to Albers Equal Area
-    # (https://proj.org/operations/projections/aea.html)
-    # For this it is necessary to provide as input the area of the world of interest, which is
-    # input herein from the bounding box of the polygon whose area we want to calculate.
-    inputstr = '+proj=aea +lat_1={} +lat_2={} +lat_0={} +lon_0={}'.format(bbox[1],
-                                                                          bbox[3],
-                                                                          (bbox[1]+bbox[3])/2.0,
-                                                                          (bbox[0]+bbox[2])/2.0)
-    proj = partial(pyproj.transform,
-                   pyproj.Proj(init='epsg:4326'),
-                   pyproj.Proj(inputstr))
-    projected_area = shapely.ops.transform(proj, polygon_in).area # in m2
     if units == 'km2':
         projected_area = projected_area / 1000000.0
     elif units != 'm2':
         print('      UNITS NOT RECOGNISED IN get_area_km2_polygon_on_Earth FUNCTION. ERROR!!!')
         projected_area= -999.9
 
-    return projected_area 
+    return projected_area
-- 
GitLab