总集篇！三种山体阴影绘制方法

用户11172986

发布于 2024-06-20 18:08:33

790

发布于 2024-06-20 18:08:33

前言

在地理信息系统（GIS）和地形分析中，山体阴影（也称为地形阴影）是一种重要的可视化技术，它通过模拟太阳光照对地形起伏产生的阴影效果，增强地形的三维感觉，使地图读者能够直观地感受到地形的高低起伏和复杂性。这种技术不仅广泛应用于地质研究、城市规划、环境评估等领域，而且因其所提供的美观、直观的视觉效果，也常见于各类地图和地理信息产品中。

山体阴影的生成通常基于数字高程模型（DEM），即一个二维数组，其中每个元素的值代表对应地理位置的高度。通过模拟特定方向和角度的光照条件，计算每个地理位置的明暗程度，从而生成整个地区的山体阴影图。这一过程涉及到光照模型、视角和高度数据的复杂计算，但幸运的是，许多GIS软件和图形库已经提供了实现这些计算的工具和函数。

在本文中，我们将介绍三种不同的山体阴影绘制方法，每种方法都使用Python编程语言，并借助于流行的Matplotlib库和Cartopy库来实现。这些方法包括：

梯度计算方法：这是最基础的山体阴影绘制方法，通过模拟单一光源（如太阳）照射到地形上，根据地形的高度和坡度计算阴影效果。 matplotlib多模式山体阴影：该方法通过使用不同的混合模式来增强山体阴影的视觉效果，例如使用hsv、overlay和soft等混合模式，可以更细致地展示地形特征和光照变化。xarray-xarrayspatial函数山体阴影：在这种方法中，代码最少。通过掌握这些技术，您将能够为您的地理信息项目或地形分析任务创建更加生动和信息丰富的地形可视化效果。接下来，我们将逐一介绍每种方法的实现步骤和代码示例，帮助您快速上手并应用于实际项目中。

方法一：matplotlib 点灯法

来源：https://matplotlib.org.cn/gallery/specialty_plots/topographic_hillshading.html 关键代码：from matplotlib.colors import LightSource

In [1]:

!pip install cnmaps -i https://pypi.mirrors.ustc.edu.cn/simple/

In [2]:

import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib.colors import LightSource
import cmaps
from cnmaps import get_adm_maps, draw_map
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
from matplotlib.colors import LightSource
import xarray as xr

In [3]:

def load_dem_data(file_path):
    """加载地形高程数据"""
    ds = nc.Dataset(file_path)
    _lon = ds.variables['LON'][:]
    _lat = ds.variables['LAT'][:]
    _dem = ds.variables['elevation'][:]
    return _lon, _lat, _dem

def process_dem_data(_lon, _lat, _dem, lon_range, lat_range):
    """处理地形高程数据"""
    lon = _lon[lon_range[0]:lon_range[1]]
    lat = _lat[lat_range[0]:lat_range[1]]
    dem = _dem[lat_range[0]:lat_range[1], lon_range[0]:lon_range[1]]
    return lon, lat, dem

file_path = '/home/mw/input/china_dem3276/cldasgrid_dem.nc'

lon, lat, dem = load_dem_data(file_path)
lon_range = (4500 ,5500)
lat_range = (1500, 2100)
lon, lat, dem = process_dem_data(lon, lat, dem, lon_range, lat_range)

In [4]:

# 创建光源对象
ls = LightSource(azdeg=315, altdeg=45)

# 设置图形大小
fig = plt.figure(figsize=(30, 15))

# 定义三种混合模式
blend_modes = ['hsv', 'soft', 'overlay']

# 循环遍历每种混合模式并创建子图
for i, blend_mode in enumerate(blend_modes, start=1):
    # 使用shade方法并指定混合模式
    rgb = ls.shade(dem, cmap=plt.cm.gist_earth, blend_mode=blend_mode,
                   vert_exag=5, dx=10, dy=10, fraction=1.5,vmin=-1500)
    
    # 创建子图
    ax = fig.add_subplot(1, 3, i, projection=ccrs.PlateCarree())
    ax.set_title(f"Blend mode: {blend_mode}")
    
    # 绘制地形图
    ax.imshow(rgb, extent=[101, 110, 27, 35], transform=ccrs.PlateCarree(), origin='lower')
plt.show()

效果不满意可调节参数

方法二梯度法

In [5]:

plt.figure(figsize=(15, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
x, y = np.gradient(dem)

slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))

# -x here because of pixel orders in the SRTM tile
aspect = np.arctan2(-x, y)

altitude = np.pi / 4.
azimuth = np.pi / 2.

shaded = np.sin(altitude) * np.sin(slope)\
    + np.cos(altitude) * np.cos(slope)\
    * np.cos((azimuth - np.pi/2.) - aspect)


plt.imshow(shaded,extent=(lon.min(), lon.max(), lat.min(), lat.max()), transform=ccrs.PlateCarree(),
           cmap=plt.cm.gist_earth, origin='lower')  
draw_map(get_adm_maps(province='江苏省', only_polygon=True, record='first'), color='w', linewidth=2)
plt.show()