MeteoInfoLab脚本示例：FY4A CTT

MeteoInfo

由于FY4A数据中空间信息缺失，需要根据数据实际情况在脚本中补充。FY卫星圆盘数据是Geostationay satellite投影，在此投影的的数据网格是矩形，x, y坐标单位是米。这个数据的分辨率是4km (4000m)，数据网格中心坐标是(0, 0)，x方向格点数是2748，那么x方向的总长度是2748*4000，再除以2就可以得到x方向坐标起始和结束值为(-5496000.0,5496000.0)。

1. fn = 'D:/Temp/FY/FY4A-_AGRI--_N_DISK_1047E_L2-_CTT-_MULT_NOM_20190209140000_20190209141459_4000M_V0001.NC'
   
2. f = addfile(fn)
   
3. x = linspace(-5496000.0,5496000.0, 2748)
   
4. y = linspace(-5496000.0,5496000.0, 2748)
   
5. data = f['CTT'][::-1,:]
   
6. data[data>1000] = nan
   
7. data[data==-999] = nan
   
8. height = f['nominal_satellite_height'][:]
   
9. 
   
10. #Plot
    
11. lon0 = 104.7
    
12. ax = axesm(proj='geos', lon_0=lon0, h=height, gridlabelloc='all', griddx=30,
    
13.     griddy=30, gridline=True, frameon=False)
    
14. geoshow('coastline', color='k')
    
15. levs = arange(160, 311, 5)
    
16. layer = imshow(x, y, data, levs, proj=ax.proj)
    
17. colorbar(layer, shrink=0.8, xshift=15)
    
18. title('FY4A Cloud Top Temperature')
    
19. 
    
20. #savefig('D:/Temp/figure/fy4a_ctt.png', 500, 400, dpi=300)

复制代码

风格而才

麻烦问一下老师，这种投影我怎么只显示我研究的某个区域呢。

风格而才

风格而才 发表于 2024-4-29 00:19
麻烦问一下老师，这种投影我怎么只显示我研究的某个区域呢。

找到了哈，axis([-122, -75, 25, 50])

风格而才

麻烦 老师能帮忙看下风云4b数据画CTT嘛

MeteoInfo

风格而才 发表于 2025-12-17 15:54
麻烦 老师能帮忙看下风云4b数据画CTT嘛

1. fn = 'D:/Temp/satellite/FY/FY4B-_AGRI--_N_DISK_1050E_L2-_CTT-_MULT_NOM_20250303000000_20250303001459_4000M_V0001.NC'
   
2. f = addfile(fn)
   
3. 
   
4. data = f.CTT[::-1]
   
5. data[data==65535] = nan
   
6. data[data==-999] = nan
   
7. lon0 = f.nominal_satellite_subpoint_lon.values[0]
   
8. height = f.nominal_satellite_height.values[0]
   
9. ny, nx = data.shape
   
10. ex = (nx * 4000) / 2.
    
11. sx = -ex
    
12. ey = (ny * 4000) / 2.
    
13. sy = -ey
    
14. x = linspace(sx, ex, nx)
    
15. y = linspace(sy, ey, ny)
    
16. 
    
17. #Plot
    
18. proj = geolib.Geostationary(central_longitude=lon0, satellite_height=height)
    
19. ax = axesm(projection=proj, frameon=False, axison=False)
    
20. geoshow('coastline', color='k')
    
21. grid(True, tickvisible=True, tickposition='all')
    
22. levs = arange(160, 320, 2)
    
23. layer = imshow(x, y, data, levs, transform=proj)
    
24. colorbar(layer, shrink=0.8, xshift=18)
    
25. title('FY4B L2 CTT')

fn = 'D:/Temp/satellite/FY/FY4B-_AGRI--_N_DISK_1050E_L2-_CTT-_MULT_NOM_20250303000000_20250303001459_4000M_V0001.NC'<br /> f = addfile(fn)<br /> <br /> data = f.CTT[::-1]<br /> data[data==65535] = nan<br /> data[data==-999] = nan<br /> lon0 = f.nominal_satellite_subpoint_lon.values[0]<br /> height = f.nominal_satellite_height.values[0]<br /> ny, nx = data.shape<br /> ex = (nx * 4000) / 2.<br /> sx = -ex<br /> ey = (ny * 4000) / 2.<br /> sy = -ey<br /> x = linspace(sx, ex, nx)<br /> y = linspace(sy, ey, ny)<br /> <br /> #Plot<br /> proj = geolib.Geostationary(central_longitude=lon0, satellite_height=height)<br /> ax = axesm(projection=proj, frameon=False, axison=False)<br /> geoshow('coastline', color='k')<br /> grid(True, tickvisible=True, tickposition='all')<br /> levs = arange(160, 320, 2)<br /> layer = imshow(x, y, data, levs, transform=proj)<br /> colorbar(layer, shrink=0.8, xshift=18)<br /> title('FY4B L2 CTT')

风格而才

MeteoInfo 发表于 2025-12-18 09:42

收到  感谢