MeteoInfoLab脚本示例：Himawari standard data

MeteoInfo

Himawari卫星数据提供两种格式：netCDF和HSD数据（Himawari standard data），数据下载需要先在此网站注册：http://www.eorc.jaxa.jp/ptree/index.html，然后用提供的FTP下载数据。HSD数据是二进制数据，数据格式的详细描述见此文档：http://www.data.jma.go.jp/mscweb ... rs_guide_en_v12.pdf。这里示例读取HSD数据并绘图：

1. import struct
   
2. 
   
3. def read_h8(fn):
   
4.     #Read data header   
   
5.     f = open(fn, 'rb')
   
6.     hlen = 0
   
7.     #1 Basic information block
   
8.     f.read(282)
   
9.     hlen += 282
   
10.     #2 Data information block
    
11.     f.read(5)
    
12.     ncol, = struct.unpack('<h', f.read(2))
    
13.     nrow, = struct.unpack('<h', f.read(2))
    
14.     f.read(41)
    
15.     hlen += 50
    
16.     #3 Projection information block
    
17.     #f.read(127)
    
18.     f.read(19)
    
19.     sx, = struct.unpack('<f', f.read(4))
    
20.     sy, = struct.unpack('<f', f.read(4))
    
21.     f.read(127 - 27)
    
22.     hlen += 127
    
23.     #4 Navigation information block
    
24.     f.read(139)
    
25.     hlen += 139
    
26.     #5 Calibration information block
    
27.     f.read(147)
    
28.     hlen += 147
    
29.     #6 Inter-calibration information block
    
30.     f.read(259)
    
31.     hlen += 259
    
32.     #7 Segment information block
    
33.     #f.read(47)
    
34.     f.read(3)
    
35.     tns, = struct.unpack('b', f.read(1))
    
36.     ssn, = struct.unpack('b', f.read(1))
    
37.     fln, = struct.unpack('<h', f.read(2))
    
38.     f.read(40)
    
39.     hlen += 47
    
40.     #8 Navigation correction information block
    
41.     f.read(1)
    
42.     blen, = struct.unpack('<h', f.read(2))
    
43.     f.read(blen - 3)
    
44.     hlen += blen
    
45.     #9 Observation time information block
    
46.     f.read(1)
    
47.     blen, = struct.unpack('<h', f.read(2))
    
48.     f.read(blen - 3)
    
49.     hlen += blen
    
50.     #10 Error information block
    
51.     f.read(1)
    
52.     blen, = struct.unpack('<h', f.read(2))
    
53.     f.read(blen - 3)
    
54.     hlen += blen
    
55.     #11 Spare block
    
56.     f.read(259)
    
57.     hlen += 259
    
58.    
    
59.     f.close()
    
60. 
    
61.     #Read data
    
62.     data = binread(fn, [nrow, ncol], 'short', skip=hlen)
    
63.     data = data.astype('float')
    
64.     data[data<0] = nan
    
65.     return data, ncol, nrow, fln
    
66. 
    
67. #Read data files
    
68. segments = range(1, 11)
    
69. for segment in segments:
    
70.     fn = 'E:/Temp/himawari8/HS_H08_20170921_0410_B16_FLDK_R20_S%02i10.DAT' % segment
    
71.     print fn
    
72.     data1,ncol,nrow1,fln1 = read_h8(fn)
    
73.     if segment == segments[0]:
    
74.         data = data1
    
75.         fln = fln1
    
76.         nrow = nrow1
    
77.     else:
    
78.         data = concatenate([data, data1], axis=0)
    
79.         nrow += nrow1
    
80. data = data[::-1,:]
    
81. 
    
82. #Plot
    
83. sx = -5500000
    
84. sy = 5500000 - segments[-1] * 550 * 2000
    
85. x = arange1(sx, ncol, 2000)
    
86. y = arange1(sy, nrow, 2000)
    
87. ax = axesm(proj='geos', lon_0=140.7, h=35785863, gridlabel=True, gridline=True, frameon=False)
    
