Support of pythonic operators

GeoUtils integrates pythonic operators for shorter, more intuitive code, and to perform arithmetic and logical operations consistently.

These operators work on Rasters much as they would on ndarrays, with some more details.

Arithmetic of Raster classes

Arithmetic operators (+, -, /, //, *, **, %) can be used on a Raster in combination with any other Raster, ndarray or number.

For an operation with another Raster, the georeferencing (crs and transform) must match. For another ndarray, the shape must match. The operation always returns a Raster.

import geoutils as gu
import rasterio as rio
import pyproj
import numpy as np

# Create a random 3 x 3 masked array
np.random.seed(42)
arr = np.random.randint(0, 255, size=(3, 3), dtype="uint8")
mask = np.random.randint(0, 2, size=(3, 3), dtype="bool")
ma = np.ma.masked_array(data=arr, mask=mask)

# Create an example raster
rast = gu.Raster.from_array(
       data = ma,
       transform = rio.transform.from_bounds(0, 0, 1, 1, 3, 3),
       crs = pyproj.CRS.from_epsg(4326),
       nodata = 255
    )

rast

Raster(
  data=[[102 -- --]
        [-- 179 61]
        [234 203 --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=255)

# Arithmetic with a number
rast + 1

Raster(
  data=[[103 -- --]
        [-- 180 62]
        [235 204 --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=255)

# Arithmetic with an array
rast / arr

Raster(
  data=[[1.0 -- --]
        [-- 1.0 1.0]
        [1.0 1.0 --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=255)

# Arithmetic with a raster
rast - (rast**0.5)

Raster(
  data=[[91.90049506163793 -- --]
        [-- 165.62091183974036 53.189750324093346]
        [218.70294145922165 188.752193151225 --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=255)

If an unmasked ndarray is passed, it will internally be cast into a MaskedArray to respect the propagation of nodata values. Additionally, the dtype are also reconciled as they would for ndarray, following standard NumPy coercion rules.

Logical comparisons cast to Mask

Logical comparison operators (==, !=, >=, >, <=, <) can be used on a Raster, also in combination with any other Raster, ndarray or number.

Those operation always return a Mask, a subclass of Raster with a boolean MaskedArray as data.

# Logical comparison with a number
mask = rast > 100
mask

Mask(
  data=[[True -- --]
        [-- True False]
        [True True --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326)

Note

A Mask’s data remains a MaskedArray. Therefore, it still maps invalid values through its mask, but has no associated nodata.

Logical bitwise operations on Mask

Logical bitwise operators (~, &, |, ^) can be used to combine a Mask with another Mask, and always output a Mask.

# Logical bitwise operation between masks
mask = (rast > 100) & ((rast % 2) == 0)
mask

Mask(
  data=[[True -- --]
        [-- False False]
        [True False --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326)

Indexing a Raster with a Mask

Finally, indexing and index assignment operations ([], []=) are both supported by Rasters.

For indexing, they can be passed either a Mask with the same georeferencing, or a boolean ndarray of the same shape. For assignment, either a Raster with the same georeferencing, or any ndarray of the same shape is expected.

When indexing, a flattened MaskedArray is returned with the indexed values of the Mask excluding those masked in its data’s MaskedArray (for instance, nodata values present during a previous logical comparison). To bypass this behaviour, simply index without the mask using Mask.data.data.

# Indexing the raster with the previous mask
rast[mask]

masked_array(data=[102, 234],
             mask=[False, False],
       fill_value=255,
            dtype=uint8)