Python GeoDataFrame – Extracting Data Within Geometry (Shape)

Question

I have a NetCDF dataset which I normally call using the xarray function. I have recently been doing some transect analysis on the dataset, which looks something like this (the area between the horizontal blue lines are the transect for analysis):

For now, I am only able to extract the data, which either falls within the vertical or horizontal transect (as they are easy to extract.

But now I want to select a certain shape within the data-set, on which I would like to analyse. And these shapes might not be symmetrical or in a straight line. Some examples are shown below:

Basically, these shapes can be irregular (with known co-ordinates of each point).

Is it possible to extract the dataset/values (with co-ordinates), specifically only for the region of interest?

These steps are easy to perform just using ArcMap's Clip function or Google Earth Engine clip function. But I am not able to use this in Python (as I want solely use Python for all steps). Can someone provide some suggestions on how to this? If anyone knows any package, that has great integration with xarray, that would be great.

Additional Info (after import rioxarray):

Properties of output raster:

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

Treecover_clipped.plot()

Answer 1

I believe what you are looking for is rioxarray. An example of what you want to do is at: https://corteva.github.io/rioxarray/stable/examples/clip_geom.html

import rioxarray

geometries = [
    {
        'type': 'Polygon',
        'coordinates': [[
            [425499.18381405267, 4615331.540546387],
            [425499.18381405267, 4615478.540546387],
            [425526.18381405267, 4615478.540546387],
            [425526.18381405267, 4615331.540546387],
            [425499.18381405267, 4615331.540546387]
        ]]
    }
]
clipped = xds.rio.clip(geometries, xds.rio.crs)

Note that you can have the input geometries in any CRS. The example above has the geometries in the same projection as the dataset. If you have latitude and longitude values, you just modify the second argument to be "epsg:4326".

clipped = xds.rio.clip(geometries, "epsg:4326")

Also, if your CRS is not able to be determined on your xarray dataset, you will need to set it with set_crs:

xds.rio.set_crs("epsg:4326")

You can check if it is able to be determined with:

xds.rio.crs

Hopefully this helps.