python处理wrf气象数据_气象编程 | Python3之WRF的投影转换
有一种WRF输出的数据采用兰伯特双标准纬线投影,那么除非刚好需要同样的投影,想对这种数据进行处理的话往往要进行投影转换,WRF应该是有一套工具可以进行相关的处理,比如wrf-python,但是作为并不熟悉wrf、仅仅是使用WRF输出数据的小白,去使用WRF系工具的话学习成本就比较高了,如何用熟悉、更通用的工具实现这一投影转换呢?
难道不是设置几个投影参数,用常见的投影相关的包就可以实现了吗?
对,问题在于这个参数怎么设置?这个坑还是很坑的,好在最终找到一篇2018年的英文博客https://fabienmaussion.info/2018/01/06/wrf-projection/,加上我自己的尝试和理解,梳理出本文
采用WRF输出的数据格式为GrADS二进制码
关键
从我的角度来看,WRF的兰伯特双标准纬线投影有两个坑:椭球体参数
投影中心参数
完全不了解WRF输出的情况下
“啪”的一下,很快啊,观察到WRF输出的GrADS二进制数据对应的ctl文件中pdef如下所示1pdef 288 288 lcc 32.318 117.203 144.500 144.500 60.00000 30.00000 117.30000 3000.000 3000.000
然后查到pdef的lcc语法如下:
那么熟悉GDAL的小伙伴应该就会很开心地想:哟,这不就是SetLCC里需要的参数吗?咱这么写(默认WGS84地理坐标系):1from osgeo import osr 2 3(isize, jsize, latref, lonref, iref, jref, Struelat, Ntruelat, slon, dx, dy) = \ 4 (288, 288, 32.318, 117.203, 144.5, 144.5, 60, 30, 117.3, 3000, 3000) 5 6lcc = osr.SpatialReference() 7lcc.SetLCC(Struelat, Ntruelat, latref, lonref, iref * dx, jref * dy) 8lcc.SetLinearUnitsAndUpdateParameters("kilometre", 3000) 9geo = lcc.CloneGeogCS()10geo2lcc = osr.CoordinateTransformation(geo, lcc)11lcc2geo = osr.CoordinateTransformation(lcc, geo)
接下来通过geo2lcc和lcc2geo就可以愉快的在投影坐标和地理坐标之间反复横跳了不是?
当然我们需要验证。WRF输出数据中有XLAT和 XLONG两个要素,描述每个格点对应的经纬度(即地理坐标),加上很自然认为投影坐标是均匀地从西南角(0, 0)到东北角(isize-1, jsize-1),如果能够用我们的lcc2geo计算出的经纬度矩阵与XLAT和 XLONG一致那就说明可以放心地左右横跳1x = np.arange(isize)2y = np.arange(jsize)3x_mesh, y_mesh = np.meshgrid(x, y)4latlon = np.array(lcc2geo.TransformPoints(np.vstack([x_mesh.ravel(), y_mesh.ravel()]).T))5lat = latlon[:, 0].reshape(shape)6lon = latlon[:, 1].reshape(shape)
然而当我们计算出咱们的经纬度坐标lat和lon后,与XLAT和 XLONG作差(欧氏距离)快速画个图一看:1d = ((lat - lat_in_wrf)**2 + (lon - lon_in_wrf)**2)**0.523import matplotlib.pyplot as plt45fig = plt.figure()6ax = fig.add_subplot(1, 1, 1)7im = ax.imshow(d, cmap="Reds")8cb = fig.colorbar(im, ax=ax)9fig.set_tight_layout(True)
这时候就有弹幕说:“我不满意!”其实我也不满意,误差有点大,但是也没大到离谱,说明应该是投影参数有一点点问题。
修正椭球体参数
于是查资料发现,WRF的椭球体不是椭球体,是【芬芳的语气词】一个半径为6370000米的正球体!那么咱这么干:1lcc = osr.SpatialReference()2lcc.ImportFromProj4("+proj=longlat +a=6370000 +b=6370000 +no_defs")3lcc.SetLCC(Struelat, Ntruelat, latref, lonref, iref * dx, jref * dy)4lcc.SetLinearUnitsAndUpdateParameters("kilometre", 3000)5geo = lcc.CloneGeogCS()6geo2lcc = osr.CoordinateTransformation(geo, lcc)7lcc2geo = osr.CoordinateTransformation(lcc, geo)
再画个图瞧瞧
【芬芳的语气词】!但是可以发现最大误差变小了,说明有所改进,最小误差变大了,说明剩下的误差应该是水平偏移的成分多一些了。
继续查资料
修正投影中心地理坐标参数并计算投影坐标
终于,找到开头提到的博客。我虽然不会WRF,但是根据我的理解:WRF大概有两个区域,大区域里嵌套小区域,投影参数是按大区域的来,最终产出的数据是小区域的。这意味着,pdef中的latref和lonref并不是投影中心地理坐标,只是一个reference(洋屁,“参考”的意思),是用来计算真正的投影坐标用的,而真正的投影中心地理坐标,是MOAD_CEN_LAT和STAND_LON,在ctl文件中最下面像注释一样的东西里1@ global String comment MOAD_CEN_LAT = 32.402@ global String comment STAND_LON = 117.30
那么咱这么干:1MOAD_CEN_LAT = 32.40 2STAND_LON = 117.30 3 4lcc = osr.SpatialReference() 5lcc.ImportFromProj4("+proj=longlat +a=6370000 +b=6370000 +no_defs") 6lcc.SetLCC(Struelat, Ntruelat, MOAD_CEN_LAT, STAND_LON, 0, 0) 7lcc.SetLinearUnitsAndUpdateParameters("kilometre", 3000) 8geo = lcc.CloneGeogCS() 9geo2lcc = osr.CoordinateTransformation(geo, lcc)10lcc2geo = osr.CoordinateTransformation(lcc, geo)1112e, n, _ = geo2lcc.TransformPoint(lonref, latref)13false_west = -(iref-1) + e14false_south = -(jref-1) + n15x, y = np.arange(isize) + false_west, np.arange(jsize) + false_south16x_mesh, y_mesh = np.meshgrid(x, y)1718lonlat = np.array(lcc2geo.TransformPoints(np.vstack([x_mesh.ravel(), y_mesh.ravel()]).T))19lon = lonlat[:, 0].reshape(shape)20lat = lonlat[:, 1].reshape(shape)
再画个图瞧瞧
误差降了两个数量级啦,我已经比较满意了
总结一下WRF投影中心参数为MOAD_CEN_LAT和STAND_LON,采用的椭球体为半径为6370000米的正球体
pdef中的latref和lonref表示输出数据第jref行第iref列的地理坐标(经纬度),不代表投影中心(jref和iref是从1开始计数的)
讨论一下我觉得投影坐标的绝对数值由于直接受单位和偏移的影响可以任意变化,真正反映出空间信息的是投影坐标之间的差值(相对数值),所以我在最后的SetLCC投影参数中没有设置false_east和false_north偏移(最后两个参数)。虽然可以求出符合WRF输出数据的false_east和false_north,再调用一次lcc.SetLCC获得“完美投影关系”,但是似乎会让人感到更复杂。(其实作为非地理专业的我,这段话就已经有点绕了,权当作讨论再绕一点:最后一段代码中lcc.SetLinearUnits也是没有必要的,之后所有距离按米计算就行)
原文采用pyproj进行投影,由于我更熟悉GDAL所以换成了GDAL实现。但是吐槽一点,TransformPoint这个方法在公开地理坐标系下如WGS84传参或返回的地理坐标都是先纬度再经度(如本文第一次调用),在自定义地理坐标系下传参或返回的地理坐标却是先经度再维度(本文最后两次调用),这也有可能是我用得不好,有了解的读者麻烦指教一下
原文的WRF输出是nc,我用的WRF输出是GrADS二进制码,不影响这个投影的逻辑
原文误差在10^-5,而本文误差在10^-4,我觉得应该是nc文件提供的各参数小数位数比GrADS的ctl文件多
查资料时发现WRF也有基于WGS84的数据,不在本文讨论范围内
python处理wrf气象数据_气象编程 | Python3之WRF的投影转换相关推荐
- 机器学习与气象数据_气象大数据与机器学习联合实验室 大数据和气象的“联姻”...
气象大数据与机器学习联合实验室 大数据和气象的"联姻" 来源:<中国科学报> 时间:2017-02-13 13:36:28 作者:沈春蕾 我们每天都在看天气预报,大家会 ...
- python处理wrf气象数据_利用python-cdo处理气象数据
如果你不喜欢命令行的操作方式,那么你可以尝试使用python-cdo,利用python脚本语言的优势来处理气象数据.命令行的方式有其优势,比如简单易操作,可扩展性更强等,利用CDO的python接口也 ...
- 输出nc数据_气象数据处理的火箭加速器—CDO
happy科研 关注 Hi-新朋友,记得点蓝字关注我们哟 科研到头秃(7) 学渣 今天碰到了一个科研难题,想想都头大 学霸 学渣 模式运行得到了360个逐月输出的文件,NCL处理起来慢爆了,想哭 学霸 ...
- 合并多个nc数据_气象数据处理的火箭加速器—CDO
happy科研 关注 Hi-新朋友,记得点蓝字关注我们哟 科研到头秃(7) 学渣 今天碰到了一个科研难题,想想都头大 学霸 学渣 模式运行得到了360个逐月输出的文件,NCL处理起来慢爆了,想哭 学霸 ...
- cimiss数据_气象现代化成果汛期应用系列报道_中国气象网
中国气象报记者刘钊 编者按:2016年12月20日,由国家气象信息中心牵头建设的全国综合气象信息共享平台(CIMISS)正式业务化运行,标志着以CIMISS为核心的国省统一数据环境正式建立,标准.统一 ...
- 利用python从网页查找数据_利用Python模拟淘宝的搜索过程并对数据进行可视化分析...
数据挖掘入门与实战 公众号: datadw 本文讲述如何利用Python模拟淘宝的搜索过程并对搜索结果进行初步的数据可视化分析. 搜索过程的模拟:淘宝的搜索页面有两种形式, 一种形式是, 2019/2 ...
- mysql气象数据分析_气象行业 - 解决方案 - MySQL分布式数据库_开源数据库解决方案_数据处理技术提供商-爱可生...
1.数据源 新一代CIMISS系统所收集的结构化信息包括了28个类别:中国地面逐小时资料(新Z自动站),中国地面逐小时资料(一体化区域站),中国地面逐小时资料(无人站),中国地面分钟资料(新Z自动站) ...
- python 删除特定行数据_怎么用 Python 做数据分析实例
01 生成数据表 第一部分是生成数据表,常见的生成方法有两种,第一种是导入外部数据,第二种是直接写入数据. Excel 中的文件菜单中提供了获取外部数据的功能,支持数据库和文本文件和页面的多种数据源导 ...
- python调用通达信数据_[python]沪深龙虎榜数据导入通达信的自选板块并标注于k线图上...
[python] 沪深龙虎榜数据导入通达信的自选板块, 并标注于 K 线图上 1 #coding=utf-8 2 3 # 读取 '[paint]' 开头的 csv 文件 4 #copyright @ ...
- python玩转大数据_【小旭学长】大数据博士教你用python玩转时空大数据
好消息!好消息!手把手教你用python玩大数据 小旭学长的python大数据教程完结撒花,共26P录制完毕,总时长4小时.每10分钟的视频的录制加剪辑时间加起来都要两小时以上,讲得很细但是节奏也很快 ...
最新文章
- router6 QoS 1 基础知识
- ADO.NET Entity Data Model入门实例
- apache iotdb_Apache-IoTDB
- 利用 portupgrade快速更新通过ports安装的软件
- shiro如何保证session不失效_请问在不加锁的情况下如何保证线程安全?
- struts基本概念(2)
- 数据结构八-Trie树
- java utf8 byte_byte以及UTF-8的转码规则
- 苹果cms V10模板 手机端模板粉红色模板带会员中心
- 诗与远方:无题(二十七)- 写给我妹妹的一首诗
- 【Flink】Flink 基于事件序列最大值 AssignerWithPeriodicWatermarks
- 论文笔记_S2D.33_2015-ICCV_使用单个多尺度卷积网络,预测深度、表面法线和语义标签
- Python tkinter改变光标样式
- 解决windows虚机系统时间与北京时间相差8小时
- Java-面向对象构造函数 -(private private)关键字
- 解读Gartner2013应用交付市场魔力象限
- 黑白照片如何变彩色?亲测好用的方法分享
- 提高企业WiFi速度的快速简便的方法—Vecloud微云
- 关于评价指标的理解(TPR,FPR,TAR,FAR,FRR,ERR)
- php创建网址打不开,php网站无法打开怎么办
热门文章
- 使用 OneDrive 对电脑内的任意文件进行备份
- 有时间的时候没钱,有钱的时候没时间_天使Emily_一起游博客_一起游_17u.com
- linux c 端口复用,Linux C++ 网络编程学习系列(1)——端口复用实现
- springboot整合数据库
- 可分离变量的微分方程
- 戴尔服务器加无线网卡用不了,戴尔笔记本无线网卡驱动如何安装?(已解决)...
- 至强CPU型号系列的变化
- 小白DIY自己的系统镜像
- C语言求乘方、幂数、取余
- 为什么站点访问慢?请收好这份 Web 服务器性能提升的总结