小白学习Basemap气象画地图的第五天（读取micaps站点数据，省级能见度分布）

这一帖子，主要介绍了三个重点：
1.micaps站点数据的读取
2.站点数据的插值
3.不均匀色标的生成
在下面的介绍中一一解释

读取micaps站点数据

这一步暂时不做解释了，因为光解释这个可以开一帖，代码是我网上找到，不是自己写的。其实大家不用纠结是否能看懂，能用就好了。

import struct
import datetime
import pickle
def create_dict(_dict, index):if index not in _dict.keys():_dict[index] = list()
corr_dtype = {1:'x', 2:'h', 3:'i', 4:'l', 5:'f', 6:'d', 7:'s'}
corr_size = {1:1, 2:2, 3:4, 4:4, 5:4, 6:8, 7:1}
buf = open('data_table.pickle', 'rb')
var_table = pickle.load(buf)
buf.close()
class MDFS_Station:def __init__(self, filepath):f = open(filepath, 'rb')if f.read(4).decode() != 'mdfs':raise ValueError('Not valid mdfs data')# Headerdtype = struct.unpack('h', f.read(2))[0]self.data_dsc = f.read(100).decode('gbk').replace('\x00', '')self.level = struct.unpack('f', f.read(4))[0]self.level_dsc = f.read(50).decode('gbk').replace('\x00', '')year, month, day, hour, min_, sec, tz = struct.unpack('7i', f.read(28))self.utc_time = datetime.datetime(year, month, day, hour, min_, sec) - datetime.timedelta(hours=tz)f.seek(100, 1)#288# Data block 1station_num = struct.unpack('i', f.read(4))[0] #292# Data block 2quantity_num = struct.unpack('h', f.read(2))[0] #294x = dict([(struct.unpack('h', f.read(2))[0], struct.unpack('h', f.read(2))[0]) for i in range(quantity_num)])# Data block 3data = {}for i in ['ID', 'Lon', 'Lat']:create_dict(data, i)for i in x.keys():create_dict(data, i)for _ in range(station_num):stid, stlon, stlat = struct.unpack('iff', f.read(12))data['ID'].append(stid)data['Lon'].append(stlon)data['Lat'].append(stlat)q_num = struct.unpack('h', f.read(2))[0]id_list = list()# iterate over q_numfor __ in range(q_num):var_id = struct.unpack('h', f.read(2))[0]if var_id % 2 == 0 and var_id not in range(22):var_info = 1else:var_info = var_table[var_id]id_list.append(var_id)var_value = struct.unpack(corr_dtype[var_info], f.read(corr_size[var_info]))if var_value and var_id % 2 != 0:var_value = var_value[0]data[var_id].append(var_value)# print(x.keys())for i in x.keys():if i not in id_list:# data[i].append(np.nan)data[i].append(None)self.data = data

主要是MDFS_Station（）类

正式开始

第一步导入类

from read_mdfs import MDFS_Station#另外写的一个类，用于读取micaps数据，当然不是我写的，抄的，忘记出处了，不好意思
import numpy as np
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.colors as mcolors
import shapefile
from mpl_toolkits.basemap import Basemap
import matplotlib as mpl

第二步读取micaps站点数据（重点）

x = MDFS_Station('20210210100000.000')
#读出的x是一个字典
#经纬度，每一个是一个list
lon = x.data['Lon']
lat = x.data['Lat']
#数据，1207代码代表了【水平能见度（人工）】，micaps数据各代码含义可以自行寻找
var = x.data[1207]

这里读取的是’20210210100000.000’是micaps站点数据。
说明：
站点数据和格点的区别：格点数据是有规律排列的比如1°*1°就是一个经度一个纬度一个数据，是网格点。
站点数据就是观测站提供的数据，观测站的分部无法很规律，类似于离散的点
画站点数据需要进行插值，插值成格点数据进行画图，下面会提到
关于micaps数据各代码含义可以留下邮箱或者网盘，作者给发micaps数据结构。

第三步插值成格点数据（重点）

插值经纬度

对经纬度进行插值，这里作者讲经纬度插值成了0.05°*0.05°的网格

olon = np.linspace(108,115,140)
olat = np.linspace(24,31,140)
olon,olat = np.meshgrid(olon,olat)

np.meshgrid()，将两个一维数组合成一个二维的网格坐标矩阵，就是把经纬度的两个list变成有关系的二维经纬网格

插值数据

