Python: использование полигонов для создания маски на заданной 2d сетке


У меня есть несколько полигонов (канадские провинции), считанных с помощью GeoPandas, и я хочу использовать их для создания маски для применения к сетчатым данным на 2-d широтно-долготной сетке (считывается из файла netcdf с помощью iris). Конечная цель будет состоять в том, чтобы иметь только данные для данной провинции, а остальные данные замаскированы. Таким образом, маска будет 1 для ячеек сетки в пределах провинции и 0 или NaN для ячеек сетки за пределами провинции.


Полигоны могут можно получить из шейп-файла здесь: https://www.dropbox.com/s/o5elu01fetwnobx/CAN_adm1.shp?dl=0

Файл netcdf, который я использую, можно загрузить здесь: https://www.dropbox.com/s/kxb2v2rq17m7lp7/t2m.20090815.nc?dl=0


Я полагаю, что здесь есть два подхода, но я борюсь с обоими:

1) Используйте полигон для создания маски на сетке широта-долгота, чтобы ее можно было применить ко множеству файлов данных за пределами python (предпочтительно)

2) Используйте полигон для маскировки данных, которые были считаны, и извлекайте только данные внутри области интереса, чтобы работать с ними в интерактивном режиме.

Мой код до сих пор:

import iris
import geopandas as gpd

#read the shapefile and extract the polygon for a single province
#(province names stored as variable 'NAME_1')
Canada=gpd.read_file('CAN_adm1.shp')
BritishColumbia=Canada[Canada['NAME_1'] == 'British Columbia']

#get the latitude-longitude grid from netcdf file
cubelist=iris.load('t2m.20090815.nc')
cube=cubelist[0]
lats=cube.coord('latitude').points
lons=cube.coord('longitude').points

#create 2d grid from lats and lons (may not be necessary?)
[lon2d,lat2d]=np.meshgrid(lons,lats)

#HELP!

Большое Спасибо за любую помощь или совет.


Обновление: следуя замечательному решению от @DPeterK ниже, мои исходные данные могут быть замаскированы, давая следующее:

данные о температуре из файла netcdf, замаскированные с помощью шейп-файла канадских провинций

2 3

2 ответа:

Похоже, вы хорошо начали! Геометрия, загруженная из шейп-файлов, предоставляет различные методы геопространственного сравнения, и в этом случае вам нужен метод contains. Вы можете использовать это, чтобы проверить каждую точку в горизонтальной сетке куба на предмет того, что она содержится в геометрии Британской Колумбии. (Обратите внимание, что это Не быстрая операция!) Вы можете использовать это сравнение для построения массива двумерных масок, который может быть применен к данным вашего куба или использован другими способами.

Я написал Функция Python для выполнения вышеописанного-она берет куб и геометрию и создает маску для (заданных) горизонтальных координат Куба, и применяет маску к данным Куба. Функция находится ниже:

def geom_to_masked_cube(cube, geometry, x_coord, y_coord,
                        mask_excludes=False):
    """
    Convert a shapefile geometry into a mask for a cube's data.

    Args:

    * cube:
        The cube to mask.
    * geometry:
        A geometry from a shapefile to define a mask.
    * x_coord: (str or coord)
        A reference to a coord describing the cube's x-axis.
    * y_coord: (str or coord)
        A reference to a coord describing the cube's y-axis.

    Kwargs:

    * mask_excludes: (bool, default False)
        If False, the mask will exclude the area of the geometry from the
        cube's data. If True, the mask will include *only* the area of the
        geometry in the cube's data.

    .. note::
        This function does *not* preserve lazy cube data.

    """
    # Get horizontal coords for masking purposes.
    lats = cube.coord(y_coord).points
    lons = cube.coord(x_coord).points
    lon2d, lat2d = np.meshgrid(lons,lats)

    # Reshape to 1D for easier iteration.
    lon2 = lon2d.reshape(-1)
    lat2 = lat2d.reshape(-1)

    mask = []
    # Iterate through all horizontal points in cube, and
    # check for containment within the specified geometry.
    for lat, lon in zip(lat2, lon2):
        this_point = gpd.geoseries.Point(lon, lat)
        res = geometry.contains(this_point)
        mask.append(res.values[0])

    mask = np.array(mask).reshape(lon2d.shape)
    if mask_excludes:
        # Invert the mask if we want to include the geometry's area.
        mask = ~mask
    # Make sure the mask is the same shape as the cube.
    dim_map = (cube.coord_dims(y_coord)[0],
               cube.coord_dims(x_coord)[0])
    cube_mask = iris.util.broadcast_to_shape(mask, cube.shape, dim_map)

    # Apply the mask to the cube's data.
    data = cube.data
    masked_data = np.ma.masked_array(data, cube_mask)
    cube.data = masked_data
    return cube

Если вам просто нужна 2D-маска, Вы можете вернуть ее до того, как вышеупомянутая функция применит ее к кубу.

Чтобы использовать эту функцию в исходном коде, добавьте в конце кода следующее:

geometry = BritishColumbia.geometry
masked_cube = geom_to_masked_cube(cube, geometry,
                                  'longitude', 'latitude',
                                  mask_excludes=True)

Если это ничего не маскирует, это может означать, что ваш куб и геометрия определяются на разных уровнях. То есть координата долготы вашего Куба находится в диапазоне от 0°-360°, и если значения долготы геометрии находятся в диапазоне от -180°-180°, то тест сдерживания никогда не вернется True. Это можно исправить, изменив экстенты Куба следующим образом:

cube = cube.intersection(longitude=(-180, 180))

Я нашел альтернативное решение отличному, опубликованному @DPeterK выше, которое дает тот же результат. Он использует matplotlib.path для проверки, содержатся ли точки в внешних координатах, описанных геометрией, загруженной из файла формы. я публикую это, потому что этот метод ~в 10 раз быстрее, чем тот, который дает @DPeterK (2:23 минуты против 25:56 минут).Я не знаю, что предпочтительнее: элегантное решение или быстрое решение с применением грубой силы. Возможно один можете иметь оба?!

Одним из осложнений этого метода является то, что некоторые геометрии являются Мультиполигонами , то есть форма состоит из нескольких меньших полигонов (в данном случае провинция Британская Колумбия включает острова у западного побережья, которые не могут быть описаны координатами материковой части Британской Колумбии полигона ). У Мультиполигона нет внешних координат, но есть отдельные полигоны, поэтому каждый из них должен рассматриваться индивидуально. Я обнаружил, что Наиболее точным решением этой проблемы было использование функции, скопированной с GitHub (https://gist.github.com/mhweber/cf36bb4e09df9deee5eb54dc6be74d26 ), который "взрывает" Мультиполигоны в список отдельных полигонов, которые затем могут быть обработаны отдельно.

Рабочий код изложен ниже вместе с моей документацией. Извините, что это не самый элегантный код - я относительно новичок в Python, и я уверен, что есть много ненужных циклов/более аккуратных способов делать вещи!

import numpy as np
import iris
import geopandas as gpd
from shapely.geometry import Point
import matplotlib.path as mpltPath
from shapely.geometry.polygon import Polygon
from shapely.geometry.multipolygon import MultiPolygon


#-----


#FIRST, read in the target data and latitude-longitude grid from netcdf file
cubelist=iris.load('t2m.20090815.minus180_180.nc')
cube=cubelist[0]
lats=cube.coord('latitude').points
lons=cube.coord('longitude').points

#create 2d grid from lats and lons
[lon2d,lat2d]=np.meshgrid(lons,lats)

#create a list of coordinates of all points within grid
points=[]

for latit in range(0,241):
    for lonit in range(0,480):
        point=(lon2d[latit,lonit],lat2d[latit,lonit])
        points.append(point)

#turn into np array for later
points=np.array(points)

#get the cube data - useful for later
fld=np.squeeze(cube.data)

#create a mask array of zeros, same shape as fld, to be modified by
#the code below
mask=np.zeros_like(fld)


#NOW, read the shapefile and extract the polygon for a single province
#(province names stored as variable 'NAME_1')
Canada=gpd.read_file('/Users/ianashpole/Computing/getting_province_outlines/CAN_adm_shp/CAN_adm1.shp')
BritishColumbia=Canada[Canada['NAME_1'] == 'British Columbia']


#BritishColumbia.geometry.type reveals this to be a 'MultiPolygon'
#i.e. several (in this case, thousands...) if individual polygons.
#I ultimately want to get the exterior coordinates of the BritishColumbia
#polygon, but a MultiPolygon is a list of polygons and therefore has no
#exterior coordinates. There are probably many ways to progress from here,
#but the method I have stumbled upon is to 'explode' the multipolygon into
#it's individual polygons and treat each individually. The function below
#to 'explode' the MultiPolygon was found here:
#https://gist.github.com/mhweber/cf36bb4e09df9deee5eb54dc6be74d26


#---define function to explode MultiPolygons

def explode_polygon(indata):
    indf = indata
    outdf = gpd.GeoDataFrame(columns=indf.columns)
    for idx, row in indf.iterrows():
        if type(row.geometry) == Polygon:
            #note: now redundant, but function originally worked on
            #a shapefile which could have combinations of individual polygons
            #and MultiPolygons
            outdf = outdf.append(row,ignore_index=True)
        if type(row.geometry) == MultiPolygon:
            multdf = gpd.GeoDataFrame(columns=indf.columns)
            recs = len(row.geometry)
            multdf = multdf.append([row]*recs,ignore_index=True)
            for geom in range(recs):
                multdf.loc[geom,'geometry'] = row.geometry[geom]
            outdf = outdf.append(multdf,ignore_index=True)
    return outdf

#-------


#Explode the BritishColumbia MultiPolygon into its constituents
EBritishColumbia=explode_polygon(BritishColumbia)


#Loop over each individual polygon and get external coordinates
for index,row in EBritishColumbia.iterrows():

    print 'working on polygon', index
    mypolygon=[]
    for pt in list(row['geometry'].exterior.coords):
        print index,', ',pt
        mypolygon.append(pt)


    #See if any of the original grid points read from the netcdf file earlier
    #lie within the exterior coordinates of this polygon
    #pth.contains_points returns a boolean array (true/false), in the
    #shape of 'points'
    path=mpltPath.Path(mypolygon)
    inside=path.contains_points(points)


    #find the results in the array that were inside the polygon ('True')
    #and set them to missing. First, must reshape the result of the search
    #('points') so that it matches the mask & original data
    #reshape the result to the main grid array
    inside=np.array(inside).reshape(lon2d.shape)
    i=np.where(inside == True)
    mask[i]=1


print 'fininshed checking for points inside all polygons'


#mask now contains 0's for points that are not within British Columbia, and
#1's for points that are. FINALLY, use this to mask the original data
#(stored as 'fld')
i=np.where(mask == 0)
fld[i]=np.nan

#Done.