diff --git a/GDE_TOOLS_world_grid.py b/GDE_TOOLS_world_grid.py index 233d35d3ab7f5aececfe4208862d07b35c56decd..86c0dacfd58544719b5fd7acd003b14bd90c18fe 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