python使用shp和grib2文件绘制广东省地区降水图（有行政边界）

小白学习python绘图的第一篇文章，不足之处敬请指出。

关于绘制行政区域，我尝试了多种方法，有pyecharts，json等等，最后发现用shp文件比较方便。

首先，我下载了广东省行政区划矢量数据的shp文件，其缩略图如下图所示：

首先，我们了解shp文件的构成，可用以下代码读取内容

file = shapefile.Reader(r"C:\Users\20719\Desktop\Guangdong-city-2020.shp",encoding='GB2312')
# 通过创建reader类的对象进行shapefile文件的读取,encoding避免utf-8报错
records = file.records()

 读取后输出records，结果如下

 

可以看到shp文件中包含多个record，即多个地级市，如果直接读取shp文件中的x，y坐标并使用matplotlib绘制，将会只出现一个地级市，因此需要用循环绘制出所有地级市，代码如下：

sf = shapereader.Reader(shp_path =r'C:\Users\20719\Desktop\Guangdong-city-2020.shp',encoding='gbk')
for record in sf.records():
    ax.add_geometries([record.geometry],ccrs.PlateCarree(),facecolor='none')
    #避免线条间首尾相连出现线乱飞的情况
ax.add_feature(cfeature.RIVERS.with_scale('50m'))# 添加河流，湖泊
ax.add_feature(cfeature.LAKES.with_scale('50m'))

注意，这里画图使用的是shapereader，上面读取文件使用的是shapefile ，个人认为这两者在使用途径上有一定区别。效果图如下所示：

 接着读取grib2文件并绘图，此处过程不一一赘述，完整代码如下：

import matplotlib.pyplot as plt
import numpy as np
import shapefile
import xarray as xr
from mpl_toolkits.basemap import Basemap
import cartopy.crs as ccrs
from cartopy.io import shapereader
import cartopy.feature as cfeature

#设置画图字体的大小
plt.rcParams.update({'font.size':20})
#解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']
#解决负号乱码问题
plt.rcParams['axes.unicode_minus'] = False
#设置画布和绘图区
fig = plt.figure(figsize=[12,18])
ax = fig.add_subplot(111,projection=ccrs.PlateCarree())
#-----------------------图像设置完毕-----------------------------

file = shapefile.Reader(r"C:\Users\20719\Desktop\Guangdong-city-2020.shp",encoding='GB2312')
# 通过创建reader类的对象进行shapefile文件的读取,encoding避免utf-8报错
records = file.records()

sf = shapereader.Reader(r'C:\Users\20719\Desktop\Guangdong-city-2020.shp',encoding='gbk')
for record in sf.records():
    ax.add_geometries([record.geometry],ccrs.PlateCarree(),facecolor='none')
ax.add_feature(cfeature.RIVERS.with_scale('50m'))# 添加河流，湖泊
ax.add_feature(cfeature.LAKES.with_scale('50m'))
#-----------读取地图边界--------------------

grib_file =r"C:\Users\20719\Desktop\data\data\BEBINCA_oneday\BEBINCA_oneday\Z_SURF_C_BABJ_20180807125215_P_CMPA_FRT_CHN_0P05_HOR-PRE-2018080712.GRB2"
ds = xr.open_dataset(grib_file, engine='cfgrib')
#-----------导入数据--------------------

# 提取降水数据（例如总降水量）
prec= ds['unknown']  #unknown表示总降水量，单位通常为mm

# 获取经度和纬度
lons = ds.longitude.data
lats = ds.latitude.data

# 广东省的大致经纬度范围：东经109°50′-117°20′，北纬20°13′-25°31′，实际应根据数据精确调整
lon_min, lon_max = 108, 118
lat_min, lat_max = 20,27

# 根据经纬度范围切片数据
subset = prec.sel(latitude=slice(lat_min, lat_max), longitude=slice(lon_min, lon_max))

# 获取调整后的经度和降水数据
lats = subset.latitude.data
lons = subset.longitude.data
data = subset.data
#-----------数据读取完毕--------------------

cbar_kwargs = {
    'orientation': 'horizontal',
    'label': 'Precipitation(mm)',
    'shrink': 0.02,
    'ticks': np.arange(0, 10 + 1, 1),
    'pad': -0.01,
    'shrink': 0.95
}
#-----------色标设置完毕--------------------

levels = np.arange(0, 18, 0.2)
cs = ax.contourf(lons,lats,data,cbar_kwargs=cbar_kwargs, cmap='Spectral_r')
cb = fig.colorbar(cs,orientation='horizontal',fraction = 0.1, pad = 0.08, shrink = 0.9)
#-----------数据绘制完毕--------------------

m = Basemap(llcrnrlon=109.0,#地图左边经度
    llcrnrlat=20.0,#地图下方纬度
    urcrnrlon=118.0,#地图右边经度
    urcrnrlat=26)#地图上方纬度
#------Basemap类设置完毕--------

#绘制纬线
parallels = np.arange(20,26,1)
m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3])
#绘制经线
meridians = np.arange(109,118,1)
m.drawmeridians(meridians,labels=[True,True,False,True],color='dimgrey',dashes=[1, 3])

#移除原来的坐标轴标签
plt.ylabel('')
plt.xlabel('')

#设置标题
plt.rcParams.update({'font.size':30})
ax.set_title(u'广东省xxxx年xx月xx日降水量图',fontsize= 25)

#出图
plt.show()

最终效果图如下所示，因为我是随便用了一个grib文件，所以没啥参考价值，主要重点在shp文件的绘制。

以上就是全部内容，小白的第一篇博客，恳请各位大佬批评指正🥰