请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 5189|回复: 11

[经验总结] 利用satpy处理风云四号B星的GHI仪器数据

[复制链接]

新浪微博达人勋

发表于 2022-11-18 17:59:17 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 iui 于 2022-11-18 17:58 编辑

satpy是一个处理气象遥感数据的一个python库,用于读取和处理气象遥感数据并将其写入各种图像和数据文件格式,具体介绍详见官方网站,https://satpy.readthedocs.io/en/latest/overview.html
      satpy在2022.11.11号更新到最新版本(0.38.0),加入了读取GHI数据的功能,详见github,https://github.com/pytroll/satpy/blob/main/CHANGELOG.md

风云四号B星上搭载的快速成像仪(GHI),数据时间分辨率最高达一分钟,空间分辨率最高达250m,可以下载从2022年6月1号-至今的数据(包括一级数据和一些大气产品),具体介绍详见官方网站,http://satellite.nsmc.org.cn/PortalSite/Data/Satellite.aspx。  

下面分享一个使用satpy处理FY4B卫星GHI数据的小程序,结果如图所示。  如有错误,请指正,谢谢。

c01-20220612080300.png
  1. import cartopy.crs as ccrs
  2. import matplotlib.pyplot as plt
  3. import scipy.io as sio
  4. import datetime
  5. from datetime import timedelta
  6. import matplotlib as mpl
  7. import os, glob
  8. from satpy.scene import Scene

  9. mpl.rc('font', family='Times new roman', size=16, weight='light')

  10. lonmin, latmin, lonmax, latmax = 115, 38, 118, 42
  11. extent = [lonmin, lonmax, latmin, latmax]
  12. chp = sio.loadmat('D:/data/prvbound.mat')#地图文件
  13. ir_channel = 'C01'
  14. path = r'D:\data\250M\202206'
  15. for index, filename in enumerate(os.listdir(path)):
  16.     l1 = os.path.join(path, filename)
  17.     filenames = glob.glob(l1)
  18.     scn = Scene(reader="ghi_l1", filenames=filenames)
  19.     scn.load([ir_channel])
  20.     crop_scn = scn.crop(ll_bbox=(lonmin, latmin, lonmax, latmax))
  21.     sat_w_tb_vis = crop_scn[ir_channel].values
  22.     area_def = crop_scn[ir_channel].attrs['area']

  23.     crs = area_def.to_cartopy_crs()
  24.     proj = ccrs.PlateCarree()

  25.     fig, ax = plt.subplots(figsize=[12, 9.6], subplot_kw={'projection': proj})
  26.     ax.plot(chp['long'], chp['lat'], 'yellow', lw=0.3)
  27.     ax.set_extent(extent, crs=proj)

  28.     img = plt.imshow(sat_w_tb_vis, origin='upper', cmap=plt.cm.binary_r, extent=crs.bounds, transform=crs)
  29.     gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, alpha=0)
  30.     gl.xlabels_top = False  # 关闭顶端的经纬度标签
  31.     gl.ylabels_right = False  # 关闭右侧的经纬度标签
  32.     # 转换世界时-北京时
  33.     utc_time = filename[44:58]
  34.     t_utc = datetime.datetime.strptime(utc_time, '%Y%m%d%H%M%S')
  35.     t_bj = t_utc + timedelta(hours=8)
  36.     ax.set_title('Channel01 ' + str(t_bj))
  37.     plt.savefig('D:/pic/satpy/c01-' + utc_time + '.png')
复制代码










密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2023-1-31 21:22:04 | 显示全部楼层
插个眼,以后有需要的时候来学习。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-2-1 01:48:46 | 显示全部楼层
谢谢分享!!!!!!!!!!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2023-2-11 13:51:08 | 显示全部楼层
正在学习读取风云4A
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-7-24 21:28:36 | 显示全部楼层
scn = Scene(filenames, reader='agri_l1') 一直识别不了是什么问题呢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-8-7 15:11:42 | 显示全部楼层
张nic 发表于 2023-7-24 21:28
scn = Scene(filenames, reader='agri_l1') 一直识别不了是什么问题呢

reader名字应该改一下,现在出了4a,4b,比如“agri_fy4b_l1”
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-8-7 17:08:08 | 显示全部楼层
太棒了!必须给楼主点个大大的赞!照着你的代码去画了FY-4B AGRI L1的全圆盘云图,成功啦!!{:5_213:}
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-8-9 10:12:31 | 显示全部楼层
谢谢分享!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2023-11-10 17:27:08 | 显示全部楼层
我用的是satpy0.44版本的,scn = Scene(filenames, reader='agri_l1'),reader='agri_fy4b_l1',都不能识别;https://satpy.readthedocs.io/en/stable/index.html#reader-table,这个网站里面的reader试了下都不识别,是怎么回事
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-11-10 17:29:29 | 显示全部楼层
之前我读fy4a的都没问题
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表