Python 创建地形图

原始地 DEM。

火山口湖 (OR) 区域的起始 DEM。数据来自 NASA

DEM 本身非常美丽，但我们先进行分层。

将自定义色彩图应用于 DEM

对于我在 ArcGIS Pro 版本中所做的初始高程样式着色，我使用了“高程 #7”。在 matplotlib 中可用的标准颜色图中，我没有看到任何接近此方案的内容，因此我决定基于它创建自定义颜色图。

首先，我截取了 ArcGIS Pro 色彩图的屏幕截图并将其保存为 PNG 格式。

ArcGIS Pro 中海拔 #7 的屏幕截图。

接下来我编写了以下代码来从中提取颜色并返回颜色图。

from PIL import Image
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.pyplot as plt


def extract_colors_from_image(image_path: str) -> np.ndarray:
    """Extract colors from the center line of an image."""
    image = Image.open(image_path).convert('RGB')
    width, height = image.size
    colors = [image.getpixel((x, int(height / 2))) for x in range(width)]
    return np.array(colors) / 255.0


colors = extract_colors_from_image('input/ramp.png')
color_ramp = LinearSegmentedColormap.from_list("custom_ramp", colors)

为了检查它，我们可以使用 matplotlib 显示颜色图。

plt.imshow([colors], aspect='auto')
plt.title('Custom Color Ramp')
plt.axis('off')
plt.show()

自定义地形色彩图。

为了查看该区域的外观，让我们加载 DEM 并预览结果。

from osgeo import gdal
import numpy as np


def normalize_array(array: np.ndarray, lower_percentile: float = 0.0, upper_percentile: float = 100.0) -> np.ndarray:
    """Normalize array values to a specified percentile range."""
    lower_val = np.percentile(array, lower_percentile)
    upper_val = np.percentile(array, upper_percentile)
    clipped = np.clip(array, lower_val, upper_val)
    return (clipped - lower_val) / (upper_val - lower_val)


dem_path = 'input/oblique-clip.tif'
dem_data = gdal.Open(dem_path).ReadAsArray()
normalized_dem = normalize_array(dem_data)

plt.imshow(normalized_dem, cmap=color_ramp)
plt.colorbar()
plt.title('DEM with Custom Color Ramp')
plt.show()

这看起来很棒。它与 ArcGIS Pro 中的外观几乎完全相同。让我们开始吧。

def apply_colormap(array: np.ndarray, cmap: LinearSegmentedColormap) -> np.ndarray:
    """Apply a matplotlib colormap to an array."""
    colormap = plt.get_cmap(cmap) if isinstance(cmap, str) else cmap
    normed_data = (array - array.min()) / (array.max() - array.min())
    colored = colormap(normed_data)
    return (colored[:, :, :3] * 255).astype('uint8')


terrain_dem = apply_colormap(normalized_dem, color_ramp)
terrain_image = Image.fromarray(terrain_dem)
terrain_image.save('output/1-terrain.png')

应用了自定义颜色图的地形。

地形变暖

在 John Nelson 的教程中，他使用 ArcGIS Pro 中的“inferno”色图作为他的变暖层。matplotlib 中的等效色图是“plasma”色图。我们可以使用它的名称将其应用于 DEM。

warm_dem = apply_colormap(normalized_dem, 'plasma')
warm_image = Image.fromarray(warm_dem)
warm_image.save('output/2-warming.png')

血浆中的 DEM。

现在我们得到了两个版本的 DEM：自定义颜色图和等离子版本。在 ArcGIS Pro 中，使用柔光混合模式将它们组合在一起。柔光结合了图层的亮度值并增强了对比度。在 Python 中，我们可以使用 PIL 库来执行此操作。

from PIL import ImageChops

warming_blended = ImageChops.soft_light(terrain_image, warm_image)
warming_blended.save('output/3-warming_blended.png')

使用柔和的光用等离子体加热海拔。

现在是时候去一些山影了。

应用传统山体阴影

我编写了两个函数来处理山体阴影生成。第一个函数生成单个山体阴影，第二个函数使用该函数生成并组合多个山体阴影（下一步将用于多方向山体阴影层）。我不确定这是否接近 ArcGIS Pro 内部发生的事情，但结果非常接近。

