⚠️ Our 0.1 release refactored several early-development functions for long-term stability, to update your code see here. ⚠️
Future changes will come with deprecation warnings! 🙂

NumPy interfacing

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")
numpy interfacing

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)")
numpy interfacing

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.

numpy interfacing

Total running time of the script: (0 minutes 0.400 seconds)

Gallery generated by Sphinx-Gallery