爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 7519|回复: 8

用Meteoinfo做邻近插值遇到了一些问题,求老师指教QAQ

[复制链接]

新浪微博达人勋

发表于 2021-7-27 11:21:27 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
参考王老师给的脚本,用邻近插值方法做格点到站点的插值,但出现了一些问题。第一次用meteoinfo,不太熟练,希望大家能给予帮助,非常感谢!脚本如下:
stdata = readtable( 'C:/Users/xsx/Desktop/hly/cra40/lat_lon.csv', delimiter=',', format='%s%f%f')
x = stdata['Latitude']
y = stdata['Longitude']
print stdata
f = addfile('C:/Users/xsx/Desktop/hly/cra40/2020ERA5day-2.nc')
psv = f['t2m']
print psv
#Get time dimension length
tn = psv.dimlen(0)
#Loop
for i in range(0, tn):
    #Get dimension array
    ps = psv[i,'19:28','100:120']
    #Interpolate to stations
    d=interp2d(ps, x, y, kind='neareast')
    #Add column to table data
    colname = 't2m' + str(i)
    stdata.addcoldata(colname, '%d', ps)
#Save table data to a file
fn = 'C:/Users/xsx/Desktop/hly/cra40/test_st.csv'
stdata.savefile(fn)
print 'Finish...'


问题1. 读取csv数据,print显示读取的第一列带有引号,如下:
Latitude        Longitude
'108.45'        22.68
'109.18'        21.67
'106.95'        23.75
...(后面的省略)

问题2. 我的要素是double格式的,是否需要转换为float?找到了一个ones(*,dtype = '*')的函数,但不知该如何使用orz,报错信息如下:



气象家园图片.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2021-7-27 12:37:09 | 显示全部楼层
stdata = readtable( 'C:/Users/xsx/Desktop/hly/cra40/lat_lon.csv', delimiter=',', format='%s%f%f')
改为:
stdata = readtable( 'C:/Users/xsx/Desktop/hly/cra40/lat_lon.csv', delimiter=',', format='%f%f%f')
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-7-27 16:36:10 | 显示全部楼层
MeteoInfo 发表于 2021-7-27 12:37
stdata = readtable( 'C:/Users/xsx/Desktop/hly/cra40/lat_lon.csv', delimiter=',', format='%s%f%f')
...

老师,脚本跑通了,但是结果不太对呢
我用双线性插值的结果和观测值误差不大,但邻近插值得到的值偏小了10℃以上,是脚本还存在问题么?能麻烦王老师再帮忙看看吗?非常感谢!
观测值如下:
纬度(度)        经度(度)        年(年)        月(月)        日(日)        平均气温(摄氏度(℃))
22.68        108.45        2020                  1                  1                19.4
21.67        109.18        2020                  1                  1                20.3
23.75        106.95        2020                  1                  1                19.5

邻近插值后结果如下:
Longitude        Latitude        t2m0
22.68        108.45        269.6554688   -3.495
21.67        109.18        269.655573     -3.495
23.75        106.95        269.6487454   -3.502

密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-7-28 08:04:22 | 显示全部楼层
hlyahs 发表于 2021-7-27 16:36
老师,脚本跑通了,但是结果不太对呢
我用双线性插值的结果和观测值误差不大,但邻近插值得到 ...

插值算法应该没有问题,你看看是不是数据读错了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-7-28 09:31:37 | 显示全部楼层
MeteoInfo 发表于 2021-7-28 08:04
插值算法应该没有问题,你看看是不是数据读错了

老师您好,我核对了数据读取,print出来的站点信息是正确的。
不过我看您的脚本里,对站点经纬度的读取,x为lon,y为lat,即(lon,lat),而我的nc数据维度排列是(time,lat,lon),请问这会造成影响吗?我是否需要用命令将nc数据的纬度顺序调换为(time,lon,lat)后,再进行插值呢?
以及,老师,我的lat为[89.875..-89.875],是否需要调换为[-89.875..89.875]呢?
非常感谢您!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-7-28 12:42:41 | 显示全部楼层
hlyahs 发表于 2021-7-28 09:31
老师您好,我核对了数据读取,print出来的站点信息是正确的。
不过我看您的脚本里,对站点经纬度的读取 ...

你看看x, y是不是取反了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-7-28 16:47:02 | 显示全部楼层
MeteoInfo 发表于 2021-7-28 12:42
你看看x, y是不是取反了

王老师您好!
我确认了x,y,并没有读取错误。

d = interp2d(ps,x,y, kind='neareast'  
d是插值结果,但我print d后,得出结果如下:
array([NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN])(总计91个站点的插值,全部为NaN)

这看起来似乎没有插值成功呢,但为何插值失败,却能够生成csv文件并且有数值呢?
这是脚本运行结束后,输出的test_st.csv里的部分内容:
Latitude        Longitude        t2m0        t2m1        t2m2        t2m3
22.68        108.45        269.6554688        273.0030005        265.370255        268.5839439
21.67        109.18        269.655573        273.0394842        265.4666761        268.6766645
23.75        106.95        269.6487454        273.0703389        265.5608039        268.7718868
23.98        109.7        269.6386863        273.0970763        265.6501368        268.8651807
22.95        111        269.6195585        273.1072917        265.723521        268.9496143 (虽然有数值,但并不是插值结果)

老师,我的站点经纬度读入是(lon,lat),NC数据是(time,lat,lon),会不会是这个原因造成的呢?
一直叨扰您实在不好意思,我会尽量把问题描述清楚,还请王老师继续指教
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-8-2 10:01:03 | 显示全部楼层
MeteoInfo 发表于 2021-7-28 12:42
你看看x, y是不是取反了

老师您好,我这两天又捣鼓了一下脚本,读取时修改为psv = f['t2m'][:,::-1,:],发现插值结果不再是NaN了(虽然值不正确),所以我还是想尝试将nc资料的维度顺序调整一下,想请问老师,用meteoinfo的话,该如何实现将nc资料(time,lat,lon),变为(time,lon,lat)呢?
非常感谢您!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-8-2 11:23:34 | 显示全部楼层
hlyahs 发表于 2021-8-2 10:01
老师您好,我这两天又捣鼓了一下脚本,读取时修改为psv = f['t2m'][:,::-1,:],发现插值结果不再是NaN了( ...

time, lat, lon是正确的顺序,如果要调换维顺序可以用数组的swapaxes方法。或者你把数据文件和脚本程序发上来,我有空了看看。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表