小白学习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()
测试数据和文件需要的请下面留言，感谢支持
239
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