def generate_hillshade(dem_path: str, azimuth: float, altitude: float, zFactor: float = 1.0) -> np.ndarray:
    """Generate a hillshade for a specific azimuth and altitude."""
    options = gdal.DEMProcessingOptions(format='GTiff', computeEdges=True, zFactor=zFactor, azimuth=azimuth, altitude=altitude)
    hillshade_path = f'output/temp_hillshade_{azimuth}_{altitude}.tif'
    gdal.DEMProcessing(hillshade_path, dem_path, 'hillshade', options=options)
    return gdal.Open(hillshade_path).ReadAsArray()


def combine_hillshades(dem_path: str, azimuths: List[float], altitudes: List[float], weights: Optional[List[float]] = None) -> np.ndarray:
    """Combine hillshades from multiple directions."""
    hillshades = [generate_hillshade(dem_path, az, alt) for az, alt in zip(azimuths, altitudes)]

    if weights is None:
        weights = np.ones(len(azimuths))

    weights = np.array(weights) / np.sum(weights)
    combined_hillshade = np.average(hillshades, axis=0, weights=weights)
    combined_hillshade = np.clip(combined_hillshade, 0, 255)
    return combined_hillshade

我没有费心清理中间的山体阴影，因为我想在玩输出时检查它们。

以下生成传统（单一方位角和高度）山体阴影，然后使用叠加混合模式将其混合到变暖的地形中。叠加会使暗区变暗，亮区变亮。

def rgb_image_from_array(array: np.ndarray) -> Image.Image:
    """Create an RGB image from an array."""
    return Image.fromarray((array * 255 / array.max()).astype('uint8')).convert("RGB")


traditional_hillshade = combine_hillshades(dem_path, [315], [45])
traditional_hillshade_image = rgb_image_from_array(traditional_hillshade)
traditional_hillshade_image.save('output/4-traditional_hillshade.png')
traditional_hillshade_blended = ImageChops.overlay(warming_blended, traditional_hillshade_image)
traditional_hillshade_blended.save('output/5-hillshade_blended.png')

传统的山体阴影。

使用叠加混合将传统的山体阴影融入到暖色图像中。

看起来不错并且仍然跟踪 ArcGIS Pro 版本。

应用多方向山体阴影

此步骤使用与上一步相同的函数，只是这次我们输入了四个不同的方位角和高度，然后一起取平均值。多方向山体阴影可以更好地模拟现实世界，从而产生更逼真的效果。

azimuths = [45, 135, 225, 315]
altitudes = [45, 45, 45, 45]
multidirectional_hillshade = combine_hillshades(dem_path, azimuths, altitudes)
multidirectional_hillshade_image = rgb_image_from_array(multidirectional_hillshade)
multidirectional_hillshade_image.save('output/6-multidirectional_hillshade.png')
multidirectional_hillshade_blended = ImageChops.multiply(traditional_hillshade_blended, multidirectional_hillshade_image)
multidirectional_hillshade_blended.save('output/7-blended-multidirectional_hillshade.png')

多方向的山体阴影。真漂亮。

使用乘法混合的多方向山体阴影。这是我最喜欢的中间步骤之一。

乘法将混合层的颜色值与基础层的颜色值相乘，从而产生颜色组合后的整体图像较暗。这会稍微放大阴影。

低光山体阴影

最后一个山体阴影是另一种传统的山体阴影，但光线角度较低。这个是用柔和的光线混合而成的。

low_light_hillshade = combine_hillshades(dem_path, [315], [25])
low_light_hillshade_image = rgb_image_from_array(low_light_hillshade)
low_light_hillshade_image.save('output/8-low_light_hillshade.png')
low_light_hillshade_blended = ImageChops.soft_light(multidirectional_hillshade_blended, low_light_hillshade_image)
low_light_hillshade_blended.save('output/9-low_light_hillshade_blended.png')

低光山体阴影。

低光山体阴影混合。

灯光

下一层是照明层，用于使山峰变亮并使山谷变暗。为此，我们可以使用 matplotlib 中现有的白色到黑色（二进制反转）颜色图。

