MeteoInfoLab脚本示例：MODIS Sinusoidal投影HDF数据

MeteoInfo

孤蓝et 发表于 2019-1-29 01:17
王老师，打算将卫星云图的圆盘图（Geostationary Satellite View）转换为其他投影方式进行绘图，
看了您 ...

数据文件里读出来的数组变量是DimArray，有坐标和投影信息。由于此数据文件本身并不包含这些信息，是我们在脚本里定义的 x, y和proj，可以将其作为数组变量（比如 data）的属性，这样用 data.project投影函数就知道数组变量的原始投影了。具体参考这个例子：

1. fn = 'D:/Temp/hdf/FY4/FY4A-_AGRI--_N_REGC_1047E_L1-_FDI-_MULT_NOM_20190123043000_20190123043416_4000M_V0001.HDF'
   
2. f = addfile(fn)
   
3. data = f['NOMChannel01'][::-1,:]
   
4. data = data.astype('float')
   
5. data[data<0] = nan
   
6. x = linspace(-5496000.0,5496000.0, 2748)
   
7. y = linspace(-5496000.0,5496000.0, 2748)
   
8. lon0 = 104.74425
   
9. hh = 35781793
   
10. proj = projinfo(proj='geos', lon_0=lon0, h=hh)
    
11. xdim = dimension(x, 'x', 'X')
    
12. ydim = dimension(y, 'y', 'Y')
    
13. data.dims[0] = ydim
    
14. data.dims[1] = xdim
    
15. data.proj = proj
    
16. data = data['20:60','75:160']
    
17. 
    
18. #Project data
    
19. toproj = projinfo(proj='lcc', lon_0=104.137, lat_0=35.946, lat_1=30.0, lat_2=60.0)
    
20. data1 = data.project(toproj=toproj)
    
21. 
    
22. #Plot
    
23. ax = axesm(projinfo=toproj, gridline=True, griddx=10, griddy=10)
    
24. geoshow('cn_province', edgecolor='gray')
    
25. geoshow('country', edgecolor='k')
    
26. layer = imshowm(data1, proj=toproj)
    
27. colorbar(layer)

fn = 'D:/Temp/hdf/FY4/FY4A-_AGRI--_N_REGC_1047E_L1-_FDI-_MULT_NOM_20190123043000_20190123043416_4000M_V0001.HDF'<br /> f = addfile(fn)<br /> data = f['NOMChannel01'][::-1,:]<br /> data = data.astype('float')<br /> data[data<0] = nan<br /> x = linspace(-5496000.0,5496000.0, 2748)<br /> y = linspace(-5496000.0,5496000.0, 2748)<br /> lon0 = 104.74425<br /> hh = 35781793<br /> proj = projinfo(proj='geos', lon_0=lon0, h=hh)<br /> xdim = dimension(x, 'x', 'X')<br /> ydim = dimension(y, 'y', 'Y')<br /> data.dims[0] = ydim<br /> data.dims[1] = xdim<br /> data.proj = proj<br /> data = data['20:60','75:160']<br /> <br /> #Project data<br /> toproj = projinfo(proj='lcc', lon_0=104.137, lat_0=35.946, lat_1=30.0, lat_2=60.0)<br /> data1 = data.project(toproj=toproj)<br /> <br /> #Plot<br /> ax = axesm(projinfo=toproj, gridline=True, griddx=10, griddy=10)<br /> geoshow('cn_province', edgecolor='gray')<br /> geoshow('country', edgecolor='k')<br /> layer = imshowm(data1, proj=toproj)<br /> colorbar(layer)