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()
    #Transform incoming longitude/latitude to the shapefile's projection

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

    # check if point exists within shape
    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.scatter(ypts[truth_arr], xpts[truth_arr], c='green', s=2)
    plt.scatter(ypts[~truth_arr], xpts[~truth_arr], c='blue', s=2)
    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:

generated scatter plot

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:

region of interest


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


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
  • Hi ambal, very nice solution but I am not sure about the correctness of this, for your solution I am getting 20% landmass but according to mine, it should be at least around 50% (I am calculating 60% landmass for my estimated solution, see the map image). Could there be something wrong with the bounding box for your solution? – Varun Govind Dec 30 '20 at 1:19
  • 1
    Possibly, but it's also possibly that your land polygon is incomplete or the projection is off. Try printing the shapely output in an iPython window: print(land_poly) and print(bbox_poly.intersection(land_poly)). That should give you an indication of what is wrong. – amball Dec 30 '20 at 1:39
  • So I think the mistake is the buffer envelope line, because the bbox is in latitude-longitude coordinates (instead of kilometers). This makes it a little difficult to get the correct intersection but I think I have gotten it to work. I kind of combined the solution from John and yours and it gives me an accurate reading. I will add what I did in another answer. – Varun Govind Dec 30 '20 at 1:54
  • OK, it should use whatever CRS your shapefile is in (I thought that was already projected). But I'm glad you got it working. – amball Dec 30 '20 at 4:14

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):

solution comparison one

Here is the speedup I got:

speedup comparison

Not the answer you're looking for? Browse other questions tagged or ask your own question.