5.1 Zpracování digitálního modelu terénu pomocí Rasterio a NumPy¶

In [1]:

import numpy as np
import rasterio

import numpy as np import rasterio

In [2]:

dem = rasterio.open('/mnt/repository/155UZPR/12/dmt100.tif')
meta = dem.meta
meta

dem = rasterio.open('/mnt/repository/155UZPR/12/dmt100.tif') meta = dem.meta meta

Out[2]:

{'driver': 'GTiff',
 'dtype': 'uint16',
 'nodata': 65535.0,
 'width': 4729,
 'height': 2921,
 'count': 1,
 'crs': CRS.from_epsg(5514),
 'transform': Affine(100.0, 0.0, -904600.0,
        0.0, -100.0, -935200.0)}

In [3]:

dem.indexes

dem.indexes

Out[3]:

(1,)

In [4]:

data = dem.read(1)
data

data = dem.read(1) data

Out[4]:

array([[65535, 65535, 65535, ..., 65535, 65535, 65535],
       [65535, 65535, 65535, ..., 65535, 65535, 65535],
       [65535, 65535, 65535, ..., 65535, 65535, 65535],
       ...,
       [65535, 65535, 65535, ..., 65535, 65535, 65535],
       [65535, 65535, 65535, ..., 65535, 65535, 65535],
       [65535, 65535, 65535, ..., 65535, 65535, 65535]], dtype=uint16)

Vizualizace¶

In [5]:

import matplotlib.pyplot as plt

# data[data == meta['nodata']] = 0 # replace no-data value with zero
# data[np.where(data == meta['nodata'])] = 0

plt.figure(figsize = (16, 12))
plt.imshow(data, cmap='terrain')
plt.colorbar()

import matplotlib.pyplot as plt # data[data == meta['nodata']] = 0 # replace no-data value with zero # data[np.where(data == meta['nodata'])] = 0 plt.figure(figsize = (16, 12)) plt.imshow(data, cmap='terrain') plt.colorbar()

Out[5]:

<matplotlib.colorbar.Colorbar at 0x74b6c88905c0>

Úkol: Upravte kód, aby buňky s prázdnými hodnotami byly vykresleny bílou barvou.

Reklasifikace¶

In [6]:

dem_copy = dem.read(1).copy()

dem_copy[np.where(dem_copy <= 200)] = 1
dem_copy[np.where((200 < dem_copy) & (dem_copy <= 400))] = 2
dem_copy[np.where((400 < dem_copy) & (dem_copy <= 600))] = 3
dem_copy[np.where((600 < dem_copy) & (dem_copy <= 800))] = 4
dem_copy[np.where(dem_copy > 800)] = 5

dem_copy = dem.read(1).copy() dem_copy[np.where(dem_copy <= 200)] = 1 dem_copy[np.where((200 < dem_copy) & (dem_copy <= 400))] = 2 dem_copy[np.where((400 < dem_copy) & (dem_copy <= 600))] = 3 dem_copy[np.where((600 < dem_copy) & (dem_copy <= 800))] = 4 dem_copy[np.where(dem_copy > 800)] = 5

In [7]:

np.unique(dem_copy)

np.unique(dem_copy)

Out[7]:

array([1, 2, 3, 4, 5], dtype=uint16)

In [8]:

plt.figure(figsize = (16, 12))
plt.imshow(dem_copy, cmap='terrain')
plt.colorbar()

plt.figure(figsize = (16, 12)) plt.imshow(dem_copy, cmap='terrain') plt.colorbar()

Out[8]:

<matplotlib.colorbar.Colorbar at 0x74b6c88eec60>

Topografické analýzy¶

Zdroj pro matematický aparát: https://www.aazuspan.dev/blog/terrain-algorithms-from-scratch/

Sklon svahu¶

In [9]:

def slope(dem):
    dy, dx = np.gradient(dem)
    slope = np.arctan(np.sqrt(dx ** 2 + dy ** 2)) * 180 / np.pi
    
    return slope

def slope(dem): dy, dx = np.gradient(dem) slope = np.arctan(np.sqrt(dx ** 2 + dy ** 2)) * 180 / np.pi return slope

In [10]:

plt.figure(figsize = (16, 12))
plt.imshow(slope(dem.read(1)), cmap='gray')
plt.colorbar()

plt.figure(figsize = (16, 12)) plt.imshow(slope(dem.read(1)), cmap='gray') plt.colorbar()

Out[10]:

<matplotlib.colorbar.Colorbar at 0x74b6c542eae0>

Orientace svahu¶

In [11]:

# counter-clockwise from East
def aspect(dem):
    dy, dx = np.gradient(dem)
    aspect = np.arctan2(-dy, -dx) * 180 / np.pi + 180
    
    return aspect

# counter-clockwise from East def aspect(dem): dy, dx = np.gradient(dem) aspect = np.arctan2(-dy, -dx) * 180 / np.pi + 180 return aspect

In [12]:

plt.figure(figsize = (16, 12))
plt.imshow(aspect(dem.read(1)), cmap='gist_rainbow')
plt.colorbar()

plt.figure(figsize = (16, 12)) plt.imshow(aspect(dem.read(1)), cmap='gist_rainbow') plt.colorbar()

Out[12]:

<matplotlib.colorbar.Colorbar at 0x74b6c54d63c0>

Zkusme si orientaci pro lepší vizualizaci překlasifikovat

In [13]:

a = aspect(dem.read(1))

a[np.where(a <= 90)] = 1 # E-N
a[np.where((90 < a) & (a <= 180))] = 2 # N-W
a[np.where((180 < a) & (a <= 270))] = 3 # W-S
a[np.where((270 < a))] = 4 # S-E

plt.figure(figsize = (16, 12))
plt.imshow(a, cmap='gist_rainbow')
plt.colorbar()

a = aspect(dem.read(1)) a[np.where(a <= 90)] = 1 # E-N a[np.where((90 < a) & (a <= 180))] = 2 # N-W a[np.where((180 < a) & (a <= 270))] = 3 # W-S a[np.where((270 < a))] = 4 # S-E plt.figure(figsize = (16, 12)) plt.imshow(a, cmap='gist_rainbow') plt.colorbar()

Out[13]:

<matplotlib.colorbar.Colorbar at 0x74b6c89225a0>

Převzorkování dat¶

In [14]:

from rasterio.enums import Resampling


upscale_factor = 0.1
resampled = dem.read(
    1,
    out_shape=(
        dem.count,
        int(dem.height * upscale_factor),
        int(dem.width * upscale_factor)
    ),
    resampling=Resampling.bilinear
)
print(dem.height, dem.width, resampled.shape)
a = aspect(resampled)

a[np.where(a <= 90)] = 1
a[np.where((90 < a) & (a <= 180))] = 2
a[np.where((180 < a) & (a <= 270))] = 3
a[np.where((270 < a))] = 4

plt.figure(figsize = (16, 12))
plt.imshow(a, cmap='gist_rainbow')
plt.colorbar()

from rasterio.enums import Resampling upscale_factor = 0.1 resampled = dem.read( 1, out_shape=( dem.count, int(dem.height * upscale_factor), int(dem.width * upscale_factor) ), resampling=Resampling.bilinear ) print(dem.height, dem.width, resampled.shape) a = aspect(resampled) a[np.where(a <= 90)] = 1 a[np.where((90 < a) & (a <= 180))] = 2 a[np.where((180 < a) & (a <= 270))] = 3 a[np.where((270 < a))] = 4 plt.figure(figsize = (16, 12)) plt.imshow(a, cmap='gist_rainbow') plt.colorbar()

2921 4729 (292, 472)

Out[14]:

<matplotlib.colorbar.Colorbar at 0x74b6c52398b0>

Stínovaný reliéf¶

In [15]:

def hillshade(array, azimuth, angle_altitude):
    azimuth = 360.0 - azimuth 
    
    y, x = np.gradient(array)
    slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
    aspect = np.arctan2(-y, -x)
    azm_rad = azimuth * np.pi/180. #azimuth in radians
    alt_rad = angle_altitude * np.pi/180. #altitude in radians
 
    shaded = np.sin(alt_rad) * np.sin(slope) + np.cos(alt_rad) * np.cos(slope) * np.cos((azm_rad - np.pi/2.) - aspect)
    
    return 255 * (shaded + 1) / 2

plt.figure(figsize = (16, 12))
plt.imshow(hillshade(dem.read(1), 0, 0), cmap='gray')

def hillshade(array, azimuth, angle_altitude): azimuth = 360.0 - azimuth y, x = np.gradient(array) slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y)) aspect = np.arctan2(-y, -x) azm_rad = azimuth * np.pi/180. #azimuth in radians alt_rad = angle_altitude * np.pi/180. #altitude in radians shaded = np.sin(alt_rad) * np.sin(slope) + np.cos(alt_rad) * np.cos(slope) * np.cos((azm_rad - np.pi/2.) - aspect) return 255 * (shaded + 1) / 2 plt.figure(figsize = (16, 12)) plt.imshow(hillshade(dem.read(1), 0, 0), cmap='gray')

Out[15]:

<matplotlib.image.AxesImage at 0x74b6c5239910>

In [16]:

fig, ax = plt.subplots(1, 1, figsize=(16, 12))
data = dem.read(1)
data[np.where(data == meta['nodata'])] = 0
im1 = plt.imshow(data, cmap='terrain')
cbar = plt.colorbar()
cbar.set_label('Elevation, m', rotation=270, labelpad=20)

im2 = plt.imshow(hillshade(dem.read(1), 0, 0), cmap='gray', alpha=0.5)
ax = plt.gca()
plt.grid('on')
plt.title('Hillshade + DEM')

fig, ax = plt.subplots(1, 1, figsize=(16, 12)) data = dem.read(1) data[np.where(data == meta['nodata'])] = 0 im1 = plt.imshow(data, cmap='terrain') cbar = plt.colorbar() cbar.set_label('Elevation, m', rotation=270, labelpad=20) im2 = plt.imshow(hillshade(dem.read(1), 0, 0), cmap='gray', alpha=0.5) ax = plt.gca() plt.grid('on') plt.title('Hillshade + DEM')

Out[16]:

Text(0.5, 1.0, 'Hillshade + DEM')