本文作者：謝金材，南京信息工程大學(xué)碩士
聯(lián)系郵箱：jincaixie@163.com
# -*- coding: utf-8 -*-'''Created on Sat Apr 15 16:02:47 2023@author: 31015'''import numpy as npimport xarray as xrimport matplotlib.pyplot as pltimport scipyimport cartopy.crs as ccrsimport cartopy.mpl.ticker as ctickerimport matplotlib.ticker as mtickerimport cartopy.io.shapereader as shpreaderimport maskoutfilename=r'H:/NCEP/monthly/hgt.mon.mean.nc'ds=xr.open_dataset(filename)z3=np.array(ds['hgt'].loc[ds.time.dt.month.isin([3])].loc['1979-01-01':'2022-12-31',500,:,:],dtype='float32').squeeze()z4=np.array(ds['hgt'].loc[ds.time.dt.month.isin([4])].loc['1979-01-01':'2022-12-31',500,:,:],dtype='float32').squeeze()lon=np.array(ds.lon,dtype='float32')lat=np.array(ds.lat,dtype='float32')r=np.empty((73,144))p=np.empty((73,144))for i in range(len(lat)): for j in range((len(lon))): r[i,j],p[i,j]=scipy.stats.pearsonr(z3[:,i,j],z4[:,i,j]) fig,ax1 = plt.subplots(figsize=(6,5), subplot_kw={'projection':ccrs.PlateCarree()})leftlon, rightlon, lowerlat, upperlat = (70,135,10,60)img_extent = [leftlon, rightlon, lowerlat, upperlat]ax1.set_extent(img_extent, crs=ccrs.PlateCarree())china = shpreader.Reader(r'H:\prediction\map\bou2_4l.dbf').geometries()ax1.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='gray',zorder = 1,lw=0.5)country=shpreader.Reader(r'H:\prediction\map\country1.dbf').geometries()ax1.add_geometries(country, ccrs.PlateCarree(),facecolor='none', edgecolor='gray',zorder = 1,lw=0.5)ax1.set_title('title',loc='left',fontsize =15)ax1.set_xticks(np.arange(leftlon,rightlon+1,10))ax1.set_yticks(np.arange(lowerlat,upperlat+1,10))lon_formatter = cticker.LongitudeFormatter()lat_formatter = cticker.LatitudeFormatter()ax1.xaxis.set_major_formatter(lon_formatter)ax1.yaxis.set_major_formatter(lat_formatter)ax1.xaxis.set_minor_locator(mticker.MultipleLocator(5))#顯示x軸副刻度ax1.yaxis.set_minor_locator(mticker.MultipleLocator(5))#顯示x軸副刻度plt.tick_params(labelsize=15)cs = ax1.contourf(lon,lat,r, zorder=0 ,levels=np.arange(-1,1.1,0.1) \ , extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.bwr)clip=maskout.shp2clip(cs,ax1,'country1','China')# ax1.scatter(lon,lat,s =2,color='black',transform=ccrs.PlateCarree())cf=ax1.contourf(lon,lat,p, [0,0.05,1],zorder=2,hatches=[ 'xxx',None],colors='none',)for collection in cf.collections: collection.set_edgecolor('green') collection.set_linewidth(0)clip=maskout.shp2clip(cf,ax1,'country1','China')axx=ax1.inset_axes([0,0,1,0.3],anchor=(1,0),projection=ccrs.PlateCarree())leftlon, rightlon, lowerlat, upperlat = (107,122,0,22)extend_southsea=[leftlon, rightlon, lowerlat, upperlat]axx.set_extent(extend_southsea, crs=ccrs.PlateCarree())china = shpreader.Reader(r'H:\prediction\map\bou2_4l.dbf').geometries()axx.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='gray',zorder = 1,lw=0.5)country=shpreader.Reader(r'H:\prediction\map\country1.dbf').geometries()axx.add_geometries(country, ccrs.PlateCarree(),facecolor='none', edgecolor='gray',zorder = 1,lw=0.5)css=axx.contourf(lon,lat,r, zorder=0 ,levels=np.arange(-1,1.1,0.1) \ , extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.bwr)clip=maskout.shp2clip2(css,axx,'country1','China')csf=axx.contourf(lon,lat,p, [0,0.05,1],zorder=2,hatches=['xxx',None],colors='none')for collection in csf.collections: collection.set_edgecolor('green') collection.set_linewidth(0)clip=maskout.shp2clip2(csf,axx,'country1','China')position=fig.add_axes([0.25, 0.05, 0.5, 0.03])cb=fig.colorbar(cs,cax=position,orientation='horizontal',format='%.1f',ticks=np.arange(-1,1.1,0.4),extendrect=True,extendfrac='auto')cb.ax.tick_params(labelsize=15)plt.savefig(r'test-ori.png',dpi=500)plt.show()
#coding=utf-8########################################################################################################################################This module enables you to maskout the unneccessary data outside the interest region on a matplotlib-plotted output instance####################in an effecient way,You can use this script for free ##################################################################################################################################################################################################USAGE: INPUT include 'originfig':the matplotlib instance### 'ax': the Axes instance# 'shapefile': the shape file used for generating a basemap A# 'region':the name of a region of on the basemap A,outside the region the data is to be maskout# OUTPUT is 'clip' :the the masked-out or clipped matplotlib instance.import shapefilefrom matplotlib.path import Pathfrom matplotlib.patches import PathPatchfrom shapely.geometry import Point as ShapelyPointfrom shapely.geometry import Polygon as ShapelyPolygonfrom collections import Iterabledef shp2clip(originfig,ax,m,shpfile,region, clabel = None, vcplot = None):sf = shapefile.Reader(shpfile)vertices = []codes = []for shape_rec in sf.shapeRecords():####這里需要找到和region匹配的唯一標識符，record[]中必有一項是對應(yīng)的。if shape_rec.record[3] == region: #####在country1.shp上，對中國以外的其他國家或地區(qū)進行maskout# if shape_rec.record[7] in region: #####在bou2_4p.shp上，對中國的某幾個省份或地區(qū)之外的部分進行maskoutpts = shape_rec.shape.pointsprt = list(shape_rec.shape.parts) + [len(pts)]for i in range(len(prt) - 1):for j in range(prt[i], prt[i+1]):vertices.append(m(pts[j][0], pts[j][1]))codes += [Path.MOVETO]codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)codes += [Path.CLOSEPOLY]clip = Path(vertices, codes)clip = PathPatch(clip, transform=ax.transData)if vcplot:if isinstance(originfig,Iterable):for ivec in originfig:ivec.set_clip_path(clip)else:originfig.set_clip_path(clip)else:for contour in originfig.collections:contour.set_clip_path(clip)if clabel:clip_map_shapely = ShapelyPolygon(vertices)for text_object in clabel:if not clip_map_shapely.contains(ShapelyPoint(text_object.get_position())):text_object.set_visible(False)return clip# 改進代碼def shp2clip2(originfig,ax,shpfile,region,clabel = None, vcplot = None):sf = shapefile.Reader(shpfile)vertices = []codes = []for shape_rec in sf.shapeRecords():if shape_rec.record[3] == region: ####這里需要找到和region匹配的唯一標識符，record[]中必有一項是對應(yīng)的。pts = shape_rec.shape.pointsprt = list(shape_rec.shape.parts) + [len(pts)]for i in range(len(prt) - 1):for j in range(prt[i], prt[i+1]):if pts[j][0]<107 or pts[j][0]>122 or pts[j][1]>22 or pts[j][1]<0:vertices.append((107, 22))else:vertices.append((pts[j][0], pts[j][1]))codes += [Path.MOVETO]codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)codes += [Path.CLOSEPOLY]clip = Path(vertices, codes)# print(clip)clip = PathPatch(clip, transform=ax.transData)for contour in originfig.collections:contour.set_clip_path(clip)if vcplot:if isinstance(originfig,Iterable):for ivec in originfig:ivec.set_clip_path(clip)else:originfig.set_clip_path(clip)else:for contour in originfig.collections:contour.set_clip_path(clip)if clabel:clip_map_shapely = ShapelyPolygon(vertices)for text_object in clabel:if not clip_map_shapely.contains(ShapelyPoint(text_object.get_position())):text_object.set_visible(False)return clip
獲取本文Python白化南海小地圖數(shù)據(jù)和代碼的途徑，
