Python 绘制山体阴影+雷达图

用户11172986

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

730

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

文章被收录于专栏：气python风雨气python风雨

前言

本文旨在利用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[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

In [3]:

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