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-bounds)*xpts) + bounds ypts = (abs(bounds-bounds)*ypts) + bounds 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: