Search code examples
pythonrasterpython-xarrayrasterio

Update raster values using Python

I'm trying to read in a raster file. It's a 32-bit float raster, with either values 1 or no data. I want to update the values of 1 to 10, and write it out again (probably as a UNIT8 data type?). Here is my attempt:

import rioxarray
import numpy
my_rast = rioxarray.open_rasterio("my_file.tif", masked=True, cache=False)
my_rast[numpy.where(my_rast == 1)] = 10
my_rast.rio.set_nodata(255)
my_rast.rio.to_raster("output.tif", compress='lzw', num_threads='all_cpus', tiled=True, 
dtype='uint8', driver="GTiff", predictor=2, windowed=True)

However the fourth line never seems to finish (maybe as it's a fairly large raster?). I'm not sure what I'm doing wrong.

Here is the result of print(my_rast):

<xarray.DataArray (band: 1, y: 1140, x: 1053)>
[1200420 values with dtype=float64]
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 9.412 9.412 9.412 9.413 ... 9.703 9.704 9.704 9.704
  * y            (y) float64 47.32 47.32 47.32 47.32 ... 47.0 47.0 47.0 47.0
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0
Solution

  • You can’t use a 3D positional index like this:

    my_rast[numpy.where(my_rast == 1)] = 10

    Instead, use xarray.DataArray.where:

    my_rast = my_rast.where(my_rast != 1, 10)

    Indexing in xarray works slightly differently from numpy when providing array-like objects as indices - it’s worth giving the docs on Vectorized Indexing and More Advanced Indexing a read.