小白学习Basemap气象画地图的第三天（中国温度分布图，mask外部）

首先还是感谢公众号（气象学家），代码和测试数据来自与他，不过这次有长进了，自己学会修改了。还是逐条向大家解释。
（和大家分享一个经验，在代码中查看某个函数或者变量的定义，可以利用快捷键快速寻找，pycharm里面是ctrl + 鼠标左键点击，这样的意义在于可以定位到最原始的构造函数，定义等位置，从而了解其用法，参数设置。要学会看__init__()里面的东西）
这次是用Basemap库画的，这个库在线安装好像已经停止了，只能下载离线包进行安装了，大家注意一下，首先效果图奉上：

导入库

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import numpy as np
import shapefile
import xarray as xr
from mpl_toolkits.basemap import Basemap

绘图区设置

#设置画图字体的大小
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)

读取地图，设置mask路径（重点）

这里是重点部分，小白容易看不懂，注释比较多，大家不要耐烦，文末作者会给出注释少的代码版本

#读取shp格式的地图，有三个文件分别为.dbf,.shp,.shx缺一不可
#这里的shp文件其实就是一堆点，把各个国家圈出来
#sf得到的是一个类，存放了全球的地理信息
sf = shapefile.Reader("country1")
#sf.shapeRecords()里面放了各个国家地区的信息
#循环的目的是单个拿出来，找到中国  shape_rec是单个
for shape_rec in sf.shapeRecords():#shape_rec.record内部存放了单一国家的很多信息，比如名字，货币等等#其中shape_rec.record[2]放的是国家名#可以print(shape_rec.record)看看if shape_rec.record[2] == 'China':codes = []#用来存放移动路径(画图动作)#shape_rec.shape是一个类，图形类#里面三个属性shapeType：图形类型  points图形边界坐标  parts图形起始索引#解释一下parts属性，一个国家的边界可能不是全连在一起的，会分为一块一块，那么就相当于多个图形，存在shp文件内就不连续，parts里面就放了每个区域的起始索引（下标）pts = shape_rec.shape.points#上文已经说过了parts的意义，设想一下，两块区域的的起始坐标中间夹的不就是一块区域的所有坐标吗，但是最后一块区域没有结束值，所有加了[len(pts)]，这就是最后一个点的索引。prt = list(shape_rec.shape.parts) + [len(pts)]#下面的循环主要作用是：建立一个绘图路径，利用区块起始点的索引生成一个画图动作for i in range(len(prt) - 1):codes += [Path.MOVETO]#点移动codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)#画线codes += [Path.CLOSEPOLY]#这块画完，循环结束，下一块clip = Path(pts, codes)#利用数据和路径生成一个画图动作clip = PathPatch(clip, transform=ax.transData)#再加入ax的变换

大致原理是读取地图文件，大家俗称的shp文件，然后利用里面的点，得到一个或者多个画图路径（一般都是多个，国家边界很容易是很多段的），利用这个路径来把路径外部区域全部画成白的。细节可以看代码块的注释。

读取NC数据

这一部分比较简单

#这里是读取NC数据部分，上一个帖子已经介绍了，不再赘述
ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
lat = ds.latitude
lon = ds.longitude
data = (ds['t2m'][0,::-1,:] - 273.15) # 把温度转换为℃   [0,::-1,:]表示第一个时次、纬度反向

画温度分部图

这一部分和之前也大同小异

cbar_kwargs = {'orientation': 'horizontal','label': 'Temperature (℃)','shrink': 0.02,'ticks': np.arange(-25, 25 + 1, 5),'pad': -0.28,'shrink': 0.95
}
levels = np.arange(-25, 25 + 1, 1)
cs = data.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, cmap='Spectral_r')

到这里，咱们画出来的全球的温度分补（别忘了plt.show()）

利用掩膜路径画中国区域

其实就一条代码，简单的很，拿来用就好了

#添加掩膜路径，白化外部的分部
for contour in cs.collections:contour.set_clip_path(clip)

但是结果出现的问题，因为坐标变换的问题，生成的结果为

利用Basemap类调整图形

