如何快速检查陆地或海洋上的经/纬度多边形像素？

Question

我有二维经度/纬度数组，正在尝试检查土地类型，如下所示：

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

创建 1000x1000 像素的高分辨率图像速度较慢。有什么改进建议吗？

Answer 1

我发现roaring_landmask包真的很快：

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

在我的笔记本电脑上花费大约 5 秒；）