python绘制反向轨迹图和潜在源贡献值分析图

最新推荐文章于 2024-11-11 10:57:52 发布

萝卜一下

最新推荐文章于 2024-11-11 10:57:52 发布

阅读量1.4k

点赞数 6

文章标签： python 开发语言

本文链接：https://blog.csdn.net/weixin_41796280/article/details/131477641

版权

python绘制反向轨迹图和潜在源贡献值分析图

github地址

一、环境配置

1. 创建conda虚拟环境，安装cartopy等必要的地图绘制工具

创建conda虚拟环境

conda create -n map_draw python=3.9
conda activate map_draw

安装cartopy时使用conda install进行安装，一步搞定

conda install cartopy

下载其他依赖库

xarray # 读取气象数据 .nc文件
matplotlib
pykrige # 调用kriging插值算法
frykit # 大神写的依赖cartopy包画中国的包
openpyxl
cnmaps
scipy
basemap #需要离线下载whl进行安装 目前已经不再更新

这里安装的frykit程序包是大佬自己写的一个画中国地图的程序包，我们在此包的基础上进行画图，修改其中一些函数进行封装。
basemap需要下载whl轮子进行离线安装，首先去官网下载对应的basemap和pyproj，官网地址为https://www.lfd.uci.edu/~gohlke/pythonlibs/，该库目前已经不更新。因为在pysplit自己的例子中使用了basemap进行画图，本例中使用cartopy进行画图，尽量避开使用basemap。

安装pysplit库

pip install pysplit

github地址：pysplit

pyPSCF介绍

pyPSCF是获取潜在源贡献值的一个GUI工具，集成了Hysplit生成轨迹部分和PSCF潜在源贡献值分析部分，本例对该库进行修改，将PSCF部分引入，用于绘制潜在源贡献值分布情况。

github地址：pyPSCF

注意：

pysplit库对basemap有依赖，安装之前需要安装basemap，但是现在basemap已经不再更新，我们可以将mapdesigner.py模块中涉及到引用basemap的地方进行注释。pysplit依赖basemap就是用来绘制反向轨迹的，我们现在用cartopy进行绘制，所以可以放心注释掉相关内容，因为在本例中不会调用对应的函数。

需要注释的部分：

from mpl_toolkits.basemap import Basemap
MapDesign类中make_basemap方法

2. 安装hysplit或Meteoinfo软件，用于计算反向轨迹。

pysplit库和pyPSCF库都依赖于Hysplit离线软件，所以需要对Hysplit进行安装，或者安装Meteoinfo软件。本例中安装了Hysplit软件。

ARL官方支持在线绘制反向轨迹图。首先进入ARL官网，点Get/Run HYSPLIT，进入Run HYSPLIT Trajectory Model，可以看到有Compute archive trajectories选项，点进去根据配置就可以进行反向轨迹的绘制。

在线绘制反向轨迹图教程

linux环境下安装hysplit

首先进入ARL官网，点Get/Run HYSPLIT，在LINUX HYSPLIT下可以看到Download Public (unregistered) Version，进入网址下载预编译的linux hysplit，根据要求填写邮箱地址，不需要注册，点I agree，即可进入预编译版linux hysplit下载对应版本的源码。下载完成后，放到用户目录下，运行tar zxvf hysplit.v5.2.3_CentOS7.9.2009_public.tar.gz直接解压缩即可使用，如果需要图形界面，则需要根据官方教程，分别安装tcl和ImageMagick，其中安装tcl可以参考linux中安装tcl、Linux 下编译安装 tcl 和 tk，安装ImageMagick参考Linux下安装和使用imagemagick，本例中不使用图形界面。

下载地形信息数据ETOPO2v2c_f4.nc，放到conf/altitude_data目录下。

二、GDAS1数据

无论是离线使用Hysplit软件还是调用pysplit库进行反向轨迹的绘制，首先需要的数据是GDAS1数据。最简单的获取方法是进入官网直接下载。首先进入ARL官网，点Get/Run HYSPLIT，进入NOAA ARL Archived Data，可以看到GDAS one-degree archive (Dec 2004 - present)选项，其中FTP DATA中包括NOAA NOMADS Server (recent files only)，此部分的数据是最近几周的数据，其中包括current7days 数据，该数据是最近7天的数据，每隔6小时对数据进行更新，即每天3：33、9：33、15：33、21：33会对数据进行更新，后续会用到该数据。而https://www.ready.noaa.gov/data/archives/gdas1.则是2004年至今的全部数据，我们同样可以看到最近7天的数据。

