import numpy as npimport matplotlib.pyplot as pltimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.io.shapereader import Readerfrom cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatterimport xarray as xrimport rioxarray as rxrimport maskoutfrom matplotlib import rcParamsconfig = {"font.family":'Times New Roman',"font.size": 12,"mathtext.fontset":'stix'}rcParams.update(config)region=[70, 150, 0, 60]filename=r'./ETOPO2022.tif'ds = rxr.open_rasterio(filename)lons,lats = np.meshgrid(ds['x'],ds['y'])dem = ds[0]# plotfig=plt.figure(dpi=600,figsize=(16,9))proj=ccrs.PlateCarree()ax=fig.add_axes([0,0,0.4,0.4],projection=proj)ax.set_extent(region, crs = proj)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_feature(cfeature.OCEAN.with_scale('50m'))ax.add_feature(cfeature.LAND.with_scale('50m'))ax.add_feature(cfeature.LAKES.with_scale('50m'))ax.add_geometries(Reader(r'./river1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='b',linewidth=1)ax.add_geometries(Reader(r'./ming.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1)ax.add_geometries(Reader(r'./henan_province.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1)ax.add_geometries(Reader(r'./nineline.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1)ax.add_geometries(Reader(r'./ne_10m_land.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',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'./henan1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1.2,zorder=1)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')c11=ax.contourf(lons,lats,dem,np.arange(0,8000,1000),extend='both',transform=ccrs.PlateCarree(),cmap='gist_rainbow',zorder=0)maskout(c11,ax,'./ming.shp',1)font3={'family':'SimHei','size':12,'color':'k'}plt.ylabel("纬度",fontdict=font3)plt.xlabel("经度",fontdict=font3)position=fig.add_axes([0.08,0.02,0.005,0.12])#位置[左,下,右,上]cbar=plt.colorbar(c11,cax=position,shrink=0.4,aspect=100,fraction=.03,pad=0.01) #aspect控制bar宽度,fraction控制大小比例,pad控制与图的距离cbar.set_ticks(np.arange(0,8000,1000)) #设置colorbar范围和刻度标记间隔cbar.ax.tick_params(labelsize=6,direction='in',right=False)font3={'family':'SimHei','size':12,'color':'k'}ax.set_title('明朝DEM地形空间分布图',fontdict=font3)ax2=fig.add_axes([0.3,0,0.4,0.4],projection=proj)region2=[110,117,31,37]ax2.set_extent(region2, crs = proj)# ax.coastlines(scale, linewidth = 0.5)ax2.set_xticks(np.arange(region2[0], region2[1] + 1, 1), crs = proj)ax2.set_yticks(np.arange(region2[-2], region2[-1] + 1, 1), crs = proj)ax2.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))ax2.yaxis.set_major_formatter(LatitudeFormatter())ax2.add_feature(cfeature.OCEAN.with_scale('50m'))ax2.add_feature(cfeature.LAND.with_scale('50m'))ax2.add_feature(cfeature.LAKES.with_scale('50m'))ax2.add_geometries(Reader(r'./river1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='b',linewidth=1)ax2.add_geometries(Reader(r'./ming.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1)ax2.add_geometries(Reader(r'./henan_province.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1)ax2.add_geometries(Reader(r'./nineline.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1)ax2.add_geometries(Reader(r'./ne_10m_land.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)ax2.add_geometries(Reader(r'./ne_50m_lakes.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)ax2.add_geometries(Reader(r'./henan1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1.2,zorder=1)ax2.add_geometries(Reader(r'./river1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='b',linewidth=0.2)ax2.add_geometries(Reader(r'./ne_50m_lakes.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)ax2.add_geometries(Reader(r'./1级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='RoyalBlue',linewidth=0.4)ax2.add_geometries(Reader(r'./2级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='DodgerBlue',linewidth=0.3)ax2.add_geometries(Reader(r'./3级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='DeepSkyBlue',linewidth=0.2)ax2.add_geometries(Reader(r'./4级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='SkyBlue',linewidth=0.15)ax2.add_geometries(Reader(r'./5级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='LightSkyBlue',linewidth=0.05)ax2.add_geometries(Reader(r'./主要湖泊.shp').geometries(),ccrs.PlateCarree(),edgecolor='none',linewidth=0,facecolor='#BEE8FF')c12=ax2.contourf(lons,lats,dem,np.arange(0,1800,200),extend='both',transform=ccrs.PlateCarree(),cmap='gist_rainbow',zorder=0)maskout(c12,ax2,'./henan_province.shp',1)cbar=plt.colorbar(c12,shrink=0.98,aspect=30,fraction=.03,pad=0.01) #aspect控制bar宽度,fraction控制大小比例,pad控制与图的距离cbar.set_ticks(np.arange(0,1800,200)) #设置colorbar范围和刻度标记间隔cbar.ax.tick_params(labelsize=12, direction='in', right=False)font3={'family':'SimHei','size':12,'color':'k'}ax2.set_title('河南DEM地形空间分布图',fontdict=font3)font3={'family':'SimHei','size':12,'color':'k'}plt.ylabel("纬度",fontdict=font3)plt.xlabel("经度",fontdict=font3)ax.plot([110,117,117,110,110],[31,31,37,37,31],lw=1,transform=ccrs.PlateCarree(),color='k')import matplotlib.patches as mpatchestrans_proj=ccrs.PlateCarree()#用于转换的原坐标target_proj=ccrs.PlateCarree()x=np.array([117,117])y=np.array([37,31])x2=np.array([110,110])for m,n,h,j in zip(x,y,x2,y):#由于移动了中心坐标轴,此处的x转化需减去这个值 xyA=target_proj.transform_point(m,n,trans_proj) xyB=[h,j] con=mpatches.ConnectionPatch( xyA=xyA,xyB=xyB,coordsA=ax.transData,coordsB=ax2.transData,color='k',ls='-') fig.add_artist(con)plt.savefig('./plot280.1.png',dpi=800,bbox_inches='tight',pad_inches=0)plt.show()

c11=ax.contourf(lons,lats,dem,np.arange(0,8000,1000),extend='both',transform=ccrs.PlateCarree(),cmap=cmaps.MPL_terrain,zorder=0)
