Python – How to Clip an Xarray to a Smaller Extent Given Lat/Lon Coordinates

clipnetcdfpythonrasterioxarray

I have a netcdf file showing elevation over Europe derived from GTopo30. Turns out that the extent of this raster file is too big to combine it with other layers, so I would like to crop it to a given extent specified in lat/lon coordinates. I have seen some examples here and here trying to make this operation I am describing, but it seems they are not working in this case.

The raster file description looks like this:

<xarray.DataArray (band: 14, y: 5520, x: 8400)>
[649152000 values with dtype=float64]
Coordinates:
  * band     (band) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14
  * y        (y) float64 71.0 70.99 70.98 70.97 70.96 ... 25.03 25.02 25.01 25.0
  * x        (x) float64 -25.0 -24.99 -24.98 -24.97 ... 44.97 44.98 44.99 45.0
Attributes:
    transform:     (0.0083333333, 0.0, -25.000139509, 0.0, -0.0083333333, 70....
    crs:           +init=epsg:4326
    res:           (0.0083333333, 0.0083333333)
    is_tiled:      0
    nodatavals:    (-1.7e+308, -1.7e+308, -1.7e+308, -1.7e+308, -1.7e+308, -1...
    scales:        (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1....
    offsets:       (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0....
    descriptions:  ('tpi', 'slope', 'aspect', 'Pvec', 'Qvec', 'alt', 'dist2co...

So, as you can see, the longitude and longitudes are stored under the "x" and "y" coordinates and they come in WGS84 (EPSG:4326) format.

The extent is this:

min_lon = -24.995
min_lat = 25.05
max_lon = 45.50
max_lat = 71.55

I tried to use xarray's filters and subsetting, but it seems I am doing something wrong, since I can't get the required slice:

# Example 1
sel1 = band.sel(x=(band.x < max_lon) | (band.x > min_lon))
sel2 = sel1.sel(y=(sel1.y < max_lat) | (sel1.y > min_lat))

# Example 2
sel1 = band.where((band.x < max_lon) | (band.x > min_lon))
sel2 = sel1.where((sel1.x < max_lat) | (sel1.x > min_lat))

Both of them produce: KeyError: "not all values found in index 'x'"

What is the correct xarray way of cropping a large raster to a given lat/lon window?
Thanks!

Best Answer

You could use rioxarray: https://corteva.github.io/rioxarray/stable/examples/clip_box.html

import rioxarray

min_lon = -24.995
min_lat = 25.05
max_lon = 45.50
max_lat = 71.55

subset = band.rio.clip_box(minx=min_lon, miny=min_lat, maxx=max_lon, maxy=max_lat)

Related Solutions

[GIS] Geotransformation for polar stereographic

Too long for a comment, this is to accompany @Matej's answer.

Add the “.grd” data into ArcGIS.
Use the “Raster to Other Format” function and convert your .grd file into an ESRII GRID format. This is important because most of the raster functions in ArcGIS are accessible for this format only, either that or it usually gets too slow when you use it on other formats.
As it already has the projection file associated to the file. Rather than projecting the new converted data, define its projection. ArcToolbox > Data Management Tools > Projections and Transformations > Define Projection. You can navigate to the pre-defined polar stereo graphic ESRII projection and see if its parameters match the one given in the metadata (it doesn’t), so you can modify it as per @Matej. Only here - rather than modifying, create a new one based on the NPS projection with the central meridian and Latitude of origin changed and save it on disk, then navigate to the new projection and use that when defining your projection. This is because your on the fly modification won't be available later on when you want to use it to set the coordinate system for your data frame, which is necessary in order to align your other boundary files to the same modified NPS system.

Open Data Cube – How to Change Specific Values in an Xarray Depending on Coordinates

If I understand correctly, you want to replace the values for bare soil, vegetation and water everywhere where yyis the value of y.

If so, your mock_array might have the shape (5,3), for each of the 5 timesteps and 3 bare soil/vegetation/water values.

mock_array = np.arange(15).reshape((5,3))

mock_array

array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11],
       [12, 13, 14]])

Values in cube1 can be updated by indexing by the yy position. If the index of yy were 0, then

cube1[:,0,:] = mock_array

cube1

<xarray.DataArray (time: 5, y: 3, x: 3)>
array([[[ 0.,  1.,  2.],
        [ 4.,  4.,  6.],
        [ 4.,  4.,  6.]],

       [[ 3.,  4.,  5.],
        [ 4.,  4.,  6.],
        [ 4.,  4.,  6.]],

       [[ 6.,  7.,  8.],
        [ 6.,  6.,  6.],
        [ 6.,  6.,  6.]],

       [[ 9., 10., 11.],
        [ 4.,  4.,  6.],
        [ 4.,  4.,  6.]],

       [[12., 13., 14.],
        [ 4.,  4.,  6.],
        [ 4.,  4.,  6.]]])
Coordinates:
  * y            (y) float64 6.607e+06 6.608e+06 6.609e+06
  * x            (x) float64 7.091e+05 7.092e+05 7.093e+05
  * time         (time) datetime64[ns] 2014-09-06 2014-09-07 ... 2014-09-10
    spatial_ref  int64 -2147483647

I don't have a more elegant suggestion to finding the positional index than

cube1.y.values.tolist().index(yy)

This will give you the index of yy. You can also put the whole statement straight into the previous one:

cube1[:,cube1.y.values.tolist().index(6607000),:] = mock_array

Best Answer

Related Solutions

[GIS] Geotransformation for polar stereographic

Open Data Cube – How to Change Specific Values in an Xarray Depending on Coordinates

Related Question