MeteoInfoLab脚本示例：Brightness temperatures

MeteoInfo

National snow & ice center的Brightness temperatures卫星遥感产品有两种投影：Polar (https://nsidc.org/data/docs/daac/nsidc0001_ssmi_tbs.gd.html) 和EASE-Grid （https://nsidc.org/data/docs/daac/nsidc0032_ssmi_ease_tbs.gd.html），数据是2个字节的二进制文件，也就是短整形，读取后需要除以10得到真实温度数据。产品分北极和南极，这里示例两个投影南极地区的数据，北极数据根据数据说明进行简单的修改即可。

数据有两种分辨率：25km和12.5km，和频率有关。下面两个例子都是12.5km的分辨率。

Polar投影：

1. #Data source: https://nsidc.org/data/docs/daac/nsidc0001_ssmi_tbs.gd.html]https://nsidc.org/data/docs/daac/nsidc0001_ssmi_tbs.gd.html
   
2. #DMSP SSM/I-SSMIS Daily Polar Gridded Brightness Temperatures, Version 4
   
3. 
   
4. fn = 'D:/Temp/binary/tb_f17_20160413_v4_s91v.bin'
   
5. 
   
6. #Set projection
   
7. proj = projinfo(proj='stere', lat_ts=-70, lat_0=-90, lon_0=0)
   
8. xn = 632              #Column number
   
9. yn = 664              #Row number
   
10. dx = 12.5 * 1000      #X resolution (m)
    
11. dy = 12.5 * 1000      #Y resolution (m)
    
12. x0 = -3950 * 1000 + dx / 2     #Lower left x corner
    
13. y0 = -3950 * 1000 + dy / 2    #Lower left y corner
    
14. x = arange1(x0, xn, dx)
    
15. y = arange1(y0, yn, dy)
    
16. 
    
17. #Read data
    
18. data = binread(fn, [yn, xn], 'short')
    
19. data = data[::-1,:].astype('float')
    
20. data = data / 10
    
21. data[data==0] = nan   #0 is missing value
    
22. 
    
23. #Plot
    
24. axesm(projinfo=proj, gridline=True, griddx=30)
    
25. lworld = shaperead('D:/Temp/Map/country1.shp')
    
26. geoshow(lworld, edgecolor='k')
    
27. layer = imshowm(x, y, data, 20, cmap='BlAqGrYeOrRe', proj=proj)
    
28. colorbar(layer)
    
29. title('Brightness Temperatures')

#Data source: <a href=\"https://nsidc.org/data/docs/daac/nsidc0001_ssmi_tbs.gd.html\" target=\"_blank\">https://nsidc.org/data/docs/daac/nsidc0001_ssmi_tbs.gd.html</a>]https://nsidc.org/data/docs/daac/nsidc0001_ssmi_tbs.gd.html<br /> #DMSP SSM/I-SSMIS Daily Polar Gridded Brightness Temperatures, Version 4<br /> <br /> fn = 'D:/Temp/binary/tb_f17_20160413_v4_s91v.bin'<br /> <br /> #Set projection<br /> proj = projinfo(proj='stere', lat_ts=-70, lat_0=-90, lon_0=0)<br /> xn = 632              #Column number<br /> yn = 664              #Row number<br /> dx = 12.5 * 1000      #X resolution (m)<br /> dy = 12.5 * 1000      #Y resolution (m)<br /> x0 = -3950 * 1000 + dx / 2     #Lower left x corner<br /> y0 = -3950 * 1000 + dy / 2    #Lower left y corner<br /> x = arange1(x0, xn, dx)<br /> y = arange1(y0, yn, dy)<br /> <br /> #Read data<br /> data = binread(fn, [yn, xn], 'short')<br /> data = data[::-1,:].astype('float')<br /> data = data / 10<br /> data[data==0] = nan   #0 is missing value<br /> <br /> #Plot<br /> axesm(projinfo=proj, gridline=True, griddx=30)<br /> lworld = shaperead('D:/Temp/Map/country1.shp')<br /> geoshow(lworld, edgecolor='k')<br /> layer = imshowm(x, y, data, 20, cmap='BlAqGrYeOrRe', proj=proj)<br /> colorbar(layer)<br /> title('Brightness Temperatures')

EASE-Grid投影：

1. #Data source: https://nsidc.org/data/docs/daac/nsidc0032_ssmi_ease_tbs.gd.html]https://nsidc.org/data/docs/daac/nsidc0032_ssmi_ease_tbs.gd.html
   
2. #DMSP SSM/I-SSMIS Pathfinder Daily EASE-Grid Brightness Temperatures, Version 2
   
3. 
   
4. fn = 'D:/Temp/binary/ID2-F17-SH2015035D-V2.91V'
   
5. 
   
6. #Set projection
   
7. proj = projinfo(proj='laea', lat_0=-90, lon_0=0)
   
8. xn = 1441             #Column number
   
9. yn = 1441             #Row number
   
10. dx = 12.5 * 1000      #X resolution (m)
    
11. dy = 12.5 * 1000      #Y resolution (m)
    
12. x0 = -9000 * 1000 + dx / 2     #Lower left x corner
    
13. y0 = -9000 * 1000 + dy / 2    #Lower left y corner
    
14. x = arange1(x0, xn, dx)
    
15. y = arange1(y0, yn, dy)
    
16. 
    
17. #Read data
    
18. data = binread(fn, [yn, xn], 'short')
    
19. data = data[::-1,:].astype('float')
    
20. data = data / 10
    
21. data[data==0] = nan   #0 is missing value
    
22. 
    
23. #Plot
    
24. axesm(projinfo=proj, gridline=True, griddx=30, griddy=30)
    
25. lworld = shaperead('D:/Temp/Map/country1.shp')
    
26. geoshow(lworld, edgecolor='k')
    
27. layer = imshowm(x, y, data, 20, cmap='BlAqGrYeOrRe', proj=proj)
    
28. colorbar(layer)
    
29. title('Brightness Temperatures')

#Data source: <a href=\"https://nsidc.org/data/docs/daac/nsidc0032_ssmi_ease_tbs.gd.html\" target=\"_blank\">https://nsidc.org/data/docs/daac/nsidc0032_ssmi_ease_tbs.gd.html</a>]https://nsidc.org/data/docs/daac/nsidc0032_ssmi_ease_tbs.gd.html<br /> #DMSP SSM/I-SSMIS Pathfinder Daily EASE-Grid Brightness Temperatures, Version 2<br /> <br /> fn = 'D:/Temp/binary/ID2-F17-SH2015035D-V2.91V'<br /> <br /> #Set projection<br /> proj = projinfo(proj='laea', lat_0=-90, lon_0=0)<br /> xn = 1441             #Column number<br /> yn = 1441             #Row number<br /> dx = 12.5 * 1000      #X resolution (m)<br /> dy = 12.5 * 1000      #Y resolution (m)<br /> x0 = -9000 * 1000 + dx / 2     #Lower left x corner<br /> y0 = -9000 * 1000 + dy / 2    #Lower left y corner<br /> x = arange1(x0, xn, dx)<br /> y = arange1(y0, yn, dy)<br /> <br /> #Read data<br /> data = binread(fn, [yn, xn], 'short')<br /> data = data[::-1,:].astype('float')<br /> data = data / 10<br /> data[data==0] = nan   #0 is missing value<br /> <br /> #Plot<br /> axesm(projinfo=proj, gridline=True, griddx=30, griddy=30)<br /> lworld = shaperead('D:/Temp/Map/country1.shp')<br /> geoshow(lworld, edgecolor='k')<br /> layer = imshowm(x, y, data, 20, cmap='BlAqGrYeOrRe', proj=proj)<br /> colorbar(layer)<br /> title('Brightness Temperatures')

freekiller

赞