pytho地理绘图

前言

• matplotlib 绘图确实方便，但是具体到每一个细节调整的时候是在太繁琐了，所以我选择了另外一个包：proplot进行后续绘图

pip install proplot

矢量裁剪

• 这部分代码来源找不到了，记得是气象家园上面的一个帖子

def shp2clip(originfig, ax, shpfile, fieldVals):
    """
    This method enables you to maskout the unneccessary data
                                            outside the interest region
    :param ax:  the Axes instance
    :param shpfile:  the shape file used for clip
    :param fieldVals:  thi features attributes value list in shape file,
                    outside the region the data is to be masked out
    :return:
    """
    sf = shapefile.Reader(shpfile)
    vertices = []
    codes = []
    for shape_rec in sf.shapeRecords():
        if shape_rec.record[0] in fieldVals:  # 注意这里需要指定你的字段的索引号，我的是第3个字段
            pts = shape_rec.shape.points
            prt = 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((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)

    for contour in originfig.collections:
        contour.set_clip_path(clip)
    return clip

绘图

包的引入

-cmaps需要自行安装，这个包参考的是NCL的配色

import os
import glob
import shapefile
import glob
import gc
import cmaps 
import proplot as plot
import cartopy.crs as ccrs
from cartopy.io.shapereader import Reader
import cartopy.mpl.ticker as mticker
from osgeo import gdal,osr
import numpy as np
from tqdm import *
from matplotlib import pyplot as plt
import matplotlib.pyplot as plt
import cartopy.io.shapereader as shapereader
from matplotlib import rcParams
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import cartopy.io.shapereader as shapereader

定义tif的范围

lons=np.expand_dims(np.arange(112,120,0.05),axis=0).squeeze()
lats=np.expand_dims(np.arange(43,35,-0.05),axis=1).squeeze()

定义坐标

proj1= ccrs.PlateCarree()

其他的一些参数

plot.register_fonts(default=True)
lons=np.expand_dims(np.arange(112,120,0.05),axis=0).squeeze()
lats=np.expand_dims(np.arange(43,35,-0.05),axis=1).squeeze()
proj1= ccrs.PlateCarree()
levels=[0,15,35,55,75,115,160]  #colorbar分级配色的参考
config = {"font.family":'Times New Roman',"font.size":16,"mathtext.fontset":'stix'} #字体及大小定义
rcParams.update(config)

绘图

• 下面cliped中的三个数字是根据你选择的区域进行确定的，为区域对应的abcode
  

fig, axs = plot.subplots(nrows=3,ncols=4,space=0,proj=proj1,figsize=(6,4))#journal 'nat1'
for i in range(11):
    ds=gdal.Open(files[i])
    arr=ds.ReadAsArray()
    t2=axs[i].contourf(lons,lats,arr,transform=ccrs.PlateCarree(),levels=levels,cmap=cmaps.BlAqGrYeOrReVi200)
    jjj = shapereader.Reader(r'京津冀.shp').geometries()
    axs[i].add_geometries(jjj, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
    cliped = shp2clip(t2, axs[i], r'京津冀.shp',['110000','120000','130000']) #根据矢量范围进行裁剪
    axs[i].format(latlim=(35,43), lonlim=(112,120),lonlines=2, latlines=2,)   #定义数据框的范围
    name='{}:00'.format(files[i].split(" ")[-1][:-4])

    axs[i].text(117,36,name,fontsize=8)
fig.colorbar(t2,shrink=0.85,orientation='vertical',extend='both',pad=0.25,aspect=30)
plt.delaxes(axs[2,3]) #删除最后一个图
plt.savefig(savepath,dpi=600,bbox_inches='tight')

效果

cmaps

• 输入下面的命令，查看不同的colorbar

help(cmaps)