#生成一个Basemap类，用来变换坐标，画出合适尺寸的图
m = Basemap(llcrnrlon=72.0,#地图左边经度llcrnrlat=10.0,#地图下方纬度urcrnrlon=137.0,#地图右边经度urcrnrlat=55.0,#地图上方纬度resolution = None,#分辨率，这里设置为无projection = 'cyl')#投影方式：默认，圆柱投影
#读取图形文件，画中国行政边界
m.readshapefile('bou2_4l','China Map',color='k',linewidth=1.2)

readshapefile()函数的作用就是读取bou2_4l文件画图，结果为

图像修饰

在进行一些修饰，让图片更好看

#画纬度
parallels = np.arange(10,55,10)
m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3])
#画经度
meridians = np.arange(70,135,10)
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' 中国区域2018年1月平均气温',color='blue',fontsize= 25)#画南海
with open('CN-border-La.dat') as src:context = src.read()blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]
sub_ax = fig.add_axes([0.758, 0.249, 0.14, 0.155],projection=ccrs.LambertConformal(central_latitude=90, central_longitude=115))
for line in borders:sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',transform=ccrs.Geodetic())
sub_ax.set_extent([106, 127, 0, 25],crs=ccrs.PlateCarree())

完整代码，注释简洁版

最后有注释详细版，怕大家烦，就放最后了

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import numpy as np
import shapefile
import xarray as xr
from mpl_toolkits.basemap import Basemap#设置画图字体的大小
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)#读取shp格式的地图
sf = shapefile.Reader("country1")
#得到mask白化路径
for shape_rec in sf.shapeRecords():if shape_rec.record[2] == 'China':codes = []pts = shape_rec.shape.points#边界点prt = list(shape_rec.shape.parts) + [len(pts)]#区块起始索引for i in range(len(prt) - 1):codes += [Path.MOVETO]#点移动codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)#画线codes += [Path.CLOSEPOLY]#这块画完，循环结束，下一块clip = Path(pts, codes)#利用数据和路径生成一个画图动作clip = PathPatch(clip, transform=ax.transData)#再加入ax的变换
#########################至此，所需国界图形读取完毕####################
#这里是读取NC数据部分，上一个帖子已经介绍了，不再赘述
ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
lat = ds.latitude
lon = ds.longitude
data = (ds['t2m'][0,::-1,:] - 273.15) # 把温度转换为℃   [0,::-1,:]表示第一个时次、纬度反向
############################至此读取温度数据结束#######################下面画
cbar_kwargs = {'orientation': 'horizontal','label': 'Temperature (℃)','shrink': 0.02,'ticks': np.arange(-25, 25 + 1, 5),'pad': -0.28,'shrink': 0.95
}levels = np.arange(-25, 25 + 1, 1)
cs = data.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, cmap='Spectral_r')
#但是画出的全球温度分布
############################至此温度画完了######################
#添加掩膜路径，白化外部的分部
for contour in cs.collections:contour.set_clip_path(clip)#生成一个Basemap类，用来变换坐标，画出合适尺寸的图
m = Basemap(llcrnrlon=72.0,#地图左边经度llcrnrlat=10.0,#地图下方纬度urcrnrlon=137.0,#地图右边经度urcrnrlat=55.0,#地图上方纬度resolution = None,#分辨率，这里设置为无projection = 'cyl')#投影方式：默认，圆柱投影#读取图形文件，画中国行政边界
m.readshapefile('bou2_4l','China Map',color='k',linewidth=1.2)
######################至此图基本画完了，下面是一些修饰过程######################
#画纬度
parallels = np.arange(10,55,10)
m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3])
#画经度
meridians = np.arange(70,135,10)
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' 中国区域2018年1月平均气温',color='blue',fontsize= 25)
#画南海
with open('CN-border-La.dat') as src:context = src.read()blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]
sub_ax = fig.add_axes([0.758, 0.249, 0.14, 0.155],projection=ccrs.LambertConformal(central_latitude=90, central_longitude=115))
for line in borders:sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',transform=ccrs.Geodetic())
sub_ax.set_extent([106, 127, 0, 25],crs=ccrs.PlateCarree())
######################修饰完毕，出图######################
plt.savefig("China_mask_T2m.png", dpi=300, bbox_inches='tight')
plt.show()

