pyQt5结合Cartopy和matplotlib在界面中画micaps欧洲数值预报

一、概述

这则帖子主要介绍了如何在pyqt5中利用Cartopy画地图，画一些想要的信息（本文画的是micaps的欧洲数值预报）。大概介绍本文涉及的知识点，供大家参考：

1. pyqt5与matplotlib的结合。
2. matplotlib和Cartopy的结合。
3. micaps欧洲数值预报的读取。
4. 数据的绘制和平滑（插值）
5. 风杆过密的解决
6. 界面 matplotlib图像的交互，以及交互过程中数据图像和地图的刷新

先上一张最终效果图吸引一下大家：

二、界面绘制

这里利用pyqt5做了一个简单的界面，因为不是本次的重点，所以界面连半成品都算不上，只是大家可以完善，主要是介绍pyqt5和matplotlib的结合，下面是主函数main.py

from PyQt5.Qt import *
#自定义的一个类，用来做pyqt5和matplotlib的结合
# from My_Class import MyDataFigure
class Window(QWidget):
    def __init__(self):
        super().__init__()
        self.setWindowTitle("给个赞呗")
        self.showMaximized()#将主窗口最大化
        self.setup_ui()
    def setup_ui(self):
        #添加了三个控件  分别是左中右  中间控件到时候用来画图，左右两边摆各种按钮
        self.left_ql = QLabel(self)
        self.mid_ql = QLabel(self)#用来画图
        self.right_wt = QWidget(self)
        #为了方便区分，左边弄成黄色，右边弄成蓝色
        self.left_ql.setStyleSheet("background-color:yellow")
        self.right_wt.setStyleSheet("background-color:blue")
        #做了一个水平的动态布局，总共将窗口分成十份，中间占八份
        main_layout = QHBoxLayout()
        main_layout.addWidget(self.left_ql,1)
        main_layout.addWidget(self.mid_ql,8)
        main_layout.addWidget(self.right_wt,1)
        main_layout.setContentsMargins(0,0,0,0)#边框为0
        main_layout.setSpacing(0)#控件间隔为0
        self.setLayout(main_layout)
if __name__ == "__main__":
    import sys
    app = QApplication(sys.argv)
    window = Window()
    window.show()
    sys.exit(app.exec_())

上面代码会跑出这样一个窗口

界面就简单做成这个样子了，重点不在于介绍怎么做界面，左右两边大家可以自己加控件。

二、pyqt5和matplotlib的结合

使pyqt5和matplotlib的结合，实质上是利用了matplotlib中的一个类FigureCanvasQTAgg，可嵌入到pyqt5的QLabel控件中展示出来，也可像matplotlib的画布一样画图，下面提供My_Class.py的代码：

#以下引入的包一个都不能少
import matplotlib
matplotlib.use("Qt5Agg")
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg
import matplotlib.pyplot as plt
from matplotlib.figure import Figure

#需要用来嵌入界面的类，继承自FigureCanvasQTAgg
class MyDataFigure(FigureCanvasQTAgg):
    #构造函数
    #前面两个参数ql ,mid_ql是两个QLabel控件，不是必须的，是作者自行添加的，一个用来显示鼠标经纬度，一个用来使鼠标按下变换形状的，同理，大家需要传入什么参数，都可以写在__init__()中
    def __init__(self, ql ,mid_ql,width=100, height=100, dpi=30):
        self.figs = Figure(figsize=(width, height), dpi=dpi)
        super(MyDataFigure, self).__init__(self.figs)#前面的代码一条都不能少，类名要保持一致
        plt.rcParams['font.sans-serif'] = ['SimHei']#解决汉字乱码
        plt.rcParams['axes.unicode_minus'] = False#解决负号不显示
        print('创建成功')#下面就可以开始画图了

设计好画图类之后，需要在主界面实例化该类
首先引入类

from My_Class import MyDataFigure

然后实例化

#实例化自定义的类MyDataFigure，传入了两个控件参数，可以自己修改
        self.canvas_data = MyDataFigure(self.left_ql,self.mid_ql)
        #将类嵌入QLabel控件
        self.hboxlayout = QHBoxLayout(self.mid_ql)
        self.hboxlayout.addWidget(self.canvas_data)
        #设置边框为0
        self.hboxlayout.setContentsMargins(0,0,0,0)
        #调节画图的区域边界
        self.canvas_data.figs.subplots_adjust(left=0.05, right=1.3, top=0.9, bottom=0.5)

至此，主窗口main.py的代码量已全部完成，会在文末给出。

三、matplotlib和Cartopy的结合

这一部分相对简单一点，我们主要的目的是利用Cartopy在matplotlib上面画地图，之后所有的代码部分就在My_Class上面完成。
注意：运行仍是main.py

1.引入需要的包

相关包作用已经注释

#画图需要的包
import matplotlib.colors as colors#颜色
import cartopy.feature as cfeature#地图加载
import cartopy.crs as ccrs#投影方式
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER#经纬度转化
import matplotlib.ticker as mticker#x,y轴刻度显示

2.画地图

基础设置分为1.设置显示地图范围；2.创建子图；3.把区域加载到子图中。

