Commit 67308692 by Cecilia Nievas

### Fixed bug in get_area_polygon_on_Earth()

parent 56afc9da
 ... ... @@ -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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!