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')