登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 jl2587t 于 2021-2-5 23:48 编辑
小白学习Basemap气象画地图的第三天（中国温度分布图，mask外部）首先还是感谢公众号（气象学家），代码和测试数据来自与他，不过这次有长进了，自己学会修改了。还是逐条向大家解释。
（和大家分享一个经验，在代码中查看某个函数或者变量的定义，可以利用快捷键快速寻找，pycharm里面是ctrl + 鼠标左键点击，这样的意义在于可以定位到最原始的构造函数，定义等位置，从而了解其用法，参数设置。要学会看__init__()里面的东西）
这次是用Basemap库画的，这个库在线安装好像已经停止了，只能下载离线包进行安装了，大家注意一下，首先效果图奉上： 看着不舒服的建议前往CSDN看，htps://blog.csdn.net/weixin_42372313/article/details/113705254
导入库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 -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())
作者准备了两个版本的源代码，一个注释简洁，一个详细。
| 
窥视卡
雷达卡
发表于 2021-2-5 23:35:07
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
