Python 绘制山体阴影+雷达图
本帖最后由 气ython风雨 于 2024-3-6 20:04 编辑前言
本文旨在利用Python编程语言，将山体阴影与雷达速度产品相结合，以探索其可视化效果
环境：python 3.9
导入库
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
import cinrad
from cinrad.io import CinradReader, StandardData
数据处理
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]
lat = _lat:lat_range]
dem = _dem:lat_range, lon_range:lon_range]
return lon, lat, dem
In :
file_path = '/home/mw/input/china_dem3276/cldasgrid_dem.nc'
f = CinradReader("/home/mw/input/data5692/Z_RADR_I_Z9250_20200612054800_O_DOR_SA_CAP.bin")
V = f.get_data(5, 230, 'VEL')
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)
绘图一：不加阴影
plt.figure(figsize=(10, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
plt.imshow(dem,extent=(lon.min(), lon.max(), lat.min(), lat.max()), transform=ccrs.PlateCarree(),
cmap=cmaps.MPL_terrain, origin='lower')
plt.contourf(V.longitude, V.latitude,V.VEL , levels=np.arange(0, 50, 5), cmap=cmaps.radar,
transform=ccrs.PlateCarree())
draw_map(get_adm_maps(province='江苏省', only_polygon=True, record='first'), color='w', linewidth=2)
plt.show()
绘图二：加阴影
plt.figure(figsize=(10, 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='Greys', origin='lower')
plt.contourf(V.longitude, V.latitude,V.VEL , levels=np.arange(0, 50, 5), cmap=cmaps.radar,
transform=ccrs.PlateCarree())
draw_map(get_adm_maps(province='江苏省', only_polygon=True, record='first'), color='w', linewidth=2)
plt.show()
完整的代码与文件
文章来源于微信公众号：气ython风雨
页:
[1]