Как найти все включенные точки (лат / лоны) заданных вершин разных областей?

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

Один только этот код не работает, поэтому я попытался заполнить пробелы, создав входные данные.

См. , как создать пример Minimal, Complete и Verifiable

Определяя регион-нарушитель, как показано ниже, работает нижняя часть вашего кода. Это говорит о том, что проблема зависит от координат для области, в том числе определены полюса. Если вы измените их, вы сможете сделать его совместимым с .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.]