气象编程 | Python基础地图构建

import numpy as np 
import matplotlib.pyplot as plt 
from mpl_toolkits.basemap import Basemap
import cmaps
import shapefile
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import os
import maskout2
from netCDF4 import Dataset #netCDFを扱うため
nc = Dataset('F:/Rpython/lp14/MERRA2_400.tavgM_2d_slv_Nx.201912.nc4.nc4', mode='r')
lon = nc.variables["lon"][:].data
lat = nc.variables["lat"][:].data
t2m0= nc.variables["T2M"][:,:,:].data
t2m = t2m0.mean(axis=0)
grid_z = t2m -273.15
fig=plt.figure(figsize=(16,9))
ax=fig.add_subplot(111)
mp=Basemap(llcrnrlon=70,llcrnrlat=15,urcrnrlon=140,urcrnrlat=55,projection='cyl')
mp.drawparallels(np.arange(15, 55 + 1, 10), labels = [1, 0, 0, 0],fontsize=14,linewidth='0.2',color='black')
mp.drawmeridians(np.arange(70, 140 + 1, 10), labels = [0, 0, 0, 1],fontsize=14,linewidth='0.2',color='black')
CHN='F:/Rpython/lp27/data/china-shapefiles-master/'
mp.readshapefile(CHN+'china_nine_dotted_line','china_nine_dotted_line',drawbounds=True,linewidth=1.5)  
mp.readshapefile(CHN+'china','china',drawbounds=True)
clevs = np.arange(-50,50,10)
cf=plt.contourf(lon,lat,grid_z,clevs,cmap='Spectral_r',extend='both')
cbar=mp.colorbar(cf,location='right',size=0.3)
clip=maskout2.shp2clip(cf,ax,r'F:/Rpython/lp27/data/china-shapefiles-master/china_country')
#添加南海，实际上就是新建一个子图覆盖在之前子图的右下角
sub_ax = fig.add_axes([0.6888,0.1183,0.25,0.25])    #[*left*, *bottom*, *width*,*height*]
mf=Basemap(llcrnrlon=107,llcrnrlat=2,urcrnrlon=122,urcrnrlat=22,projection='cyl')
clevs2 = np.arange(-50,50,10)
cf2=sub_ax.contourf(lon, lat, grid_z,clevs2,cmap='Spectral_r')
CHN='F:/Rpython/lp27/data/china-shapefiles-master/'
mf.readshapefile(CHN+'china','china',drawbounds=True)
CHN='F:/Rpython/lp27/data/china-shapefiles-master/'
mf.readshapefile(CHN+'china_nine_dotted_line','china_nine_dotted_line',drawbounds=True,linewidth=1.5)
clip=maskout2.shp2clip(cf2,sub_ax,r'F:/Rpython/lp27/data/china-shapefiles-master/china_country')
plt.savefig(r"F:/Rpython/lp27/plot42.3.png",dpi=600)
plt.show()

import maskout2
import os
import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeat
import matplotlib.colors as mcolors
from cartopy.io.shapereader import Reader
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
from matplotlib import rcParams
config = {"font.family":'Times New Roman', "font.size": 16, "mathtext.fontset":'stix'}
rcParams.update(config)
SHP = 'F:/Rpython/lp12/es'
# 数据读取
ds = xr.open_dataset(r'F:/Rpython/lp27/data/pr_Amon_FGOALS-f3-L_ssp370_r1i1p1f1_gr_201501-210012.nc')
pr = ds['pr']
#取一段时间内的降水作为绘图的数据yuan
pr_0=pr[793:1032]
pr_ave=pr_0.mean(dim='time')*86400*30*12
lon=ds.lon
lat=ds.lat
lon_range=lon[(lon>70)&(lon<140)]
lat_range=lat[(lat>0)&(lat<60)]
pr_region=pr_ave.sel(lon=lon_range,lat=lat_range)
# 创建画图空间
proj = ccrs.PlateCarree()  #创建投影
fig = plt.figure(figsize=(16,12))  #创建页面
ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  #主图
# 设置网格点属性
region=[70, 140, 0, 60]
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.grid(linestyle = '--')
ax.grid(linewidth=1.2,color='r',alpha=0.2,linestyle='--')
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
ax.yaxis.set_major_formatter(LatitudeFormatter())
# 画图
levels = np.arange(0,3000,500)
font2={'family':'SimHei','size':20,'color':'k'}
ax.set_title('2080-2100年均降水总量unit:mm',loc='left',fontdict=font2)
cs=ax.contourf(lon_range,lat_range,pr_region,levels=levels,cmap='Spectral_r',transform=ccrs.PlateCarree(),extend='both')
b=plt.colorbar(cs,shrink=0.88,orientation='vertical',extend='both',pad=0.035,aspect=20)
#白化
ax.add_geometries(Reader(os.path.join(SHP,'bou2_4l.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1)
clip=maskout2.shp2clip(cs,ax,r'F:/Rpython/lp3/hls/china0')
#添加南海，实际上就是新建一个子图覆盖在之前子图的右下角
f2_ax2=fig.add_axes([0.6566,0.133,0.1,0.17],projection=proj)
f2_ax2.set_extent([105,125,0,25],crs=ccrs.PlateCarree())
f2_ax2.add_feature(cfeature.COASTLINE.with_scale('50m'))
china=shpreader.Reader('F:/Rpython/lp12/es/bou2_4l.dbf').geometries()
f2_ax2.add_geometries(china,ccrs.PlateCarree(),facecolor='none',edgecolor='black',zorder=1)
cs3=f2_ax2.contourf(lon_range,lat_range,pr_region,range(0,3000,500),cmap='Spectral_r')
#白化
f2_ax2.add_geometries(Reader(os.path.join(SHP,'bou2_4l.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1)
clip=maskout2.shp2clip(cs3,f2_ax2,r'F:/Rpython/lp3/hls/china0')
plt.savefig('F:/Rpython/lp27/plot8.7.png',dpi=300)
plt.show()#显示图片