Python arithmetic

This example demonstrates arithmetic operations using raster arithmetic on Rasters. See Support of pythonic operators for more details.

We open a raster

import geoutils as gu

filename_rast = gu.examples.get_path("everest_landsat_b4")
rast = gu.Raster(filename_rast)
rast

Raster(
  data=not_loaded; shape on disk (1, 655, 800); will load (655, 800)
  transform=| 30.00, 0.00, 478000.00|
            | 0.00,-30.00, 3108140.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:32645
  nodata=None)

rast.plot(cmap="Greys_r")

Performing arithmetic operations implicitly loads the data.

rast = (rast + 1.0) ** 0.5 / 5
rast.plot(cmap="Greys_r")

Important

Arithmetic operations cast to new dtypes automatically following NumPy coercion rules. If we had written (rast + 1), this calculation would have conserved the original numpy.uint8 dtype of the raster.

Logical comparison operations will naturally cast to a Mask.

mask = rast == 200
mask

Mask(
  data=[[False False False ... False False False]
        [False False False ... False False False]
        [False False False ... False False False]
        ...
        [False False False ... False False False]
        [False False False ... False False False]
        [False False False ... False False False]]
  transform=| 30.00, 0.00, 478000.00|
            | 0.00,-30.00, 3108140.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:32645)

Masks support python logical operators to be combined together

mask = (rast >= 3) | (rast % 2 == 0) & (rast != 80)
mask.plot()

Finally, Masks can be used for indexing and assigning to a Rasters

values = rast[mask]