如何在给定不同区域的顶点的情况下找到所有包含的点(纬度/纬度)?

Question

北极有16个区域，我得到了每个区域的纬度/纬度顶点。我需要找到每个区域内存在的所有点。数据分辨率为 1.5x1.5 度，0 至 358.5 经度，90 至 -90 纬度。

我尝试使用 matplotlib.path，我已经能够处理经度从 -180 到 +180 的交叉点，但是北极中部(即北极)似乎有问题。沿着~80度向北绘制似乎无法识别出有点一直通向90。

首先列出"问题区域"(北极中部)的代码，然后列出未越过-180/180线的其他区域的代码。

# regions[i] represents a class of 16 regions
# Problem region!
x, y = np.meshgrid(lons, lats) 
x, y = x.flatten(), y.flatten()
points = np.vstack((x,y)).T
lons_temp = regions[i].lons - 180
lons_temp = np.append(lons_temp, lons_temp[0])
lats_temp = np.append(regions[i].lats, regions[i].lats[0])
cross = np.where(np.diff(np.signbit(lons_temp)))[0]
lons_temp = lons_temp + 180
for c in cross:
cross_point = lons_temp[c]
if cross_point<90:
mid_val = 0
elif cross_point >= 90 and cross_point < 270: 
mid_val = 180
elif cross_point >=270:
mid_val = 360
interp = np.interp(mid_val, lons_temp[c:c+2], lats_temp[c:c+2])
lons_temp = np.insert(lons_temp, c+1, mid_val)
lats_temp = np.insert(lats_temp, c+1, interp)
lons_neg = lons_temp[np.where(lons_temp <= 180)]
lats_neg = lats_temp[np.where(lons_temp <= 180)]
lons_pos = lons_temp[np.where(lons_temp >= 180)]
lats_pos = lats_temp[np.where(lons_temp >= 180)]
gg_neg = np.array([lons_neg, lats_neg])
gg_pos = np.array([lons_pos, lats_pos])
pp_neg = gg_neg.T # lon lat pair
pp_pos = gg_pos.T
p_neg=Path(pp_neg, closed=False)
p_pos=Path(pp_pos, closed=False)
grid_neg = p_neg.contains_points(points)
grid_pos = p_pos.contains_points(points)
grid = np.logical_or(grid_neg, grid_pos)
grid = np.where(points.T[1]>=82.5, True, grid)
latslons1 = np.where(grid==True)[0]
regions[i].included_points = latslons1

# Other regions that do not cross over -180/180
x, y = np.meshgrid(lons, lats) 
x, y = x.flatten(), y.flatten()
points = np.vstack((x,y)).T
gg = np.array([regions[i].lons, regions[i].lats])
pp = gg.T
p=Path(pp)
grid = p.contains_points(points)
latslons1 = np.where(grid==True)[0]
regions[i].included_points = latslons1

我希望有一种方法可以在一个代码块中包含所有区域，它可以处理围绕 -180/180 线的包装以及解决北极中部的问题。下面的链接显示了我试图重现的图像(但是，我确实需要数据点正确进行分析)。https://nsidc.org/data/masie/browse_regions

Answer 1

仅此代码不起作用，因此我尝试通过创建输入数据来填补空白。

请参阅如何创建最小、完整和可验证的示例

通过定义如下违规区域，代码的下半部分就可以工作了。这表明问题取决于定义包括极点在内的区域的坐标。如果您修改这些内容，您应该能够使其符合.contains_points

import numpy as np
from matplotlib import path
# create missing data
lons = np.arange(0,360,1.5)
lats = np.arange(-90,90.1,1.5,)
class polyg:
lons = []
lats=[]
# unproblematic region
region = polyg()
region.lats = np.array([80.,85.,80.])
region.lons = np.array([100.,100.,110.])
regions=[region]

# Including North pole
region = polyg()
region.lats = np.array([80.,80.,90.,90.])
region.lons = np.array([0.,360.,360.,0.])
regions.append(region)
i=1
x, y = np.meshgrid(lons, lats) 
x, y = x.flatten(), y.flatten()
points = np.vstack((x,y)).T
gg = np.array([regions[i].lons, regions[i].lats])
pp = gg.T
p=path.Path(pp)
grid = p.contains_points(points)
latslons1 = np.where(grid==True)[0]
regions[i].included_points = latslons1

print(x[latslons1])
print(y[latslons1])
[  0.    1.5   3.  ... 355.5 357.  358.5]
[81. 81. 81. ... 90. 90. 90.]