Commit cb1bcdc9 authored by Cecilia Nievas's avatar Cecilia Nievas
Browse files

Merge branch 'bugfix/area_function' into 'master'

Fixed bug in get_area_polygon_on_Earth()

Closes #12

See merge request !21
parents 56afc9da 67308692
...@@ -349,48 +349,38 @@ def get_cells_in_bbox_as_pandas_with_WKT(sw_lon, sw_lat, ne_lon, ne_lat, zoomlev ...@@ -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 return out_df
def get_area_polygon_on_Earth(polygon_in,units='m2'): def get_area_polygon_on_Earth(input_polygon, units='m2'):
"""This function returns the area in m2 or km2 of a polygon on Earth whose geometry is """This function returns the area (in m2 or km2) of a polygon on Earth calculated using the
defined in epsg:4326 ("normal" latitude and longitude). Albers equal area projection. The polygon's geometry is assumed to be defined in ESPG:4326.
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)
Args: 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). units (str): The units of the output, can be m2 or km2 (default: m2).
Returns: Returns:
projected_area: Area of the polygon in m2 or km2. 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': if units == 'km2':
projected_area = projected_area / 1000000.0 projected_area = projected_area / 1000000.0
elif units != 'm2': elif units != 'm2':
print(' UNITS NOT RECOGNISED IN get_area_km2_polygon_on_Earth FUNCTION. ERROR!!!') print(' UNITS NOT RECOGNISED IN get_area_km2_polygon_on_Earth FUNCTION. ERROR!!!')
projected_area= -999.9 projected_area= -999.9
return projected_area return projected_area
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment