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

In [1]:

import numpy as np
import rasterio
import matplotlib.pyplot as plt

import numpy as np import rasterio import matplotlib.pyplot as plt

In [2]:

data = 'data/dem.tif'

data = 'data/dem.tif'

In [3]:

dem = rasterio.open(data)

dem = rasterio.open(data)

In [4]:

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

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

Out[4]:

<matplotlib.colorbar.Colorbar at 0x7fde2227fc40>

Reklasifikace¶

In [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

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 [6]:

np.unique(dem_copy)

np.unique(dem_copy)

Out[6]:

array([1., 2., 3., 4.], dtype=float32)

In [7]:

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[7]:

<matplotlib.colorbar.Colorbar at 0x7fde3a7ac190>

Sklon svahu¶

In [8]:

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 [9]:

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[9]:

<matplotlib.colorbar.Colorbar at 0x7fde3ac81120>

Orientace svahu¶

In [10]:

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

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

In [11]:

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[11]:

<matplotlib.colorbar.Colorbar at 0x7fde3ab4a890>

Zkusme si orientaci pro lepší vizualizaci překlasifikovat

In [12]:

a = aspect(dem.read(1))

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()

a = aspect(dem.read(1)) 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()

Out[12]:

<matplotlib.colorbar.Colorbar at 0x7fde3aa34430>

In [13]:

from rasterio.enums import Resampling


upscale_factor = 0.02
resampled = dem.read(
    1,
    out_shape=(
        dem.count,
        int(dem.height * upscale_factor),
        int(dem.width * upscale_factor)
    ),
    resampling=Resampling.bilinear
)

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.02 resampled = dem.read( 1, out_shape=( dem.count, int(dem.height * upscale_factor), int(dem.width * upscale_factor) ), resampling=Resampling.bilinear ) 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()

Out[13]:

<matplotlib.colorbar.Colorbar at 0x7fde3aaf13c0>

Stínovaný reliéf¶

In [14]:

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[14]:

<matplotlib.image.AxesImage at 0x7fde3a9b5180>

In [15]:

fig, ax = plt.subplots(1, 1, figsize=(16, 12))
im1 = plt.imshow(dem.read(1), 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)) im1 = plt.imshow(dem.read(1), 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[15]:

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