gdas1数据命名格式：

gdas1表示数据类型，之后三个字母表示月份缩写，之后的两个数字表示年份，最后的w1~w5表示该月的第几周。

例如：gdas1.apr05.w1表示2005年四月份第一周的数据，gdas1.mar23.w3表示2023年三月份第三周的数据。

也可以通过FTP形式进行数据传输，HYSPLIT模型GDAS1数据的下载，

使用python调用ftplib进行数据下载，使用Python下载NOAA-GDAS1气象数据。

本例中对代码进行简单修改，直接运行python ./download/download_gdas1.py，可以下载最近7天的数据。同时修改了current7days的文件名称为当前时刻的gdas1数据命名格式。

如果需要下载指定日期的文件，可以修改data_downloader()函数的参数，例如 data_downloader(year=['23'],month=['jun'],week=['1'])

三、获取轨迹文件

将主要参数写入配置文件中，./conf/parameters/localParamBackTraj.json中，storage_dir表示轨迹文件生成目录，为该项目的相对路径；working_dir为Hysplit软件的工作目录，一般安装完成Hysplit软件后不修改；meteo_dir为gdas1数据的存放路径，本例中使用./download/download_gdas1.py进行数据下载，路径默认为./conf/gdas1_data；location为生成后向轨迹和PSCF的经纬度坐标；altitudes为海拔高度；runtime表示后向轨迹的时间；basename为生成轨迹文件的一个基本名称；years、month、hours分别为年、月、小时，是一个列表形式，可以输入多个值，当然也需要必要的gdas1文件；date表示取每个月的哪几天。其他参数是BackTrajHysplit.py模块所需参数，对比了两个生成轨迹文件的模块，发现bulk_trajgen.py模块生成轨迹文件速度更快，本例中暂不使用BackTrajHysplit.py模块生成轨迹文件。

配置好参数，bulk_trajgen.py模块直接调用get_trajectories()函数即可生成对应的轨迹文件，调用delate_file('trajectory')函数可对生成的轨迹文件进行删除，也可用于删除其他数据。

date 表示读取的污染源数据的时间，会根据这个时间寻找轨迹文件，如果不存在某时刻的轨迹文件，则不会画该时刻的PSCF

四、PSCF类说明

在运行run_PSCF之前，需要先配置好配置文件conf/parameters/localParamPSCF.json。

关键字	内容	格式	举例
cityname	城市名称	字符串	“石家庄市”
citycode	城市编号	字符串	“13010000”
location	经纬度坐标	列表数组	[38.03, 114.48]
altitudes	轨迹的海拔高度	int数字	1000
dirBackTraj	存放轨迹的路径	字符串	“trajectory”
prefix	生成轨迹名称关键字	字符串	“colgate”
backTraj	反向轨迹小时数	int	24或48或72
dateMin	最小北京时间	字符串	“2023-06-27 05:00:00”
dateMax	最大北京时间	字符串	“2023-06-30 05:00:00”
species	污染源类型	列表	[‘pm25’,‘pm10’]
threshold	污染源阈值	列表	[1]
wF	是否加权重	bool	一般默认true
wFmanual	是否自定义权重	bool	true
wFlim	权重设置	列表	[“0.30”,“0.6”,“0.85”]
wFval	权重设置	列表	[“0.08”,“0.35”,“0.725”,“1.0”]
smoothplot	平滑处理	bool	true
plotBT	是否画轨迹图	bool	true
plotPolar	是否画PSCF图	bool	true
cutWithRain	是否添加下雨影响	bool	false
mapMinMax	地图范围设置	字典	{“latmin”: 30.0,“latmax”: 56.0,“lonmin”: 90.0,“lonmax”: 125.0}

配置好配置文件，运行py_PSCF/run_PSCF.py中的pscf_run()函数即可生成对应的反向轨迹图和PSCF分布图。

五、时间对应

在生成轨迹文件之前必须下载对应的gdas1数据，get_trajectories()批量生成轨迹文件更快，conf/parameters/localParamBackTraj.json配置文件中对应的时间为UTC时间，所以后续调用PSCF时，对输入的污染源类型数据对应的时间进行了处理，conf/parameters/localParamPSCF.json配置文件中直接输入的为北京时间。