Python画地形图（以新疆地区为例）

import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import my_class
filename = r"D:\lunwen\lunwenshuju\DEM\cldasgrid_dem.nc"
f = xr.open_dataset(filename)
dem = f["elevation"]
lon = f["LON"]
lat = f["LAT"]
# 筛选区域
dim1 = (lon>=75)&(lon <= 96)
dim2 = (lat>=35)&(lat<=51)
dem = dem[dim2, dim1]
min_dem = np.min(dem)
max_dem = np.max(dem)
var = (max_dem-min_dem)/10
lon = lon[dim1]
lat = lat[dim2]
lon, lat = np.meshgrid(lon, lat)
plt.rcParams.update({'font.size': 13})
#  解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['KaiTi']
#  解决负号乱码问题
plt.rcParams['axes.unicode_minus'] = False
fig = plt.figure(figsize=[15, 7])
ax1 = fig.add_subplot(111, projection=ccrs.PlateCarree())
ax1.set_xticks(np.arange(75, 95 +2.5, 2.5), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(35, 55 + 2.5, 2.5), crs=ccrs.PlateCarree())
xtick_str = ['75', ' ', '80', ' ', '85', ' ', '90', ' ', '95']
ytick_str = ['35', ' ', '40', ' ', '45', ' ', '50', ' ', '55']
ax1.set_yticklabels(ytick_str, fontsize=14)
ax1.set_xticklabels(xtick_str, fontsize=14)
my_class.readshapefile(r"D:\lunwen\lunwenshuju\Shp\china_basic_map\bou2_4l", linewidth=1, ax=ax1)
my_class.readshapefile(r'D:\lunwen\lunwenshuju\Shp\vn_shp\VNM_adm0', linewidth=1, ax=ax1)
levels = np.arange(min_dem+var, max_dem-var, var)
DEM_tu = ax1.contourf(lon, lat, dem, cmap='RdYlGn', extend='both', levels=levels)
ax1.set_title("新疆塔克拉玛沙漠", fontsize=20)
T_cb = fig.colorbar(DEM_tu, orientation='vertical')
plt.show()

画出来的图片如下：

之前代码是用Basemap进行画图的，原来库不更新了，还复杂。

用cartopy画图很简单，因为DEM也就是个nc数据，只是变量变成DEM了很多卫星数据貌似也有高度数据。

这篇代码只是想分享下python也可以画地形数据，很好看有纹理，如果想知道Basemap怎么画也可私聊。

my_class库百度云盘自己提取吧，放到自己代码文件夹下，调用即可。

my_class

dem数据我每记错的话应该是这个ETOPO高程资料，很早之前下的。

根据自己需求下载合适的分辨率dem数据，地址下载链接：

https://www.ncei.noaa.gov/products/etopo-global-relief-model

shp文件大家自行下载,以后可能会更新一些，百度网盘链接：

shp文件