IDL编程学习之风雷达与GPS数据同化

15195775117

读研时有上“数据同化”的课，天书一般，老师叫马Xun林，每次机关炮一样把课匆匆上完，对学生没什么耐心
我的理解，把不同源的数据一起用，大概就叫“数据同化”了。
这个帖子，我准备画个风场图，手头有风雷达数据，记录的是：时间-高度-水平风速-方位角-垂直风速；
还有GPS数据，记录的是：时间-经纬度
由于风雷达数据里没有经纬度，所以必须从GPS中取。两者的纽带是时间，但两者的时间间隔又不一样，怎么匹配是个问题。
当然，首先得分别对风雷达和GPS数据进行观察，而后读取，数据在附件中。

15195775117

对风雷达数据的说明和观察

打开风雷达数据，我们看到第一行是高度（单位m）
接着speed是水平风速，angle是方位角，vspeed是垂直风速（向上为正），这每3个量是一组，对应的时间也是一样的
通过对时间的统计，发现风雷达每10秒采集一次数据
数据中的异常值为999.9和99.99
比较讨厌的是speed后比angle和vspeed多了个冒号，这在读取时就复杂些了

15195775117

对GPS数据的观察

日期-时间-经度-纬度
这里也有异常值，

一开始我觉得是小数点差了2位，但即便小数点左移2位，跟其他数据也不连续，所以只能删除

15195775117

因为风雷达数据中speed后多了个冒号，跟后面时间hh:mm:ss中的冒号混淆，用程序读取会复杂些，眼下赶时间，我就不在这个问题上绕了，直接在表格中删除第一列
读取GPS数据，代码：
  a=read_ascii('C:\Users\YGDY\Desktop\time_to_time\gps.txt',DELIMITER=' ',data_start=0)
  x=a.field1
  h1=x[1,*]
  a=read_ascii('C:\Users\YGDY\Desktop\time_to_time\gps.txt',DELIMITER=':',data_start=0)
  x=a.field1
  m1=x[1,*]
  s1=x[2,*]
  a=read_ascii('C:\Users\YGDY\Desktop\time_to_time\gps.txt',DELIMITER=';',data_start=0)
  x=a.field1
  ;lat纬度lon经度
  lon=x[1,*]
  lat=x[2,*]
  print,[H1,M1,S1,LON,LAT]
发现个令人悲伤的现象，精度缩水！原本是小数点后6位的精度，四舍五入缩得还有点随意

估计需要设置输入格式（format关键字什么的），暂时不管这个问题了，勉强算是读取了

15195775117

读取风雷达数据，
代码：
a=read_ascii('C:\Users\YGDY\Desktop\time_to_time\wind.csv',DELIMITER=',',data_start=1)
x=a.field01
h2=x[0,*]
w=x[1:-1,*]
a=read_ascii('C:\Users\YGDY\Desktop\time_to_time\wind.csv',DELIMITER=':',data_start=1)
x=a.field1
m2=x[1,*]
s2=x[2,*]


W               FLOAT     = Array[60, 4899]
H2              FLOAT     = Array[1, 4899]
M2              FLOAT     = Array[1, 4899]
S2              FLOAT     = Array[1, 4899]

15195775117

现在，要把风和GPS的时间匹配上
风的时间间隔是10s，GPS的1~6秒都有，我想可以这样：
对于风的每一个时刻，找到GPS中与其时间差最近的来匹配经纬度，但是超过10s的就算缺失。
这个算法首先要计算时间差，我得使用julday的逆函数caldat。
caldat使用示例：
IDL> x1=julday(6,29,2018,0,0,0)
IDL> x2=julday(6,29,2018,0,0,1)
IDL> caldat,x2-x1,a,b,c,d,e,f
IDL> a,b,c,d,e,f
           1
           1
       -4713
          12
           0
      0.99999457597732544
用这个就可以求时间差了！

qq469015280

不明觉厉

赵家清欢

之后要学习！{:eb502:}