4

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

example1

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:

example2

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

enter image description here

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/html/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). – Ep1c1aN Jul 11 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. – snowman2 Jul 11 at 11:37
  • Side note, I think geoviews should be able to handle plotting larger datasets. I would recommend looking into that. – snowman2 Jul 11 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. – snowman2 Jul 11 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 – snowman2 Jul 11 at 11:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.