lighting_dem = apply_colormap(normalized_dem, 'binary_r')
lighting_image = Image.fromarray(lighting_dem)
lighting_image.save('output/10-lighting.png')
lighting_blended = ImageChops.soft_light(low_light_hillshade_blended, lighting_image)
lighting_blended.save('output/11-lighting_blended.png')

照明层。黑暗的山谷，明亮的山峰。就像在现实生活中一样。

灯光与柔和的灯光融为一体。

雾时！

接下来的两个步骤是我最喜欢的。首先，喷雾。

雾层需要另一个自定义颜色图。它需要在低海拔处呈半透明白色，在海拔值过高之前迅速转变为完全透明的白色。要获得恰到好处的过渡需要进行一些实验。

colors = [
    (1.0, 1.0, 1.0, 0.85),
    (1.0, 1.0, 1.0, 0.45),
    (1.0, 1.0, 1.0, 0.15),
    (1.0, 1.0, 1.0, 0.0),
    (1.0, 1.0, 1.0, 0.0)
]
positions = [0.0, 0.15, 0.25, .35, 1.0]
mist_custom_cmap = LinearSegmentedColormap.from_list("custom_cmap", list(zip(positions, colors)))

colors 列表中的每个元素都是一个元组，包含 RGBA（红色、绿色、蓝色、Alpha）颜色分量的值。所有颜色都有三个 1，即白色。唯一的变化是 alpha 分量，它决定透明度。第一个记录 0.85 表示 85% 不透明（轻微透明度）。

位置列表指定了颜色列表中每种颜色在颜色图中的确切位置，范围从 0 到 1。每个位置值对应于渐变中的特定点，确定颜色之间的过渡发生的位置。在这个自定义颜色图中，我们逐渐将白色的不透明度从开始时的 85% 更改为完全透明，每次更改 35%，当我们从低海拔到高海拔范围的三分之一时，它实际上就完全透明了。

让我们应用它。

mist_dem = mist_custom_cmap(normalized_dem)
mist_image = Image.fromarray((mist_dem * 255).astype(np.uint8), 'RGBA')
mist_image.save("output/12-dem_with_transparent_mist.png")

terrain_image = Image.open('output/11-lighting_blended.png').convert('RGBA')
combined_terrain_and_mist = Image.alpha_composite(terrain_image, mist_image).convert('RGB')
combined_terrain_and_mist.save('output/13-combined_terrain_and_mist.png')

黑暗背景上弥漫着薄雾。

地形上空雾气弥漫。

最后…山体阴影的坡度

最后一步，我们生成坡度栅格并应用“cividis”色彩图来突出显示陡峭区域。

def generate_slope_array(dem_path: str) -> Tuple[np.ndarray, Optional[float]]:
    """Generate a slope array from a DEM."""
    dem_data = gdal.Open(dem_path, gdal.GA_ReadOnly)
    slope_ds = gdal.DEMProcessing('', dem_data, 'slope', format='MEM', scale=1)
    slope_band = slope_ds.GetRasterBand(1)
    slope_array = slope_band.ReadAsArray()
    no_data_value = slope_band.GetNoDataValue()
    
    if no_data_value is not None:
        mask = slope_array == no_data_value
        slope_array = np.ma.masked_where(mask, slope_array)

    min_val = np.min(slope_array)
    max_val = np.max(slope_array)
    norm_slope_array = (slope_array - min_val) / (max_val - min_val)
    norm_slope_array = np.ma.filled(norm_slope_array, 0)
    
    return norm_slope_array

cividis_cmap = plt.get_cmap('cividis')
slope_dem = cividis_cmap(generate_slope_array(dem_path))
slope_image = rgb_image_from_array(slope_dem)
slope_image.save('output/14-slope.png', 'png')

应用了 cividis 色彩图的斜率。

最终结果

最后一步是用柔和的光线将斜坡融入雾蒙蒙的地形。

slope_blended = ImageChops.soft_light(combined_terrain_and_mist, slope_image)
slope_blended.save('output/15-slope_blended.png')

最终结果。

我对结果非常满意。仅使用几个 Python 库，我就能获得非常接近 ArcGIS Pro 版本的结果。

以下是一些特写镜头。

雾霭！金脊从坡层而来。

高地。