Calculating proportion of area that is land within bounding box using Python [duplicate]

Question

I am trying to find the proportion of land in an area specified within a bounding box (in python). The bounding box is given by min/max latitude & longitude coordinates.

The initial solution I came up with was to sample N random points within the bounding box and check (for each point) whether it lies within a landmass or not. Specifically, I checked the point within a landmass shapefile (from here: https://www.naturalearthdata.com/downloads/10m-physical-vectors/10m-land/). It doesn't give an exact area answer but an estimate.

This solution works but it is really slow (even with multithreading) so I was wondering if there is a better way to do this? I was thinking that you could isolate a landmass polygon within the bounding box from the shapefile and calculate the area within that. But I am very new to the osgeo package and shapefiles in general so I am not sure where to start with this problem.

Is there another way that I could do this?

Here is my reference code and some pictures of my initial solution (some code is missing like the bounding box calculation code, I can add it at request if needed):

import numpy as np
from osgeo import ogr
from itertools import repeat
import multiprocessing as mp
import matplotlib.pyplot as plt

def is_pt_land(lat, lon):
    shape = ogr.Open('./ne_10m_land_shape/ne_10m_land.shp', 0)
    shape_layer = shape.GetLayer()
    geo_ref = shape_layer.GetSpatialRef()
    point_ref = ogr.osr.SpatialReference()
    point_ref.ImportFromEPSG(4326)
    ctran=ogr.osr.CoordinateTransformation(point_ref,geo_ref)
    #Transform incoming longitude/latitude to the shapefile's projection
    [lon,lat,z]=ctran.TransformPoint(lon,lat)

    #Create a point
    pt = ogr.Geometry(ogr.wkbPoint)
    pt.SetPoint_2D(0, lon, lat)

    # check if point exists within shape
    shape_layer.SetSpatialFilter(pt) 
    return len(shape_layer)>0

def is_land_thresh(bounds, iters=5000, thresh=0.80):
    xpts = np.random.random_sample((iters,))
    ypts = np.random.random_sample((iters,))
    xpts = (abs(bounds[0][0]-bounds[1][0])*xpts) + bounds[0][0]
    ypts = (abs(bounds[0][1]-bounds[1][1])*ypts) + bounds[0][1]

    with mp.Pool() as pool: # applying the is land function to all the random points
        truth_arr = pool.starmap(is_pt_land, zip(xpts, ypts))

    truth_arr = np.asarray(truth_arr) # creating a visual of the truth array to generate an estimate
    plt.figure()
    plt.scatter(ypts[truth_arr], xpts[truth_arr], c='green', s=2)
    plt.scatter(ypts[~truth_arr], xpts[~truth_arr], c='blue', s=2)
    plt.show()
    return (np.sum(truth_arr)/iters) >= thresh
        
# Testing bbox and land threshold calc
loc = GeoLocation.from_degrees(26.062951, -80.238853) # create center point
bbox = loc.get_bbox(50)   # get bounding box with 50km from center point (100km x 100km roi)
is_land_thresh(bbox)      # determine if area is mostly landmass

Here is a photo of the scatter plot calculation:

Here is the landmass bounding box denoted by the min lat/long as the SW coord., the max lat/long as the NE coord., and the center coordinate:

Answer 1

Please see this answer, I think it is what your are trying to do.

Returning percentage of area of polygon intersecting another polygon using shapely

Answer 2

The simplest way is to calculate the area of the intersection between your bounding box polygon and your land polygon. Shapely and geopandas are a little more user friendly, but you could do it in ogr too.

import geopandas as gpd
import pyproj
from shapely.geometry import Polygon, Point
from shapely.ops import transform

shapeDf = gpd.read_file('./ne_10m_land_shape/ne_10m_land.shp')
land_poly = shapeDf.unary_union        # dissolve to single land polygon

crs = shapeDf.crs
project = pyproj.Transformer.from_crs(pyproj.CRS('EPSG:4326'), crs)
loc = transform(project.transform, Point(26.062951, -80.238853))

bbox_poly = loc.buffer(50).envelope    # create bounding box of circular buffer

bbox_area = bbox_poly.area
land_area = (bbox_poly.intersection(land_poly)).area   # land area within bbox
pc_land = land_area / bbox_area * 100
print(pc_land)

Answer 3

I found a solution based on John's link and Amball's answer. I think Amball's answer is correct if the projection is different.

The combined solution I used in code was this:

import fiona
import pyproj
from shapely.geometry import Polygon, Point, shape
from shapely.ops import transform, unary_union, cascaded_union    

# the new function
def is_land(bbox, shp_file, landmass, landmass_ratio=0.80):
    project = pyproj.Transformer.from_crs(pyproj.CRS('EPSG:4326'), shp_file.crs)
    # transform to coordinate system of shape file
    SW = transform(project.transform, Point(bbox[0][0], bbox[0][1]))
    NW = transform(project.transform, Point(bbox[0][0], bbox[1][1]))
    NE = transform(project.transform, Point(bbox[1][0], bbox[1][1]))
    SE = transform(project.transform, Point(bbox[1][0], bbox[0][1]))

    bbox_poly = Polygon([SW, NW, NE, SE, SW])

    area_intersection = (landmass.intersection(bbox_poly)).area/bbox_poly.area
    print('Landmass to ocean ratio: ', area_intersection)
    return area_intersection>=landmass_ratio

# Load shape file and combine into one landmass
shp_file = fiona.open('./ne_10m_land_shape/ne_10m_land.shp')
landmass_geom = [shape(feat["geometry"]) for feat in shp_file]
landmass_geom = unary_union(landmass_geom)

# create lat/lon bounding box using standard units (KM)
center = GeoLocation.from_degrees(26.062951, -80.238853)
bbox = center.get_bbox(50) #bbox in latitude and long coordinates with 50km from center location

# function call
is_land(bbox, shp_file, landmass_geom, landmass_ratio=0.80)

Here is the new solution compared with the old one in terms of accuracy (It is very close):

Here is the speedup I got: