Извлечение данных внутри геометрии (фигуры) - PullRequest
Новая
       49

Извлечение данных внутри геометрии (фигуры)

1 голос
/ 05 июля 2019

У меня есть набор данных NetCDF, который я обычно вызываю с помощью функции xarray. Недавно я проводил некоторый анализ разреза для набора данных, который выглядит примерно так (область между горизонтальными синими линиями - это трансект для анализа):

enter image description here

На данный момент я могу извлечь только те данные, которые попадают в вертикальный или горизонтальный разрез (поскольку их легко извлечь.

Но теперь я хочу выбрать определенную фигуру в наборе данных, которую я хотел бы проанализировать. И эти формы могут быть не симметричными или прямолинейными. Некоторые примеры приведены ниже: enter image description here

В принципе, эти формы могут быть неправильными (с известными координатами каждой точки).

Можно ли извлечь набор / значения данных (с координатами), конкретно только для интересующей области?

Эти шаги легко выполнить, используя функцию клипа ARCMaps или функцию клипа Google Earth Engine. Но я не могу использовать это в Python (так как я хочу использовать Python исключительно для всех шагов). Может кто-нибудь дать несколько советов, как это сделать? Если кто-нибудь знает какой-либо пакет, который имеет отличную интеграцию с xarray, это было бы здорово.

1 Ответ

1 голос
/ 11 июля 2019

Я использовал пакет rioxarray, который, кажется, отлично интегрируется с xarray и действительно прост в реализации. 

import xarray as xr
import rioxarray as rx

Treecover  = xr.open_rasterio('/home/chandra/data/Treecover_MOD44B_2000_250m_AMAZON.tif')
[Output]: 
<xarray.DataArray (band: 1, y: 32093, x: 20818)>
[668112074 values with dtype=float64]
Coordinates:
  * band     (band) int64 1
  * y        (y) float64 13.71 13.71 13.71 13.71 ... -58.35 -58.35 -58.36 -58.36
  * x        (x) float64 -81.38 -81.37 -81.37 -81.37 ... -34.63 -34.63 -34.62
Attributes:
    transform:   (0.002245788210298804, 0.0, -81.37613580017715, 0.0, -0.0022...
    crs:         +init=epsg:4326
    res:         (0.002245788210298804, 0.002245788210298804)
    is_tiled:    0
    nodatavals:  (nan,)

geometries = [
    {
        'type': 'Polygon',
        'coordinates': [[
            [-46.23140155225633, -21.53505449239459],
          [-44.91304217725633, -20.221175092759253],
          [-70.22554217725633, 1.5816072875439455],
          [-71.36812030225633, 0.5271132528460204]
        ]]
    }
]

Treecover_clipped = Treecover.rio.clip(geometries, Treecover.rio.crs)
[Output]: 
<xarray.DataArray (band: 1, y: 10293, x: 11779)>
array([[[nan, nan, ..., nan, nan],
        [nan, nan, ..., nan, nan],
        ...,
        [nan, nan, ..., nan, nan],
        [nan, nan, ..., nan, nan]]])
Coordinates:
  * band         (band) int64 1
  * y            (y) float64 1.58 1.578 1.575 1.573 ... -21.53 -21.53 -21.53
  * x            (x) float64 -71.37 -71.36 -71.36 ... -44.92 -44.92 -44.91
    spatial_ref  int64 0
Attributes:
    transform:     (0.0022457882102988043, 0.0, -71.36665774687539, 0.0, -0.0...
    _FillValue:    nan
    grid_mapping:  spatial_ref
...