python 风云4号云顶温度处理

根据网上的一些关于风云4号卫星的数据处理方法与本人自己的修改

风云4号卫星云顶温度出图代码如下：

对于色斑的选择下一步还需要继续修改和精进

#!usr/bin/env python
# -*- coding:utf-8 -*-
"""
@author: Suyue
@file: CTT_map.py
@time: 2024/06/26
@desc:
"""
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib as mpl
from pyresample import geometry, image
import cartopy.feature as cfeature
import cartopy.crs as ccrs
import warnings
warnings.filterwarnings('ignore')



# 读取FY4A的4km全圆盘的经纬度对照表
def get_lonlat_fulldisk(tbl):
    dim = 2748
    sz = np.fromfile(tbl,dtype=float,count=dim*dim*2)
    latlon = np.reshape(sz,(dim,dim,2))
    lats_tbl = latlon[:,:,0]
    lons_tbl = latlon[:,:,1]
    return lons_tbl, lats_tbl

# 使用geometry生成FY4A的GEOS投影
def get_area_fulldisk():
    # NOMCCenterLon标称数据中心经度
    sublon = 104.7
    # NOMSatHeight标称卫星高度
    height = 35786000
    # dEA*1000地球赤道半径
    long_axis = 6378140
    # d0bRecFlat地球扁率的倒数
    rf = 298.257223563
    area_id = 'FullDisk'
    description = 'geos:4KM'
    proj_id = 'geos'
    proj_dict = {'proj':'geos','lon_0':sublon,'h':height,'a':long_axis,'rf':rf,'units':'m'}
    width = 2748
    height = 2748
    radius = 5496005.4
    area_extent = (-radius,-radius,radius,radius)
    area_def_fulldisk = geometry.AreaDefinition(area_id,proj_id,description,proj_dict,width,height,area_extent)
    return area_def_fulldisk

# 检查生成的GEOS投影与经纬度对照表的经纬度是否相近
def check_lation_fulldisk(tbl,area_def_fulldisk):
    lons_fd_test,lats_fd_test = area_def_fulldisk.get_lonlats()
    lons_tbl,lats_tbl = get_lonlat_fulldisk(tbl)
    lons_fd_test = np.where(lons_tbl==999999.9999,1.e30,lons_tbl)
    lats_fd_test = np.where(lats_tbl==999999.9999,1.e30,lats_tbl)
    # ~50m
    assert (np.abs(lons_fd_test-lons_fd_test)<0.0005).all()
    # ~50m
    assert (np.abs(lats_fd_test-lats_fd_test)<0.005).all()
    return

# 定义关注区域的投影方式，这里以兰伯特投影为例进行说明
def get_area_china():
    area_id = 'China'
    description = 'lambert:10KM'
    proj_id = 'lambert'
    proj_dict = {'proj':'lcc','lat_0':35,'lon_0':105,'lat_1':20,'lat_2':50,'units':'m'}
    width = 803
    height = 603
    lenx,leny = 4.015E6,3.015E6
    area_extent = (-lenx,leny,lenx,-leny)
    area_def_china = geometry.AreaDefinition(area_id,proj_id,description,proj_dict,width,height,area_extent)
    return area_def_china

# 读取FY4A的L2级产品，这里以CTT为例进行说明
def get_data_ctt(filename_ctt_l2):
    fy4a_l2 = xr.open_dataset(filename_ctt_l2)
    # space value取CTT值
    ctt = np.ma.masked_equal(fy4a_l2['CTT'].values,65535)
    # fill value填充值
    ctt = np.ma.masked_equal(ctt,-1)
    ctt = xr.DataArray(ctt,dims=('y','x'))
    fulldisk = get_area_fulldisk()
    x2d,y2d = fulldisk.get_proj_coords()
    ctt.coords['x'],ctt.coords['y'] = x2d[0,:],y2d[:,0]
    print('CTT:min value is {},max value is {}'.format(np.nanmin(ctt),np.nanmax(ctt)))
    return ctt

filename_ctt_l2 = 'G:/Z_SATE_C_BAWX_20230703004255_P_FY4A-_AGRI--_N_DISK_1047E_L2-_CTT-_MULT_NOM_20230703001500_20230703002959_4000M_V0001.NC'
tbl = 'G:/FullMask_Grid_4000.raw'
check_latlon = True
date = '2023-07-3_00:15:00'

if check_latlon:
    area_def_fulldisk = get_area_fulldisk()
    check_lation_fulldisk(tbl, area_def_fulldisk)
ctt = get_data_ctt(filename_ctt_l2)

### FY4A L2绘图函数
# data: 绘图数据
# crs_from/crs_to: 数据投影方式/图片投影方式
# nrow/ncol/index: Panel中图片的总行数/总列数/图片位置
# bounds: 等值线
# title: 图片标题
# cbar/unit: 是否展示colorbar/colorbar名称
def plot_var(data, ctt_from, ctt_to, figure, nrow, ncol, index, bounds, title, cbar, unit):
    ax = figure.add_subplot(nrow,ncol,index,projection=ctt_to)
    ### please define your own cmap
    colors = np.array([[40,84,255],[71,169,255],[128,224,255],[191,250,255],[255,237,169],
                       [255,188,125],[248,122,98],[220,47,55]])/255.
    cmap = mpl.colors.ListedColormap(colors).with_extremes(under=np.array([41,10,216])/255.,
                                                           over=np.array([165,0,33])/255.)
    assert len(bounds)==cmap.N+1
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

    img = data.plot.imshow(cmap=cmap, norm=norm, extend='both', ax=ax, origin='upper', transform=ctt_from, extent=ctt_from.bounds, add_colorbar=False)
    bodr = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_boundary_lines_land', scale='10m', facecolor='none')
    ax.add_feature(bodr, edgecolor='k', alpha=1, linewidth=1.) #添加国界线
    ax.coastlines(resolution='10m', alpha=1, linewidth=1.) #添加海岸线
    ax.gridlines(linestyle='--', alpha=0.5, linewidth=1., color='k')
    ax.set_title(title, fontsize=16)
    plt.axis('off')
    if cbar:
        clb = plt.colorbar(img, ax=ax, orientation="vertical", pad=0.03, shrink=0.6)
        clb.ax.set_title(unit)
    return ax


figure = plt.figure(figsize=(12,12))
bounds = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
fulldisk = get_area_fulldisk()
plot_var(ctt, fulldisk.to_cartopy_crs(), fulldisk.to_cartopy_crs(), figure, 1, 1, 1, bounds, 'FY4A L2 Cloud Fraction', True, '')

plt.savefig('G:/pic.jpg')