# Extracting data within geometry (shape)

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.

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()

• Check this comment on a similar issue posted on xarray GitHub page. Commented Jul 8, 2019 at 15:56

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.

• Thanks, @snowman2. The code worked great, I just have 2 questions/clarifications. (1) The code is good for a tiff but does not work for a netCDF (which has 'lat' and 'lon' in place of x and y). Do you happen to know any way around that (as i was not able to do it)? (2) My data is really high resolution (250 m), that's why I avoid to plot it with lots of data points (near to half a billion). But even after clipping the dataset, I am still not able to plot it, as the plot just shows the colorbar with an empty plot. But the analysis works fine. An example is pasted above (in the question). Commented Jul 11, 2019 at 6:20
• I will create an issue to allow users to define their xy coordinates. Currently you can get around it by renaming your coordinates and dimensions from lat lon to latitude and longitude. Commented Jul 11, 2019 at 11:37
• Side note, I think geoviews should be able to handle plotting larger datasets. I would recommend looking into that. Commented Jul 11, 2019 at 11:38
• Your dimensions are still very large on the clipped dataset - I don't think matplotlib can handle it. I think geoviews (geoviews.org) should be able to handle it. Commented Jul 11, 2019 at 11:45
• >>> xds = xarray.Dataset(coords={"lon": [1, 2], "lat": [3, 4]}) >>> xds <xarray.Dataset> Dimensions: (lat: 2, lon: 2) Coordinates: * lon (lon) int64 1 2 * lat (lat) int64 3 4 Data variables: empty >>> xds.rename({"lon": "longitude", "lat": "latitude"}) <xarray.Dataset> Dimensions: (latitude: 2, longitude: 2) Coordinates: * longitude (longitude) int64 1 2 * latitude (latitude) int64 3 4 Data variables: empty Commented Jul 11, 2019 at 11:48