生成一个插值规则

func = Rbf(lon,lat,var,function='cubic')

Rbf()是一种插值函数，好像是径向插值什么的，字面理解的意思应该就是以这一点开始，向外插值。
function=‘cubic’，就是插值的方式，列举了三种linear、cubic、gaussian。
gaussian作者的电脑比较烂，跑崩了，大家可以试一试。
另外两种作者下面有效果，自己看。
其他的就自己去摸索啦。

#根据插值规则生成新的数据，得到格点数据
var_new1 = func(olon,olat)

#筛选数据，把因为插值造成的小于0，大于30的数据变回0和30
var_new1[var_new1 < 0] = 0
var_new1[var_new1 > 30] = 30

第四步画图

图像设置，读取地图文件在这里不再介绍，可以看作者之前的帖子，有详细的说明
主要解释不均匀色板的生成
设置不均匀的色板的原因在于，如果使用均匀的色板，从能见度200m到30km均匀等分，图像不美观，能见度差的区域也不够凸显

设置层次

#设置层次，仔细看可以发现间隔是不均匀的（抄的micaps）
clevs = [0.00, 0.20,0.50,1.00,2.00,3.00,5.00,10.00,20.00,30.00]

建立对应的颜色列表

#建立对应的颜色列表（抄的micaps）
cdict = ['#6f2800','#9c00ff','#fe0100','#fc5506','#fbb249','#ffff00','#75ff30','#92e0f8','#c9ecff']
#利用颜色列表生成一个cmap类，用于画图
my_cmap = mcolors.ListedColormap(cdict)

归一化颜色

#BoundaryNorm是数据分组中数据归一化比较好的方法
#意思就是把层次全部归一化到0-1，然后映射到颜色上，比较抽象，如果不理解也没事，照葫芦画瓢，照做就好了
norm = mcolors.BoundaryNorm(clevs,my_cmap.N)

画图画色板

#画图
cs = ax.contourf(olon,olat,var_new1,levels=clevs, cmap=my_cmap,norm =norm )
#画色板
plt.rcParams.update({'font.size':10})
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=my_cmap), ax = ax, label='能见度（km）')

第五步画掩膜，美化图片

这一步也不做详细解释了，前期的帖子也有，大家可以去看

#mask掩膜
for contour in cs.collections:contour.set_clip_path(clip)
#创建Basemap类
m = Basemap(llcrnrlon=108.0,#地图左边经度llcrnrlat=24.0,#地图下方纬度urcrnrlon=115.0,#地图右边经度urcrnrlat=31.0,#地图上方纬度resolution = None,#分辨率，这里设置为无projection = 'cyl')#投影方式：默认，圆柱投影
#画省界和地级市界
m.readshapefile('Hunan_Province','Map',color='k',linewidth=1.2)
m.readshapefile('Hunan_city','city Map',color='k',linewidth=1.2)#画纬度
plt.rcParams.update({'font.size':10})
parallels = np.arange(23,32,1)
m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3])
#画经度
meridians = np.arange(108.5,115,1)
m.drawmeridians(meridians,labels=[True,True,False,True],color='dimgrey',dashes=[1, 3])
#移除原来的坐标轴标签
plt.ylabel('')
plt.xlabel('')
#设置标题
plt.rcParams.update({'font.size':20})
ax.set_title(u'湖南省2021年2月10日10时能见度',color='blue',fontsize= 25)bill2 =  113.0
tip2  = 28.21
plt.scatter(bill2, tip2,marker='*',s=150 ,color ="orange" )#画个点plt.rcParams.update({'font.size':20})
plt.text(bill2+0.2, tip2, u"长沙市", fontsize= 20 )#在斜上方写个字plt.show()

出图

function=‘cubic’

function=‘linear’

完整代码