self.extent = [70,140,20,60]#显示东经70-140，北纬20-60的区域
self.axes_map = self.figs.add_axes([0.03, 0, 0.94, 0.95], projection=ccrs.PlateCarree())#创建子图，投影普通投影
self.axes_map.set_extent(self.extent,crs=ccrs.PlateCarree())#设置范围，投影普通投影
#利用cartopy自带地图画海
self.axes_map.add_feature(cfeature.OCEAN.with_scale('110m'))
#利用cartopy自带地图陆地
self.axes_map.add_feature(cfeature.LAND.with_scale('110m'))
# #利用cartopy自带地图河流
self.axes_map.add_feature(cfeature.RIVERS.with_scale('110m'))
#利用cartopy自带地图湖泊
self.axes_map.add_feature(cfeature.LAKES.with_scale('110m'))

这一部分代码在后期会修改，为了刷新地图，会放到一个方法内，下文会介绍，这里是方便大家理解。
加载的地图分辨率为110m，是为了降低对计算机速度的要求

3.画国界

这个不多做介绍了，在以前的文章中有介绍，直接上代码

#读取CN-border-La.dat文件
        with open('CN-border-La.dat') as src:
            context = src.read()
            blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
            self.borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]
        # 画国界
        for line in self.borders:
            self.axes_map.plot(line[0::2], line[1::2], '-', color='gray',transform=ccrs.PlateCarree())

刷新地图时只需要执行画国界这部分就可以，不需要反复读取CN-border-La.dat。
注意：实际上作者是用了shp文件画的国界，为了交互，那样更快一些，但是毕竟是发表出来的，还是尽量不出错，就用了比较正式的，国界这种情况，一旦错了，就太敏感了

四、micaps欧洲数值预报的读取

1.read_mdfs.py

micaps欧洲数值预报是格点数据，这里只介绍一个，数据读取部分作者是在气象家园抄的别人的，在这里公布给大家read_mdfs.py

import struct
import datetime
import numpy as np

class MDFS_Grid:
    def __init__(self, filepath):
        f = open(filepath, 'rb')
        if f.read(4).decode() != 'mdfs':
            raise ValueError('Not valid mdfs data')
        self.datatype = struct.unpack('h', f.read(2))[0]
        self.model_name = f.read(20).decode('gbk').replace('\x00', '')
        self.element = f.read(50).decode('gbk').replace('\x00', '')
        self.data_dsc = f.read(30).decode('gbk').replace('\x00', '')
        self.level = struct.unpack('f', f.read(4))
        year, month, day, hour, tz = struct.unpack('5i', f.read(20))
        self.utc_time = datetime.datetime(year, month, day, hour) - datetime.timedelta(hours=tz)
        self.period = struct.unpack('i', f.read(4))
        start_lon, end_lon, lon_spacing, lon_number = struct.unpack('3fi', f.read(16))
        start_lat, end_lat, lat_spacing, lat_number = struct.unpack('3fi', f.read(16))
        lon_array = np.arange(start_lon, end_lon + lon_spacing, lon_spacing)
        lat_array = np.arange(start_lat, end_lat + lat_spacing, lat_spacing)
        isoline_start_value, isoline_end_value, isoline_space = struct.unpack('3f', f.read(12))
        f.seek(100, 1)
        block_num = lat_number * lon_number
        data = {
   }
        data['Lon'] = lon_array
        data['Lat'] = lat_array
        if self.datatype == 4:
            # Grid form
            grid = struct.unpack('{}f'.format(block_num), f.read(block_num * 4))
            grid_array = np.array(grid).reshape(lat_number, lon_number)
            data['Grid'] = grid_array
        elif self.datatype == 11:
            # Vector form
            norm = struct.unpack('{}f'.format(block_num), f.read(block_num * 4))
            angle = struct.unpack('{}f'.format(block_num), f.read(block_num * 4))
            norm_array = np.array(norm).reshape(lat_number, lon_number)
            angle_array = np.array(angle).reshape(lat_number, lon_number)
            # Convert stupid self-defined angle into correct direction angle
            corr_angle_array = 270 - angle_array
            corr_angle_array[corr_angle_array < 0] += 360
            data['Norm'] = norm_array
            data['Direction'] = corr_angle_array
        self.data = data

2.读取

看不懂没有关系，拿来直接用就可以了，接下来在My_Class.py中引入

from read_mdfs import MDFS_Grid#读取格点数据使用

简单应用

a = MDFS_Grid('ECMWF_HR_HGT_500_21010308.003')
self.lon = a.data['Lon']#经度
self.lat = a.data['Lat']#纬度
self.var = a.data['Grid']#数据

上述的代码是为了给大家介绍使用原理：及创建一个MDFS_Grid类，传入参数即文件路径（这里是21年1月3日08时起报的欧洲数值预报，第3小时预报，500hpa的高度场）

实际实现

第一步，在My_Class.py创建一个标识，用来判断所需数据是否成功读入

self.flg = {
   
            'HGT':False,
            'RH':False,
            'TMP':False,
            'UGRD':False,
            'VGRD':False
        }

当数据成功读入以后，相应部分改为Ture。
第二部分，要读取的不是一个文件，所以创建一个读取方法
这里没有读取经纬度，因为经纬度都一样，我们可以自己生产一个。减少程序冗余

#数据读取方法，第一个参数是路径，第二个是数据类型
    def read_data(self,filepath,data_type):
        if data_type == 'HGT':#高度
            a = MDFS_Grid(filepath)
            self.data_hgt = a.data['Grid']
            self.flg['HGT'] = True#读取成功
        if data_type == 'RH':#湿度
            a = MDFS_Grid(filepath)
            self.data_rh = a.data['Grid']
            self.flg['RH']