NCL读取MICAPS的4类数据

柿子柿子柿子

本想直接用《NCL代码分享-MICAPS第四类数据转NetCDF文件》里的方法的，毕竟“不要重复发明轮子”，结果发现不知道为什么下载下来的数据里面带有大量的空行，因而上面提到的帖子里的内容就无法直接使用了，不得已又“重复发明轮子”了。

1. function readMicaps3ModelScalar(filename:string)
   
2. begin
   
3.   a = asciiread(filename, -1, "string")
   
4.   ; 去头
   
5.   asciiwrite("tmp.tmp", a(1:))
   
6.   delete(a)
   
7. 
   
8.   a = asciiread("tmp.tmp", -1, "float")
   
9.   ; 删除临时文件
   
10.   system("rm -rf tmp.tmp")
    
11. 
    
12.   ; 时间
    
13.   tunits = "hours since 1900-01-01 00:00:00"
    
14.   t = cd_inv_calendar(a(0), a(1), a(2), a(3), 0, 0, tunits, 0)
    
15.   period = (/a(4)/)
    
16. 
    
17.   ; 层次
    
18.   lvl = (/a(5)/)
    
19. 
    
20.   ; 经纬度
    
21.   m = fspan(a(8), a(9), toint(a(12)))
    
22.   m@units = "degrees_east"
    
23.   n = fspan(a(10), a(11), toint(a(13)))
    
24.   n@units = "degrees_north"
    
25. 
    
26.   ; 数据
    
27.   dat = onedtond(a(19:), (/1, 1, 1, toint(a(13)), toint(a(12))/))
    
28.   dat!0 = "fcst"
    
29.   dat!1 = "period"
    
30.   dat!2 = "level"
    
31.   dat!3 = "lat"
    
32.   dat!4 = "lon"
    
33.   dat&fcst = t
    
34.   dat&period = period
    
35.   dat&level = lvl
    
36.   dat&lat = n
    
37.   dat&lon = m
    
38. 
    
39.   return dat
    
40. end
    
41. 
    
42. 
    
43. 
    
44. begin
    
45.   dir = systemfunc("ls dat/")
    
46. 
    
47.   do i = 0, dimsizes(dir) - 1
    
48.     a = readMicaps3ModelScalar("dat/" + dir(i))
    
49.     printVarSummary(a)
    
50.     delete(a)
    
51.   end do
    
52. end

复制代码

function的功能是传入数据文件的路径名，然后返回一个该4类数据对应的数组变量。

注1：考虑了万一以后会把多个起报时次、多个预报时效和多个预报层次拼接在一起的问题，返回的变量是个5维数组，依次是（起报时间，预报时效，层次，纬度，经度）。

注2：function中的“去头”是考虑到通常第一行会有一串数据说明的string，这串东西里面会有若干个数字，但是又暂时没找到什么方法跳过不读第一行，所以就暂时用了个临时文件tmp.tmp。

落然

谢谢，很有帮助！！！

woshiwzm

多的不说，肯定会用得到的