In [1]:

#mamba install -c conda-forge hvplot datashader geoviews xarray-spatial

#mamba install -c conda-forge hvplot datashader geoviews xarray-spatial

In [24]:

import time

import datashader.colors as dscolors
import rioxarray as rxr
import hvplot.xarray
import xarray as xr
import xrspatial

import time import datashader.colors as dscolors import rioxarray as rxr import hvplot.xarray import xarray as xr import xrspatial

Mounting data¶

In [28]:

%%time

# Define the URL of the DEM file
dem_url = (
    "/vsizip/vsicurl/"
    "https://data.hydrosheds.org/file/hydrosheds-v1-dem/hyd_na_dem_15s.zip/"
    "hyd_na_dem_15s.tif")

# Open the DEM file with rioxarray
dem_da = rxr.open_rasterio(dem_url, masked=True).squeeze()

# Print out details of the raster object
print("Raster details:")
print("Shape:", dem_da.shape)
print("CRS:", dem_da.rio.crs)
dem_da

%%time # Define the URL of the DEM file dem_url = ( "/vsizip/vsicurl/" "https://data.hydrosheds.org/file/hydrosheds-v1-dem/hyd_na_dem_15s.zip/" "hyd_na_dem_15s.tif") # Open the DEM file with rioxarray dem_da = rxr.open_rasterio(dem_url, masked=True).squeeze() # Print out details of the raster object print("Raster details:") print("Shape:", dem_da.shape) print("CRS:", dem_da.rio.crs) dem_da

Raster details:
Shape: (13920, 20640)
CRS: EPSG:4326
CPU times: user 22.7 ms, sys: 146 µs, total: 22.8 ms
Wall time: 19.4 ms

Out[28]:

<xarray.DataArray (y: 13920, x: 20640)>
[287308800 values with dtype=float32]
Coordinates:
    band         int64 1
  * x            (x) float64 -138.0 -138.0 -138.0 -138.0 ... -52.01 -52.01 -52.0
  * y            (y) float64 63.0 62.99 62.99 62.99 ... 5.015 5.01 5.006 5.002
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:       Area
    RepresentationType:  THEMATIC
    scale_factor:        1.0
    add_offset:          0.0
    long_name:           Band_1

xarray.DataArray

• y: 13920
• x: 20640

• ...

  [287308800 values with dtype=float32]

• Coordinates: (4)
  • band
    ()
    int64
    1

    array(1)

  • x
    (x)
    float64
    -138.0 -138.0 ... -52.01 -52.0

    array([-137.997917, -137.99375 , -137.989583, ...,  -52.010417,  -52.00625 ,
            -52.002083])

  • y
    (y)
    float64
    63.0 62.99 62.99 ... 5.006 5.002

    array([62.997917, 62.99375 , 62.989583, ...,  5.010417,  5.00625 ,  5.002083])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    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

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -138.0 0.00416666666666667 0.0 63.00000000000005 0.0 -0.00416666666666667

    array(0)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([-137.99791666666667,          -137.99375, -137.98958333333334,
           -137.98541666666668,          -137.98125, -137.97708333333333,
           -137.97291666666666,          -137.96875, -137.96458333333334,
           -137.96041666666667,
           ...
            -52.03958333333327,  -52.03541666666659,  -52.03124999999993,
           -52.027083333333266,   -52.0229166666666, -52.018749999999926,
            -52.01458333333326,   -52.0104166666666,  -52.00624999999994,
            -52.00208333333326],
          dtype='float64', name='x', length=20640))

  • y
    PandasIndex

    PandasIndex(Index([ 62.99791666666672,  62.99375000000005, 62.989583333333385,
           62.985416666666715,  62.98125000000005,  62.97708333333338,
            62.97291666666672,  62.96875000000005,  62.96458333333339,
            62.96041666666672,
           ...
             5.03958333333334,   5.03541666666667,  5.031250000000007,
            5.027083333333337,  5.022916666666674,  5.018750000000004,
            5.014583333333334,  5.010416666666671,  5.006250000000001,
           5.0020833333333385],
          dtype='float64', name='y', length=13920))

• Attributes: (5)

  AREA_OR_POINT :

  Area

  RepresentationType :

  THEMATIC

  scale_factor :

  1.0

  add_offset :

  0.0

  long_name :

  Band_1

In [26]:

start_time = time.time()
dem_da.hvplot(x='x', y='y', geo=True, rasterize=True, cmap=dscolors.Elevation)

start_time = time.time() dem_da.hvplot(x='x', y='y', geo=True, rasterize=True, cmap=dscolors.Elevation)

Out[26]:

In [17]:

f'Elapsed time: {time.time() - start_time} s'

f'Elapsed time: {time.time() - start_time} s'

Out[17]:

'Elapsed time: 55.768351793289185 s'

Terrain calculations¶

Calculate slope from data¶

In [25]:

%%time
# THIS SLOPE IS INCORRECT BECAUSE THE DATA ARE NOT PROJECTED (units issue with xrspatial)
slope_da = xrspatial.slope(dem_da)
slope_da.hvplot(x='x', y='y', geo=True, rasterize=True, )

%%time # THIS SLOPE IS INCORRECT BECAUSE THE DATA ARE NOT PROJECTED (units issue with xrspatial) slope_da = xrspatial.slope(dem_da) slope_da.hvplot(x='x', y='y', geo=True, rasterize=True, )

CPU times: user 11.3 s, sys: 564 ms, total: 11.8 s
Wall time: 11.8 s

Out[25]:

Calculate aspect from data¶

In [19]:

%%time
aspect_da = xrspatial.aspect(dem_da)
aspect_da.hvplot(x='x', y='y', geo=True, rasterize=True, cmap='colorwheel')

%%time aspect_da = xrspatial.aspect(dem_da) aspect_da.hvplot(x='x', y='y', geo=True, rasterize=True, cmap='colorwheel')

CPU times: user 9.87 s, sys: 421 ms, total: 10.3 s
Wall time: 10.3 s

Out[19]:

Create a cube from those layers!¶

In [22]:

data_cube = xr.Dataset(dict(
    elevation=dem_da,
    slope=slope_da,
    aspect=aspect_da
))
data_cube

data_cube = xr.Dataset(dict( elevation=dem_da, slope=slope_da, aspect=aspect_da )) data_cube

Out[22]:

<xarray.Dataset>
Dimensions:      (x: 20640, y: 13920)
Coordinates:
    band         int64 1
  * x            (x) float64 -138.0 -138.0 -138.0 -138.0 ... -52.01 -52.01 -52.0
  * y            (y) float64 63.0 62.99 62.99 62.99 ... 5.015 5.01 5.006 5.002
    spatial_ref  int64 0
Data variables:
    elevation    (y, x) float32 nan nan nan nan nan nan ... nan nan nan nan nan
    slope        (y, x) float32 nan nan nan nan nan nan ... nan nan nan nan nan
    aspect       (y, x) float32 nan nan nan nan nan nan ... nan nan nan nan nan

xarray.Dataset

• Dimensions:
  • x: 20640
  • y: 13920
• Coordinates: (4)
  • band
    ()
    int64
    1

    array(1)

  • x
    (x)
    float64
    -138.0 -138.0 ... -52.01 -52.0

    array([-137.997917, -137.99375 , -137.989583, ...,  -52.010417,  -52.00625 ,
            -52.002083])

  • y
    (y)
    float64
    63.0 62.99 62.99 ... 5.006 5.002

    array([62.997917, 62.99375 , 62.989583, ...,  5.010417,  5.00625 ,  5.002083])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    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

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -138.0 0.00416666666666667 0.0 63.00000000000005 0.0 -0.00416666666666667

    array(0)

• Data variables: (3)
  • elevation
    (y, x)
    float32
    nan nan nan nan ... nan nan nan nan

    AREA_OR_POINT :

    Area

    RepresentationType :

    THEMATIC

    scale_factor :

    1.0

    add_offset :

    0.0

    long_name :

    Band_1

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

  • slope
    (y, x)
    float32
    nan nan nan nan ... nan nan nan nan

    AREA_OR_POINT :

    Area

    RepresentationType :

    THEMATIC

    scale_factor :

    1.0

    add_offset :

    0.0

    long_name :

    Band_1

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

  • aspect
    (y, x)
    float32
    nan nan nan nan ... nan nan nan nan

    AREA_OR_POINT :

    Area

    RepresentationType :

    THEMATIC

    scale_factor :

    1.0

    add_offset :

    0.0

    long_name :

    Band_1

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([-137.99791666666667,          -137.99375, -137.98958333333334,
           -137.98541666666668,          -137.98125, -137.97708333333333,
           -137.97291666666666,          -137.96875, -137.96458333333334,
           -137.96041666666667,
           ...
            -52.03958333333327,  -52.03541666666659,  -52.03124999999993,
           -52.027083333333266,   -52.0229166666666, -52.018749999999926,
            -52.01458333333326,   -52.0104166666666,  -52.00624999999994,
            -52.00208333333326],
          dtype='float64', name='x', length=20640))

  • y
    PandasIndex

    PandasIndex(Index([ 62.99791666666672,  62.99375000000005, 62.989583333333385,
           62.985416666666715,  62.98125000000005,  62.97708333333338,
            62.97291666666672,  62.96875000000005,  62.96458333333339,
            62.96041666666672,
           ...
             5.03958333333334,   5.03541666666667,  5.031250000000007,
            5.027083333333337,  5.022916666666674,  5.018750000000004,
            5.014583333333334,  5.010416666666671,  5.006250000000001,
           5.0020833333333385],
          dtype='float64', name='y', length=13920))

• Attributes: (0)

In [23]:

data_cube.mean()

data_cube.mean()

Out[23]:

<xarray.Dataset>
Dimensions:      ()
Coordinates:
    band         int64 1
    spatial_ref  int64 0
Data variables:
    elevation    float32 664.2
    slope        float32 88.27
    aspect       float32 171.0

xarray.Dataset

• Dimensions:
• Coordinates: (2)
  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    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

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -138.0 0.00416666666666667 0.0 63.00000000000005 0.0 -0.00416666666666667

    array(0)

• Data variables: (3)
  • elevation
    ()
    float32
    664.2

    array(664.2231, dtype=float32)

  • slope
    ()
    float32
    88.27

    array(88.27206, dtype=float32)

  • aspect
    ()
    float32
    171.0

    array(170.96733, dtype=float32)

• Indexes: (0)
• Attributes: (0)

---

Last update: 2023-11-16