88. geoshow('country', edgecolor='b')
    
89. cmap = 'MPL_gist_gray'
    
90. levs = arange(800, 2001, 50)
    
91. layer = imshowm(x, y, data, levs, cmap=cmap, proj=ax.proj)
    
92. colorbar(layer, shrink=0.8)

import struct<br /> <br /> def read_h8(fn):<br />     #Read data header    <br />     f = open(fn, 'rb')<br />     hlen = 0<br />     #1 Basic information block<br />     f.read(282)<br />     hlen += 282<br />     #2 Data information block<br />     f.read(5)<br />     ncol, = struct.unpack('<h', f.read(2))<br />     nrow, = struct.unpack('<h', f.read(2))<br />     f.read(41)<br />     hlen += 50<br />     #3 Projection information block<br />     #f.read(127)<br />     f.read(19)<br />     sx, = struct.unpack('<f', f.read(4))<br />     sy, = struct.unpack('<f', f.read(4))<br />     f.read(127 - 27)<br />     hlen += 127<br />     #4 Navigation information block<br />     f.read(139)<br />     hlen += 139<br />     #5 Calibration information block<br />     f.read(147)<br />     hlen += 147<br />     #6 Inter-calibration information block<br />     f.read(259)<br />     hlen += 259<br />     #7 Segment information block<br />     #f.read(47)<br />     f.read(3)<br />     tns, = struct.unpack('b', f.read(1))<br />     ssn, = struct.unpack('b', f.read(1))<br />     fln, = struct.unpack('<h', f.read(2))<br />     f.read(40)<br />     hlen += 47<br />     #8 Navigation correction information block<br />     f.read(1)<br />     blen, = struct.unpack('<h', f.read(2))<br />     f.read(blen - 3)<br />     hlen += blen<br />     #9 Observation time information block<br />     f.read(1)<br />     blen, = struct.unpack('<h', f.read(2))<br />     f.read(blen - 3)<br />     hlen += blen<br />     #10 Error information block<br />     f.read(1)<br />     blen, = struct.unpack('<h', f.read(2))<br />     f.read(blen - 3)<br />     hlen += blen<br />     #11 Spare block<br />     f.read(259)<br />     hlen += 259<br />     <br />     f.close()<br /> <br />     #Read data<br />     data = binread(fn, [nrow, ncol], 'short', skip=hlen)<br />     data = data.astype('float')<br />     data[data<0] = nan<br />     return data, ncol, nrow, fln<br /> <br /> #Read data files<br /> segments = range(1, 11)<br /> for segment in segments:<br />     fn = 'E:/Temp/himawari8/HS_H08_20170921_0410_B16_FLDK_R20_S%02i10.DAT' % segment<br />     print fn<br />     data1,ncol,nrow1,fln1 = read_h8(fn)<br />     if segment == segments[0]:<br />         data = data1<br />         fln = fln1<br />         nrow = nrow1<br />     else:<br />         data = concatenate([data, data1], axis=0)<br />         nrow += nrow1<br /> data = data[::-1,:]<br /> <br /> #Plot<br /> sx = -5500000<br /> sy = 5500000 - segments[-1] * 550 * 2000<br /> x = arange1(sx, ncol, 2000)<br /> y = arange1(sy, nrow, 2000)<br /> ax = axesm(proj='geos', lon_0=140.7, h=35785863, gridlabel=True, gridline=True, frameon=False)<br /> geoshow('country', edgecolor='b')<br /> cmap = 'MPL_gist_gray'<br /> levs = arange(800, 2001, 50)<br /> layer = imshowm(x, y, data, levs, cmap=cmap, proj=ax.proj)<br /> colorbar(layer, shrink=0.8)

随缘

紧跟业务需求

wuwei2163

王老师真给力!!!

kmdqlb

王老师威武！

Flash

多谢大神分享

rjys1219

膜拜{:eb502:}

小渔007

王老师威武！

小渔007

王老师，请问下载了FileZilla软件但是为什么连不上ftp呢，是什么原因？