NumPy interfacing

This example demonstrates NumPy interfacing with rasters on xarray.DataArray objects with the rst accessor. See NumPy functions and reductions for more details.

We open a raster.

import geoutils as gu

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

rast.rst.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.values)
# Estimate the orientation in degrees casting to 2D
aspect_values = np.arctan2(-gradient_x, gradient_y)
aspect = rast.copy(data=(aspect_values * 180 / np.pi) + np.pi)

aspect.rst.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 boolean xarray.DataArray.

mask = np.logical_and(np.logical_and(aspect > -45, aspect < 45), rast > 3000)
mask

<xarray.DataArray (y: 618, x: 539)> Size: 333kB
array([[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]], shape=(618, 539))
Coordinates:
  * y            (y) float64 5kB 4.852e+06 4.852e+06 ... 4.834e+06 4.834e+06
  * x            (x) float64 4kB 6.272e+05 6.272e+05 ... 6.433e+05 6.433e+05
    band         int64 8B 1
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:             Area
    scale_factor:              1.0
    add_offset:                0.0
    overview_factors:          []
    overview_selected_level:   -1
    overview_selected_factor:  1
    block_shape:               (3, 539)

xarray.DataArray

• y: 618
• x: 539

• False False False False False False ... False False False False False

  array([[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]], shape=(618, 539))

• Coordinates: (4)
  • y
    (y)
    float64
    4.852e+06 4.852e+06 ... 4.834e+06

    array([4852070., 4852040., 4852010., ..., 4833620., 4833590., 4833560.],
          shape=(618,))

  • x
    (x)
    float64
    6.272e+05 6.272e+05 ... 6.433e+05

    array([627190., 627220., 627250., ..., 643270., 643300., 643330.], shape=(539,))

  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 18S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32718"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 18S

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -75.0

    false_easting :

    500000.0

    false_northing :

    10000000.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 18S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32718"]]

    GeoTransform :

    627175.0 30.0 0.0 4852085.0 0.0 -30.0

    array(0)

• Attributes: (7)

  AREA_OR_POINT :

  Area

  scale_factor :

  1.0

  add_offset :

  0.0

  overview_factors :

  []

  overview_selected_level :

  -1

  overview_selected_factor :

  1

  block_shape :

  (3, 539)

We plot the mask.

mask.rst.plot()