MeteoInfoLab脚本示例：读取NOAA 16 AVHRR局地投影数据

MeteoInfo

局地投影文件是对极轨资料原始文件进行辐射订正、投影变换、几何纠正等预处理工作后生成的文件格式。相对于原始文件来说，局地文件具备已完成投影、几何纠正和文件存储空间小等优势，是气象部门常用的存储格式，常见的文件后缀包括LD2和LDF。局地文件由头文件和通道数据2部分组成，其中头文件占128字节，通道数据依照BSQ（band sequential）格式排列，按波段依次存储数据。具体的数据格式可以参考下面的文章：


谭明艳，李程，李辉，2011，基于ENVI/IDL的气象卫星局地投影数据的提取，广东气象，33（2）：62-64.


这里给出利用MeteoInfoLab脚本读取某个波段的数据并绘图的例子。首先利用Python读取二进制文件的模块（struct）将数据中的一些参数读取出来，然后利用MeteoInfoLab中的binread()函数读取某个波段的数据至一个二维数组并绘图。binread()函数的第一个参数是数据文件名，第二个参数是数据维的设置，第三个参数是数据类型，此外还有一些可选参数，skip参数指定了读取数据时跳过的字节数，byteorder指定了字节顺序（big_endian或者little_endian，缺省为little_endian）。局地投影数据每个数据是2个字节的整数（即short），数据读取后要除以10。


示例脚本：

1. import struct
   
2. fn = 'E:/Temp/jingjing/NOAA16_AVHRR_20050601.MOS.ldf'
   
3. 
   
4. #Read data header parameters
   
5. f = open(fn, 'rb')
   
6. fid, = struct.unpack('2s', f.read(2))
   
7. sid, = struct.unpack('<h', f.read(2))
   
8. f.read(4)
   
9. year, = struct.unpack('<h', f.read(2))
   
10. month, = struct.unpack('<h', f.read(2))
    
11. day, = struct.unpack('<h', f.read(2))
    
12. hour, = struct.unpack('<h', f.read(2))
    
13. minute, = struct.unpack('<h', f.read(2))
    
14. f.read(2)
    
15. channel_n, = struct.unpack('<h', f.read(2))
    
16. proj, = struct.unpack('<h', f.read(2))
    
17. xn, = struct.unpack('<h', f.read(2))
    
18. yn, = struct.unpack('<h', f.read(2))
    
19. xdelta, = struct.unpack('<f', f.read(4))
    
20. ydelta, = struct.unpack('<f', f.read(4))
    
21. lon_1, = struct.unpack('<f', f.read(4))
    
22. lon_2, = struct.unpack('<f', f.read(4))
    
23. raduis, = struct.unpack('<f', f.read(4))
    
24. lat_min, = struct.unpack('<f', f.read(4))
    
25. lat_max, = struct.unpack('<f', f.read(4))
    
26. lon_min, = struct.unpack('<f', f.read(4))
    
27. lon_max, = struct.unpack('<f', f.read(4))
    
28. f.close()
    
29. 
    
30. #Read one channel data
    
31. cn = 3    #4th channel
    
32. skipn = 128 + (xn * yn * 2) * cn
    
33. data = binread(fn, [yn, xn], 'short', skip=skipn)
    
34. data = data[::-1,:]
    
35. data = data.astype('float')
    
36. data[data<=0] = nan
    
37. data = data / 10
    
38. lon = arange1(lon_min, xn, xdelta)
    
39. lat = arange1(lat_min, yn, ydelta)
    
40. #Plot
    
41. axesm()
    
42. lworld = shaperead('D:/Temp/Map/country1.shp')
    
43. lchina = shaperead('D:/Temp/Map/bou2_4p.shp')
    
44. geoshow(lchina, edgecolor='b')
    
45. geoshow(lworld, edgecolor='b')
    
46. cmap = 'MPL_gist_gray'
    
47. levs = arange(2200, 3001, 50)
    
48. layer = imshowm(lon, lat, data, levs, cmap=cmap)
    
49. colorbar(layer)
    
50. t = datetime.datetime(year, month, day, hour, minute)
    
51. title('AVHRR (' + t.strftime('%Y-%m-%d') + ')')

import struct<br /> fn = 'E:/Temp/jingjing/NOAA16_AVHRR_20050601.MOS.ldf'<br /> <br /> #Read data header parameters<br /> f = open(fn, 'rb')<br /> fid, = struct.unpack('2s', f.read(2))<br /> sid, = struct.unpack('<h', f.read(2))<br /> f.read(4)<br /> year, = struct.unpack('<h', f.read(2))<br /> month, = struct.unpack('<h', f.read(2))<br /> day, = struct.unpack('<h', f.read(2))<br /> hour, = struct.unpack('<h', f.read(2))<br /> minute, = struct.unpack('<h', f.read(2))<br /> f.read(2)<br /> channel_n, = struct.unpack('<h', f.read(2))<br /> proj, = struct.unpack('<h', f.read(2))<br /> xn, = struct.unpack('<h', f.read(2))<br /> yn, = struct.unpack('<h', f.read(2))<br /> xdelta, = struct.unpack('<f', f.read(4))<br /> ydelta, = struct.unpack('<f', f.read(4))<br /> lon_1, = struct.unpack('<f', f.read(4))<br /> lon_2, = struct.unpack('<f', f.read(4))<br /> raduis, = struct.unpack('<f', f.read(4))<br /> lat_min, = struct.unpack('<f', f.read(4))<br /> lat_max, = struct.unpack('<f', f.read(4))<br /> lon_min, = struct.unpack('<f', f.read(4))<br /> lon_max, = struct.unpack('<f', f.read(4))<br /> f.close()<br /> <br /> #Read one channel data<br /> cn = 3    #4th channel<br /> skipn = 128 + (xn * yn * 2) * cn<br /> data = binread(fn, [yn, xn], 'short', skip=skipn)<br /> data = data[::-1,:]<br /> data = data.astype('float')<br /> data[data<=0] = nan<br /> data = data / 10<br /> lon = arange1(lon_min, xn, xdelta)<br /> lat = arange1(lat_min, yn, ydelta)<br /> #Plot<br /> axesm()<br /> lworld = shaperead('D:/Temp/Map/country1.shp')<br /> lchina = shaperead('D:/Temp/Map/bou2_4p.shp')<br /> geoshow(lchina, edgecolor='b')<br /> geoshow(lworld, edgecolor='b')<br /> cmap = 'MPL_gist_gray'<br /> levs = arange(2200, 3001, 50)<br /> layer = imshowm(lon, lat, data, levs, cmap=cmap)<br /> colorbar(layer)<br /> t = datetime.datetime(year, month, day, hour, minute)<br /> title('AVHRR (' + t.strftime('%Y-%m-%d') + ')')

schliezer

老师您好，我遇到了类似的问题，二进制文件是C#通过读取数据库产生的站点地面要素文件，写成了为grads准备的格式，看到上面的程序里基本是逐要素读取。现在关键是我不清楚每个要素所占的字节，以及要素之间的空格有多少，因为那个C#程序并不是我写的，请问还有更方便的读取方式吗？pickle.load试过也不行