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

jl2587t

小白学习Basemap气象画地图的第五天（读取micaps站点数据，省级能见度分布）
这一帖子，主要介绍了三个重点：
1.micaps站点数据的读取
2.站点数据的插值
3.不均匀色标的生成
在下面的介绍中一一解释


读取micaps站点数据

这一步暂时不做解释了，因为光解释这个可以开一帖，代码是我网上找到，不是自己写的。其实大家不用纠结是否能看懂，能用就好了。主要是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）')


第五步画掩膜，美化图片

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

出图

function=‘cubic’

function=‘linear’

详细的代码解释大家可以下载附件，或者去CSDN论坛看：htps://blog.csdn.net/weixin_42372313/article/details/113784707

oydl1985

这图画的好像一个骷颅头！看来霾真的很“致命”

jl2587t

oydl1985 发表于 2021-2-10 19:32
这图画的好像一个骷颅头！看来霾真的很“致命”

湖南那个应该是雾，今天南方下雨

noisy

楼主，赞，学习中~~~

jl2587t

noisy 发表于 2021-2-11 15:26
楼主，赞，学习中~~~

{:eb348:}

花容天下重莲

楼主很厉害啊

jl2587t

花容天下重莲 发表于 2021-2-17 09:20
楼主很厉害啊

哈哈，我不厉害，我也是刚开始学，还在拾人牙慧

qinhantang10

请问这句是啥意思 buf = open('data_table.pickle', 'rb')，没找到路径。

jl2587t

qinhantang10 发表于 2021-02-17 16:17
请问这句是啥意思 buf = open('data_table.pickle', 'rb')，没找到路径。

哦哦，附件里少发了一个文件，留下邮箱，我发给你

来自安卓APP客户端

qinhantang10

78498062@qq.com,谢谢！请问这个20210210100000.000这是什么数据？micaps里gds的吗？