<jupyter_start><jupyter_text>MIT LicenseCopyright (c) 2019 <NAME>, https://orcid.org/0000-0001-9626-8615 (ORCID)Build Super-resolution DEM 0.5m from DEM 1m enhanced by one orthophoto image 0.5mSource dataset: Elevation1 - Mashhad, IranElevation1 DSM + Pléiades Ortho 0.5m pan-sharpened (Orthoimage included)https://www.intelligence-airbusds.com/en/9317-sample-imagery-detail?product=18896&keyword=&type=366<jupyter_code>from osgeo import gdal
import os
import numpy as np
import xarray as xr
import pandas as pd
from scipy.ndimage.filters import gaussian_filter
from scipy.spatial import cKDTree

import matplotlib.pyplot as plt
%matplotlib inline
# select work area
def crop_area(raster):
    return raster.sel(x=slice(730000,730500),y=slice(4012000, 4011500))
def crop_sample(raster):
    return raster.sel(x=slice(730000,730100),y=slice(4011700,4011600))<jupyter_output><empty_output><jupyter_text>Load DEM<jupyter_code># DEM image
dem = crop_area(xr.open_rasterio("data/Mashhad-DEM.tif")[0])
del dem.attrs['units']
dem<jupyter_output><empty_output><jupyter_text>Load orthophoto image<jupyter_code># orthophoto image 0.5m for the same area as DEM data above
img = crop_area(xr.open_rasterio("data/7289-40126_Mashhad.tif")[0])<jupyter_output><empty_output><jupyter_text>Regrid DEM 1m on the same grid as orthophoto image 0.5mThis is nearest neighbor interpolation without any data quality enhancement<jupyter_code># define source values
df_dem = dem.to_dataframe(name='dem').dropna().reset_index()
# target grid 
df_grid = img.to_dataframe(name='_').reset_index()[['y','x']]
# nearest neighbor interpolation
tree = cKDTree(list(zip(df_dem.x, df_dem.y)))
distance, indices = tree.query(list(zip(df_grid.x, df_grid.y)), k = 1)
values = df_dem.dem.values[indices]
dem2x = xr.DataArray(values.reshape(img.shape),
                          coords=[img.y,img.x])<jupyter_output><empty_output><jupyter_text>Enhance DEM by orthophoto imageAs explained in the article we need to transfer spatial components 0-5mWith pixel size 0.5m the required filter radius is equal to 10 pixels because 10*0.5m = 5m<jupyter_code># low-pass filter
def raster_filter(src, gamma):
    dst = src.copy()
    dst.values = gaussian_filter(dst.values.astype(np.float32),gamma,mode='nearest')
    return dst

# define spectrum components to transfer
radius = 5/img.res[0]

img_lowpass = raster_filter(img, radius)
img_hipass  = img - img_lowpass

dem2x_lowpass = raster_filter(dem2x, radius)
dem2x_hipass  = dem2x - dem2x_lowpass

# caclulate approximate scale factor for short wavelenghs
scale = (img_hipass.max()-img_hipass.min())/(dem2x_hipass.max()-dem2x_hipass.min())
scale

# super-resolution dataset
dem2x_hires = dem2x_lowpass + img_hipass/scale
fig, ax = plt.subplots(2,3,figsize=(16,9))
((ax1,ax2,ax3),(ax4,ax5,ax6)) = ax

dem.plot(vmin=1270,vmax=1400,cmap='RdBu_r',ax=ax1)
ax1.set_title('DEM 1m', fontsize=18)
img.plot(vmin=0,vmax=160,cmap='RdBu_r',ax=ax2)
ax2.set_title('Ortho 0.5m [Red band]', fontsize=18)

crop_sample(dem).plot(vmin=1270,vmax=1400,cmap='RdBu_r',ax=ax4)
ax4.set_title('')
crop_sample(img).plot(vmin=0,vmax=160,cmap='RdBu_r',ax=ax5)
ax5.set_title('')

title = """Build Super-resolution DEM 0.5m from
Elevation 1m DSM (actually, ~10m) + Pléiades Ortho 0.5m for Mashhad, Iran:
"""

dem2x_hires.plot(vmin=1270,vmax=1400,cmap='RdBu_r',ax=ax3)
ax3.set_title('Input DEM 0.5m', fontsize=18)

crop_sample(dem2x_hires).plot(vmin=1270,vmax=1400,cmap='RdBu_r',ax=ax6)
ax6.set_title('')

for _ax in ax:
    for __ax in _ax:
        __ax.ticklabel_format(useOffset=False, style='plain')
        __ax.set_xlabel('')
        __ax.set_ylabel('')
        
plt.suptitle(title, fontsize=20)
fig.tight_layout(rect=[0.03, 0.03, 0.97, 1])
#plt.savefig('Super-resolution DEM.jpg', dpi=150)
plt.show()
# compare original and target DEM
float((dem2x_hires-dem2x).mean()),float((dem2x_hires-dem2x).std())<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13664278403275332, 'fraction_numerical': 0.060644831115660186, 'mean_word_length': 6.320224719101123, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '348015', 'license': 'MIT', 'n_tokens_mistral': 1700, 'n_tokens_neox': 1542, 'n_words': 287, 'path': 'Elevation1 - Mashhad, Iran/dem-super-resolution.ipynb', 'repo_name': 'mobigroup/gis-snippets', 'size': '3918'}