from read_mdfs import MDFS_Station#另外写的一个类，用于读取micaps数据，当然不是我写的，抄的，忘记出处了，不好意思
import numpy as np
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.colors as mcolors
import shapefile
from mpl_toolkits.basemap import Basemap
import matplotlib as mpl#-------------读取micaps数据，读取的站点数据-----------
#站点数据和格点的区别：格点数据是有规律排列的比如1°*1°就是一个经度一个纬度一个数据，是网格点
#站点数据就是观测站提供的数据，观测站的分部无法很规律，类似于离散的点
#画站点数据需要进行插值，插值成格点数据进行画图，下面会提到
x = MDFS_Station('20210210100000.000')
#读出的x是一个字典
#经纬度，每一个是一个list
lon = x.data['Lon']
lat = x.data['Lat']
#数据，1207代码代表了【水平能见度（人工）】，micaps数据各代码含义可以自行寻找
var = x.data[1207]#-------------插值成格点数据-（这里是重点）----------
#对经纬度进行插值，这里作者讲经纬度插值成了0.05°*0.05°的网格
olon = np.linspace(108,115,140)
olat = np.linspace(24,31,140)
#np.meshgrid()，将两个一维数组合成一个二维的网格坐标矩阵，就是把经纬度的两个list变成有关系的二维经纬网格
olon,olat = np.meshgrid(olon,olat)#Rbf()是一种插值函数，好像是径向插值什么的，字面理解的意思应该就是以这一点开始，向外插值
#function='cubic'，就是插值的方式，列举了三种linear、cubic、gaussian
#gaussian作者的电脑比较烂，跑崩了，大家可以试一试
#另外两种作者下面有效果，自己看
#其他的就自己去摸索啦
func = Rbf(lon,lat,var,function='cubic')
#对var进行插值，得到格点var_new1
var_new1 = func(olon,olat)
#筛选数据，把因为插值造成的小于0，大于30的数据变回0和30
var_new1[var_new1 < 0] = 0
var_new1[var_new1 > 30] = 30#-------------开始画图-----------
#解决中文乱码问题
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('Hunan_Province')#读取shp地图文件
shapes = sf.shapes()#获取shapefile的类
codes = []#用来存放移动路径(画图动作)
pts = shapes[0].points#边界点，这里只有一个类所以用了shapes[0]，没做循环
prt = list(shapes[0].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的变换#这里是重点，设置不均匀的色板（原因在于，如果使用均匀的色板，从能见度200m到30km均匀等分，图像不美观，能见度差的区域也不够凸显）
#设置层次，仔细看可以发现间隔是不均匀的（抄的micaps）
clevs = [0.00, 0.20,0.50,1.00,2.00,3.00,5.00,10.00,20.00,30.00]
#建立对应的颜色列表（抄的micaps）
cdict = ['#6f2800','#9c00ff','#fe0100','#fc5506','#fbb249','#ffff00','#75ff30','#92e0f8','#c9ecff']
#利用颜色列表生成一个cmap类，用于画图
my_cmap = mcolors.ListedColormap(cdict)
#BoundaryNorm是数据分组中数据归一化比较好的方法
#意思就是把层次全部归一化到0-1，然后映射到颜色上，比较抽象，如果不理解也没事，照葫芦画瓢，照做就好了
norm = mcolors.BoundaryNorm(clevs,my_cmap.N)
#画图
cs = ax.contourf(olon,olat,var_new1,levels=clevs, cmap=my_cmap,norm =norm )
#画色板
plt.rcParams.update({'font.size':10})
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=my_cmap), ax = ax, label='能见度（km）')#mask掩膜
for contour in cs.collections:contour.set_clip_path(clip)
#创建Basemap类
m = Basemap(llcrnrlon=108.0,#地图左边经度llcrnrlat=24.0,#地图下方纬度urcrnrlon=115.0,#地图右边经度urcrnrlat=31.0,#地图上方纬度resolution = None,#分辨率，这里设置为无projection = 'cyl')#投影方式：默认，圆柱投影
#画省界和地级市界
m.readshapefile('Hunan_Province','Map',color='k',linewidth=1.2)
m.readshapefile('Hunan_city','city Map',color='k',linewidth=1.2)#画纬度
plt.rcParams.update({'font.size':10})
parallels = np.arange(23,32,1)
m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3])
#画经度
meridians = np.arange(108.5,115,1)
m.drawmeridians(meridians,labels=[True,True,False,True],color='dimgrey',dashes=[1, 3])
#移除原来的坐标轴标签
plt.ylabel('')
plt.xlabel('')
#设置标题
plt.rcParams.update({'font.size':20})
ax.set_title(u'湖南省2021年2月10日10时能见度',color='blue',fontsize= 25)bill2 =  113.0
tip2  = 28.21
plt.scatter(bill2, tip2,marker='*',s=150 ,color ="orange" )#画个点plt.rcParams.update({'font.size':20})
plt.text(bill2+0.2, tip2, u"长沙市", fontsize= 20 )#在斜上方写个字plt.show()