How to check lon/lat polygon pixels over land or ocean quickly?

Question

I have 2d lon/lat arrays and am trying to check the land type like this:

import numpy as np
from shapely.geometry import Polygon
import cartopy.io.shapereader as shpreader
from shapely.ops import unary_union

lon = np.arange(-180, 181, .1)
lat = np.arange(-90, 91, .1) 
lons, lats = np.meshgrid(lon, lat)

land_shp_fname = shpreader.natural_earth(resolution='110m',  
                                         category='physical', name='land')

land_geom = unary_union(list(shpreader.Reader(land_shp_fname).geometries()))

grid_names = np.empty_like(lons, dtype=int)

for i in range(len(lon)-1):
    for j in range(len(lat)-1):
        poly = Polygon([(lon[i], lat[j]), (lon[i+1], lat[j]),  
                        (lon[i+1], lat[j+1]), (lon[i], lat[j+1])])
        if poly.intersects(land_geom):
            grid_names[j,i] = 1 # Land
        else:
            grid_names[j,i] = 0 # Ocean

The speed is slow for creating the high-resolution one for 1000x1000 pixels. Any suggestions for improvement?

Answer 1

I find the roaring_landmask package is really fast:

import numpy as np
from roaring_landmask import RoaringLandmask

lon = np.arange(-180, 180, .1)
lat = np.arange(-90, 90, .1) 
lons, lats = np.meshgrid(lon, lat)

l = RoaringLandmask.new()
mask = l.contains_many(lons.ravel(), lats.ravel())

It costs around 5s on my laptop ;)

Update

I find it's better to use geopandas and cartopy do this, because roaring_landmask doesn't consider lake as not_land type.

Here's an example of checking points. It should be similar for polygons.

import numpy as np
import geopandas as gpd
from geopandas import GeoSeries
import cartopy.feature as cfeature

def random_lon_lat(n=1, lat_min=-90., lat_max=90., lon_min=-180., lon_max=180.):
    """
    this code produces an array with pairs lat, lon
    """
    lat = np.random.uniform(lat_min, lat_max, n)
    lon = np.random.uniform(lon_min, lon_max, n)

    points = GeoSeries(gpd.points_from_xy(lon, lat))
    points_gdf = gpd.GeoDataFrame(geometry=points, crs="EPSG:4326")

    return points_gdf

points_gdf = random_lon_lat(int(1e6))

land_50m = cfeature.NaturalEarthFeature('physical','land','50m')
land_polygons_cartopy = list(land_50m.geometries())

land_gdf = gpd.GeoDataFrame(crs='epsg:4326', geometry=land_polygons_cartopy)

# Spatially join the points with the land polygons
joined = gpd.sjoin(points_gdf, land_gdf, how='left', op='within')

# Check if each point is within a land polygon
is_within_land = joined['index_right'].notnull()

It takes less than 1 second to check all point types.