如何快速判断经纬度网格多边形像素属于陆地还是海洋？

原代码的核心性能瓶颈在于Python层面的双重循环逐格创建多边形并判断相交——对于高分辨率网格，需要执行数十万甚至数百万次几何运算，Python循环本身效率极低，且每次创建Polygon对象的开销也很大。以下是几种针对性的优化方案：

方案1：使用栅格化工具批量处理（最推荐） #

直接调用底层优化的栅格化函数，一次性将陆地多边形转换为栅格数据，速度比逐循环快几个数量级。推荐使用rasterio库的rasterize函数：

import numpy as np
import cartopy.io.shapereader as shpreader
from rasterio.features import rasterize
from affine import Affine

# 定义网格参数
lon_step = 0.1
lat_step = 0.1
lon = np.arange(-180, 181, lon_step)
lat = np.arange(-90, 91, lat_step)
# 栅格的实际尺寸（网格单元数量）
grid_height = len(lat) - 1
grid_width = len(lon) - 1

# 读取陆地矢量数据
land_shp_fname = shpreader.natural_earth(resolution='110m', category='physical', name='land')
land_geoms = list(shpreader.Reader(land_shp_fname).geometries())

# 定义栅格的仿射变换参数（关联地理坐标与栅格像素位置）
transform = Affine(lon_step, 0, lon[0],
                   0, -lat_step, lat[-1])  # 纬度从上到下递减，故y方向步长为负

# 栅格化陆地：1代表陆地，0代表海洋
grid_names = rasterize(
    land_geoms,
    out_shape=(grid_height, grid_width),
    transform=transform,
    fill=0,
    default_value=1,
    dtype=np.int32
)

优化原理 #

rasterize是基于C实现的批量栅格化算法，直接在底层处理所有多边形与栅格的映射，完全避开了Python循环的开销，处理1000×1000的网格仅需几秒。

方案2：预处理陆地几何，建立空间索引 #

如果不想引入新库，可以通过shapely.prepared模块为陆地几何建立空间索引，大幅加速相交判断：

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

lon = np.arange(-180, 181, 0.1)
lat = np.arange(-90, 91, 0.1) 

# 读取并合并陆地几何
land_shp_fname = shpreader.natural_earth(resolution='110m', category='physical', name='land')
land_geom = unary_union(list(shpreader.Reader(land_shp_fname).geometries()))

# 预处理陆地几何，建立空间索引
prepared_land = prep(land_geom)

grid_names = np.empty((len(lat)-1, len(lon)-1), dtype=int)

# 循环判断，但使用预处理后的几何加速相交检测
for j in range(len(lat)-1):
    for i in range(len(lon)-1):
        poly = Polygon([(lon[i], lat[j]), (lon[i+1], lat[j]),  
                        (lon[i+1], lat[j+1]), (lon[i], lat[j+1])])
        # 预处理后的intersects比原生方法快数倍
        grid_names[j,i] = 1 if prepared_land.intersects(poly) else 0

优化原理 #

prep函数会为几何对象构建空间索引，后续的相交判断会先通过索引过滤掉不可能相交的区域，减少不必要的几何计算，能将原代码的速度提升3-5倍。

方案3：并行化循环处理 #

利用Python的多进程库，将网格拆分为多个块并行处理，充分利用多核CPU资源：

import numpy as np
from shapely.geometry import Polygon
import cartopy.io.shapereader as shpreader
from shapely.ops import unary_union
from shapely.prepared import prep
from multiprocessing import Pool

def process_chunk(j_range):
    """处理指定行范围的网格单元"""
    chunk_result = np.empty((len(j_range), len(lon)-1), dtype=int)
    for idx, j in enumerate(j_range):
        for i in range(len(lon)-1):
            poly = Polygon([(lon[i], lat[j]), (lon[i+1], lat[j]),  
                            (lon[i+1], lat[j+1]), (lon[i], lat[j+1])])
            chunk_result[idx,i] = 1 if prepared_land.intersects(poly) else 0
    return chunk_result

# 初始化全局变量（多进程需要）
lon = np.arange(-180, 181, 0.1)
lat = np.arange(-90, 91, 0.1) 
land_shp_fname = shpreader.natural_earth(resolution='110m', category='physical', name='land')
land_geom = unary_union(list(shpreader.Reader(land_shp_fname).geometries()))
prepared_land = prep(land_geom)

# 拆分行范围为4个块（根据CPU核心数调整）
num_chunks = 4
j_chunks = np.array_split(range(len(lat)-1), num_chunks)

# 并行处理
with Pool(num_chunks) as pool:
    results = pool.map(process_chunk, j_chunks)

# 合并结果
grid_names = np.vstack(results)

注意事项 #

• 需根据CPU核心数调整num_chunks，避免进程过多导致调度开销
• 多进程的优势在网格越大时越明显，1000×1000的网格可提升2-3倍速度

额外优化建议 #

• 如果对精度要求不高，可以用land_geom.simplify(tolerance=0.1)简化陆地几何，减少几何复杂度，进一步加速判断
• 避免使用unary_union合并所有陆地多边形（如果不需要），直接传入多个多边形给栅格化函数，能减少合并的开销

内容的提问来源于stack exchange，提问作者zxdawn