NumPy interfacing

This example demonstrates NumPy interfacing with rasters on Rasters. See Masked-array NumPy interface for more details.

We open a raster.

import geoutils as gu

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

rast.plot(cmap="terrain")

The NumPy interface allows to use almost any NumPy function directly on the raster.

import numpy as np

# Get the x and y gradient as 1D arrays
gradient_y, gradient_x = np.gradient(rast)
# Estimate the orientation in degrees casting to 2D
aspect = np.arctan2(-gradient_x, gradient_y)
aspect = (aspect * 180 / np.pi) + np.pi
# We copy the new array into a raster
asp = rast.copy(new_array=aspect)

asp.plot(cmap="twilight", cbar_title="Aspect (degrees)")

Important

For rigorous slope and aspect calculation (matching that of GDAL), check-out our sister package xDEM.

We use NumPy logical operations to isolate the terrain oriented South and above three thousand meters. The rasters will be logically cast to a Mask.

mask = np.logical_and.reduce((asp > -45, asp < 45, rast > 3000))
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, 627175.00|
            | 0.00,-30.00, 4852085.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:32718)

We plot the mask.

mask.plot()