完整代码，注释详细版

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import numpy as np
import shapefile
import xarray as xr
from mpl_toolkits.basemap import Basemap#设置画图字体的大小
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)#读取shp格式的地图，有三个文件分别为.dbf,.shp,.shx缺一不可
#这里的shp文件其实就是一堆点，把各个国家圈出来
#sf得到的是一个类，存放了全球的地理信息
sf = shapefile.Reader("country1")
#sf.shapeRecords()里面放了各个国家地区的信息
#循环的目的是单个拿出来，找到中国  shape_rec是单个
for shape_rec in sf.shapeRecords():#shape_rec.record内部存放了单一国家的很多信息，比如名字，货币等等#其中shape_rec.record[2]放的是国家名#可以print(shape_rec.record)看看if shape_rec.record[2] == 'China':codes = []#用来存放移动路径(画图动作)#shape_rec.shape是一个类，图形类#里面三个属性shapeType：图形类型  points图形边界坐标  parts图形起始索引#解释一下parts属性，一个国家的边界可能不是全连在一起的，会分为一块一块，那么就相当于多个图形，存在shp文件内就不连续，parts里面就放了每个区域的起始索引（下标）pts = shape_rec.shape.points#上文已经说过了parts的意义，设想一下，两块区域的的起始坐标中间夹的不就是一块区域的所有坐标吗，但是最后一块区域没有结束值，所有加了[len(pts)]，这就是最后一个点的索引。prt = list(shape_rec.shape.parts) + [len(pts)]#下面的循环主要作用是：建立一个绘图路径，利用区块起始点的索引生成一个画图动作for i in range(len(prt) - 1):codes += [Path.MOVETO]#点移动codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)#画线codes += [Path.CLOSEPOLY]#这块画完，循环结束，下一块clip = Path(pts, codes)#利用数据和路径生成一个画图动作clip = PathPatch(clip, transform=ax.transData)#再加入ax的变换#########################至此，所需国界图形读取完毕#####################这里是读取NC数据部分，上一个帖子已经介绍了，不再赘述
ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
lat = ds.latitude
lon = ds.longitude
data = (ds['t2m'][0,::-1,:] - 273.15) # 把温度转换为℃   [0,::-1,:]表示第一个时次、纬度反向
############################至此读取温度数据结束#######################下面画
cbar_kwargs = {'orientation': 'horizontal','label': 'Temperature (℃)','shrink': 0.02,'ticks': np.arange(-25, 25 + 1, 5),'pad': -0.28,'shrink': 0.95
}levels = np.arange(-25, 25 + 1, 1)
cs = data.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, cmap='Spectral_r')
#但是画出的全球温度分布
############################至此温度画完了######################
#添加掩膜路径，白化外部的分部
for contour in cs.collections:contour.set_clip_path(clip)#生成一个Basemap类，用来变换坐标，画出合适尺寸的图
m = Basemap(llcrnrlon=72.0,#地图左边经度llcrnrlat=10.0,#地图下方纬度urcrnrlon=137.0,#地图右边经度urcrnrlat=55.0,#地图上方纬度resolution = None,#分辨率，这里设置为无projection = 'cyl')#投影方式：默认，圆柱投影#读取图形文件，画中国行政边界
m.readshapefile('bou2_4l','China Map',color='k',linewidth=1.2)
######################至此图基本画完了，下面是一些修饰过程######################
#画纬度
parallels = np.arange(10,55,10)
m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3])
#画经度
meridians = np.arange(70,135,10)
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' 中国区域2018年1月平均气温',color='blue',fontsize= 25)#画南海
with open('CN-border-La.dat') as src:context = src.read()blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]
sub_ax = fig.add_axes([0.758, 0.249, 0.14, 0.155],projection=ccrs.LambertConformal(central_latitude=90, central_longitude=115))
for line in borders:sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',transform=ccrs.Geodetic())
sub_ax.set_extent([106, 127, 0, 25],crs=ccrs.PlateCarree())######################修饰完毕，出图######################
plt.savefig("China_mask_T2m.png", dpi=300, bbox_inches='tight')
plt.show()

测试数据和文件需要的请下面留言，感谢支持