绘制science研究区域概况图

Python绘制研究区域概况图：

import numpy as np
import matplotlib.pyplot as plt
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
from matplotlib import pyplot as plt
from matplotlib.colors import LightSource, ListedColormap
from mpl_toolkits.axes_grid1 import make_axes_locatable
from osgeo import gdal
from osgeo import osr
from matplotlib import rcParams
config = {"font.family":'Times New Roman',"font.size":18,"mathtext.fontset":'stix'}
rcParams.update(config)
# 设置DEM路径
dem_file = "./china_dem.tif"
# 使用GDAL驱动并读取DEM数据
ds = gdal.Open(dem_file)
tif_proj =ds.GetProjection()
tif_geotrans =ds.GetGeoTransform()
#确保TIFF的坐标系统信息
tif_crs = osr.SpatialReference(wkt=tif_proj)
# 读取TIFF数据到numpy数组
band = ds.GetRasterBand(1)
data = band.ReadAsArray()
# 设置无效值
nodata = ds.GetRasterBand(1).GetNoDataValue()
if np.any(data == nodata):
    data = np.ma.masked_equal(data, nodata)
ls=LightSource(azdeg=360,altdeg=30)
# 自定义色带
colors = ["#2B6743","#80A363","#F0EF9B","#C79F4B","#923100","#692610", "#470900"]
cmap = ListedColormap(colors)
rgb=ls.shade(data,cmap=cmap,blend_mode='overlay',vert_exag=2,dx=10,dy=10,fraction=1.05)
region=[70, 140, 15, 55]
fig = plt.figure(figsize=(16,9))
ax = fig.add_subplot(111,projection=ccrs.PlateCarree())
proj=ccrs.PlateCarree()
ax.set_extent(region, crs = proj)
ax.stock_img()
extent=(tif_geotrans[0],tif_geotrans[0]+tif_geotrans[1]*data.shape[1],tif_geotrans[3]+tif_geotrans[5]*data.shape[0],tif_geotrans[3])
img = ax.imshow(rgb,aspect='auto',origin='upper',cmap=cmap,extent=extent,transform=ccrs.PlateCarree())
# 设置colorbar
img2 = ax.imshow(data,cmap=cmap)
cbar=plt.colorbar(img2,shrink=0.85,orientation='vertical',extend='both',pad=0.01,aspect=30)
cbar.set_label('DEM/m')
cbar.ax.tick_params(labelsize=12, direction='in', right=False)
ax.set_xticks(np.arange(region[0], region[1] + 1, 10), crs = proj)
ax.set_yticks(np.arange(region[-2], region[-1] + 1, 10), crs = proj)
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.add_geometries(Reader(r'./China_Province.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='black',linewidth=0.8)
ax.add_geometries(Reader(r'./10line.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='black',linewidth=0.8)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'),linewidth=0.5,zorder=2,color='k')# 添加海岸线
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_geometries(Reader(r'./river1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='b',linewidth=0.2)
ax.add_geometries(Reader(r'./ne_50m_lakes.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)
ax.add_geometries(Reader(r'./1级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='RoyalBlue',linewidth=0.4)
ax.add_geometries(Reader(r'./2级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='DodgerBlue',linewidth=0.3)
ax.add_geometries(Reader(r'./3级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='DeepSkyBlue',linewidth=0.2)
ax.add_geometries(Reader(r'./4级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='SkyBlue',linewidth=0.15)
ax.add_geometries(Reader(r'./5级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='LightSkyBlue',linewidth=0.05)
ax.add_geometries(Reader(r'./主要湖泊.shp').geometries(),ccrs.PlateCarree(),edgecolor='none',linewidth=0,facecolor='#BEE8FF')
ax.set_title('中国立体DEM地形图【审图号GS(2024)0650号】',{'family':'simhei','size':18,'color':'k'})
plt.savefig('./plot330.8.2.png',dpi=200,bbox_inches='tight',pad_inches=0)
plt.show()