FY4A标称投影数据处理
MeteoInfo
编辑于 2024年10月08日 09:17
收录于文集
共51篇

风云四号气象卫星是我国第二代静止气象卫星，风云四号A星于2016年12月11日在西昌卫星发射中心成功发射。2017年9月25日，风云四号A星正式交付使用，标志着中国静止轨道气象卫星观测系统实现了更新换代。风云四号A星装载了先进的静止轨道辐射成像仪、干涉式大气垂直探测仪、闪电成像仪和空间环境监测仪器。

国家卫星气象中心发布的FY4A数据主要有NetCDF和HDF格式，投影方面主要是标称投影（NOM），也就是Geostationary投影，主要的投影参数有中央经度和卫星高度。NetCDF和HDF格式都是自描述的，除数据外可以放入很多辅助信息，遗憾的是FY4A数据中并没有包含完整的地理定位信息，如何读出数据对应的地理位置是首先要考虑的问题。官网提供了数据格点对应经纬度的查找表文件，可以用来对数据进行地理定位。

这里利用一个非常简洁的方法来处理地理定位问题，经过投影后的数据地理坐标单位为米，比如FY4A就有空间分辨率为4000米的数据，而且数据在投影后的坐标系中是规则的矩形，在全圆盘数据中，数据中心点坐标为0。利用上述特性可以快速重构出FY4A数据在Geostationary投影下的地理坐标，达到地理定位的目的。

先以一个全圆盘数据为例，MeteoInfoLab的addfile打开数据文件生成文件对象f：

代码块
Python
自动换行
复制代码
fn = 'D:/Temp/FY/FY4A-_AGRI--_N_DISK_1047E_L2-_CTT-_MULT_NOM_20190209140000_20190209141459_4000M_V0001.NC'
f = addfile(fn)
复制成功

在Console中查看文件对象的信息：

代码块
Python
自动换行
复制代码
>>> f
File Name: D:/Temp/FY/FY4A-_AGRI--_N_DISK_1047E_L2-_CTT-_MULT_NOM_20190209140000_20190209141459_4000M_V0001.NC
File type: Hierarchical Data Format, version 5 (NetCDF-4)
Dimensions: 3
  y = 2748;
  x = 2748;
  o = 1;
X Dimension: Xmin = 0.0; Xmax = 2747.0; Xsize = 2748; Xdelta = 1.0
Y Dimension: Ymin = 0.0; Ymax = 0.0; Ysize = 1; Ydelta = 1.0
Global Attributes: 
  : :dataset_name = "CTT"
  : :naming_authority = "NSMC CMA"
  : :Institution = "NSMC"
  : :Project = "NOM"
  : :Conventions = "CF-1.7"
  : :Metadata_Conventions = "Unidata Dataset Discovery v1.0"
  : :standard_name_vocabulary = "CF Standard Name Table (v25, 05 July 2013)"
  : :Title = "FY4A AGRI L2 Cloud Top Temperature"
  : :Summary = "Cloud Top Temperature"
  : :platform_ID = "FY4A"
  : :instrument_type = "FY4A Advanced Geosynchronous Radiation Imager"
  : :instrument_ID = "AGRI"
  : :processing_level = "L2"
  : :date_created = "2019-02-09T14:25:48Z"
  : :production_site = "NSMC"
  : :production_environment = "UNIX"
  : :Version_Of_Software = "V1.0"
  : :Software_Revision_Date = "2017-09-12"
  : :scene_id = "Full Disk possible values are Full Disk, Southern HEMisphere, the Northern HEMisphere, Regional, China Regional"
  : :spatial_resolution = "4km at nadir"
  : :time_coverage_start = "2019-02-09T14:00:04.568Z"
  : :time_coverage_end = "2019-02-09T14:12:42.257Z"
  : :Data_Quality = 0US
  : :L0QualityFlag = "0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000"
  : :PosQualityFlag = "0 0 0 0 0 0 0 0 0 0 0 0 0 0"
  : :CalQualityFlag = "0 0 0 0 0 0 0 0 0 0 0 0 0 0"
Variations: 11
  float y(y);
    y: :long_name = "FY4A fixed grid projection y-coordinate"
    y: :_ChunkSizes = 2748U
  float x(x);
    x: :long_name = "FY4A fixed grid projection x-coordinate"
    x: :_ChunkSizes = 2748U
  float CTT(x,y);
    CTT: :long_name = "FY4A PGS AGRI L2 Cloud Top Temperature"
    CTT: :standard_name = "Cloud Top Temperature"
    CTT: :_Unsigned = "TRUE"
    CTT: :FillValue = -999.0f
    CTT: :valid_range = 160.0f, 320.0f
    CTT: :scale_factor = 1.0f
    CTT: :add_offset = 0.0f
    CTT: :units = "K"
    CTT: :resolution = "4KM"
    CTT: :coordinates = "x y"
    CTT: :ancillary_variables = "DQF"
    CTT: :Description = "65535:Space, -999:FillValue"
    CTT: :_ChunkSizes = 916U, 916U
  byte DQF(x,y);
    DQF: :long_name = "CTT data quality flags"
    DQF: :standard_name = "status_flag"
    DQF: :_Unsigned = "TRUE"
    DQF: :FillValue = 127B
    DQF: :valid_range = 0B, 3B
    DQF: :units = "NULL"
    DQF: :coordinates = "x y"
    DQF: :flag_values = 0B, 1B, 2B, 3B
    DQF: :flag_meanings = "0:good_pixel,1:conditionally_usable_pixel,2:out_of_range_pixel,3:no_value_pixel"
    DQF: :number_of_qf_values = 4B
    DQF: :Description = 
    DQF: :_ChunkSizes = 1374U, 1374U
  float nominal_satellite_subpoint_lat(o);
    nominal_satellite_subpoint_lat: :long_name = "nominal satellite subpoint latitude (platform latitude)"
    nominal_satellite_subpoint_lat: :standard_name = "Latitude"
    nominal_satellite_subpoint_lat: :units = "degrees_north"
    nominal_satellite_subpoint_lat: :_ChunkSizes = 1U
  float nominal_satellite_subpoint_lon(o);
    nominal_satellite_subpoint_lon: :units = "degrees_east"
    nominal_satellite_subpoint_lon: :long_name = "nominal satellite subpoint longitude (platformlongitude)"
    nominal_satellite_subpoint_lon: :standard_name = "Longitude"
    nominal_satellite_subpoint_lon: :_ChunkSizes = 1U
  float nominal_satellite_height(o);
    nominal_satellite_height: :units = "km"
    nominal_satellite_height: :long_name = "nominal satellite height above GRS 80 ellipsoid(platform altitude)"
    nominal_satellite_height: :standard_name = "height_above_reference_ellipsoid"
    nominal_satellite_height: :_ChunkSizes = 1U
  float geospatial_lat_lon_extent(o);
    geospatial_lat_lon_extent: :long_name = "geospatial latitude and longitude references"
    geospatial_lat_lon_extent: :begin_line_number = 0US
    geospatial_lat_lon_extent: :end_line_number = 2747US
    geospatial_lat_lon_extent: :begin_pixel_number = 0US
    geospatial_lat_lon_extent: :end_pixel_number = 2747US
    geospatial_lat_lon_extent: :RegCenterLat = 65535.0f
    geospatial_lat_lon_extent: :RegCenterLon = 65535.0f
    geospatial_lat_lon_extent: :RegLength = 2748.0f
    geospatial_lat_lon_extent: :RegWidth = 2748.0f
    geospatial_lat_lon_extent: :geospatial_lat_units = "degrees_north"
    geospatial_lat_lon_extent: :geospatial_lon_units = "degrees_east"
    geospatial_lat_lon_extent: :_ChunkSizes = 1U
  int OBIType(o);
    OBIType: :OBIType_meanings = "0:Full_disk_observation 1:Southern_hemisphere_observation 2:Northern _hemisphere_observation 3:Regional observation"
    OBIType: :long_name = "Observing Type"
    OBIType: :standard_name = "OBIType"
    OBIType: :OBIType_values = 0, 1, 2, 3
    OBIType: :_ChunkSizes = 1U
  int processing_parm_version_container(o);
    processing_parm_version_container: :algorithm_version = "2016-10-10"
    processing_parm_version_container: :long_name = "container for processing parameter package filename and product version"
    processing_parm_version_container: :product_version = "2016-10-10"
    processing_parm_version_container: :_ChunkSizes = 1U
  int algorithm_product_version_container(o);
    algorithm_product_version_container: :long_name = "container for algorithm package filename and product version"
    algorithm_product_version_container: :product_version = "2016-10-10"
    algorithm_product_version_container: :algorithm_version = "2016-10-10"
    algorithm_product_version_container: :_ChunkSizes = 1U
复制成功

可以看出缺省读出的X Dimension和Y Dimension没有什么参考价值，CTT变量的维是(x, y)，这个写法也不太规范，实际数据的第一维是Y方向，第二维是X方向，写成(y, x)更合理。我们将文件中CTT变量的数据读出来，注意卫星数据通常Y方向是从上往下存储，与坐标系Y方向的顺序相反，因此读取数据时对Y轴进行反向处理，并将缺测数据处理为NaN。将Geostationary投影中的参数中央经度和卫星高度也读取出来：

代码块
Python
自动换行
复制代码
data = f['CTT'][::-1,:]
data[data>1000] = nan
data[data==-999] = nan
lon0 = f['nominal_satellite_subpoint_lon'][:]
height = f['nominal_satellite_height'][:]
复制成功

数据中X, Y方向的格点数均为2748，分辨率为4000米，数据中心位置坐标为(0, 0)。通过这些信息可以很简单的重建出在Geostationary投影中数据的地理坐标，X方向格点数是2748，那么x方向的总长度是(2748 - 1) * 4000，再除以2就可以得到X方向坐标起始和结束值为(-5494000.0,5494000.0)，Y方向也是类似。然后通过linspace函数生成X, Y坐标一维数组。

代码块
Python
自动换行
复制代码
ny, nx = data.shape
ex = (nx - 1) * 4000 / 2.
sx = -ex
ey = (ny - 1) * 4000 / 2.
sy = -ey
x = linspace(sx, ex, nx)
y = linspace(sy, ey, ny)
复制成功

可以通过绘图进行验证。利用projinfo函数创建Geostationary地图投影：

代码块
JavaScript
自动换行
复制代码
proj = projinfo(proj='geos', lon_0=lon0, h=height)
复制成功

然后创建该投影的地图坐标系，在将地图数据和卫星CTT数据绘制在地图坐标系中：

代码块
Python
自动换行
复制代码
ax = axesm(projection=proj, frameon=False, axison=False)
geoshow('coastline', color='k')
grid(True, tickvisible=True, tickposition='all')
xticks(arange(60, 151, 30))
yticks(arange(-60, 61, 30))
levs = arange(160, 311, 5)
layer = imshow(x, y, data, levs, proj=ax.proj, cmap='matlab_jet')
colorbar(layer, shrink=0.8, xshift=15)
title('FY4A Cloud Top Temperature')
复制成功

完整的脚本程序代码如下：

代码块
Python
自动换行
复制代码
fn = 'D:/Temp/FY/FY4A-_AGRI--_N_DISK_1047E_L2-_CTT-_MULT_NOM_20190209140000_20190209141459_4000M_V0001.NC'
f = addfile(fn)
data = f['CTT'][::-1,:]
data[data>1000] = nan
data[data==-999] = nan
lon0 = f['nominal_satellite_subpoint_lon'][:]
height = f['nominal_satellite_height'][:]
ny, nx = data.shape
ex = (nx - 1) * 4000 / 2.
sx = -ex
ey = (ny - 1) * 4000 / 2.
sy = -ey
x = linspace(sx, ex, nx)
y = linspace(sy, ey, ny)

#Plot
proj = projinfo(proj='geos', lon_0=lon0, h=height)
ax = axesm(projection=proj, frameon=False, axison=False)
geoshow('coastline', color='k')
grid(True, tickvisible=True, tickposition='all')
xticks(arange(60, 151, 30))
yticks(arange(-60, 61, 30))
levs = arange(160, 311, 5)
layer = imshow(x, y, data, levs, proj=ax.proj, cmap='matlab_jet')
colorbar(layer, shrink=0.8, xshift=15)
title('FY4A Cloud Top Temperature')
复制成功

生成的图形：

再看一个中国区域的FY4A标称投影数据：

代码块
Python
自动换行
复制代码
fn = 'D:/Temp/FY/FY4A-_AGRI--_N_REGC_1047E_L2-_TBB-_MULT_NOM_20190728123000_20190728123416_4000M_V0001.NC'
f = addfile(fn)
复制成功

数据文件信息如下：

代码块
Python
自动换行
复制代码
>>> f
File Name: D:/Temp/FY/FY4A-_AGRI--_N_REGC_1047E_L2-_TBB-_MULT_NOM_20190728123000_20190728123416_4000M_V0001.NC
File type: Hierarchical Data Format, version 5 (NetCDF-4)
Dimensions: 3
  o = 1;
  x = 1108;
  y = 2748;
X Dimension: Xmin = 0.0; Xmax = 1107.0; Xsize = 1108; Xdelta = 1.0
Y Dimension: Ymin = 0.0; Ymax = 0.0; Ysize = 2748; Ydelta = 1.0
Global Attributes: 
  : :dataset_name = "refer to filename conventions for L2 products in Appendix A."
  : :naming_authority = "NSMC CMA"
  : :Institution = "NSMC"
  : :Project = "NOM"
  : :Conventions = "CF-1.7"
  : :Metadata_Conv = "Unidata Dataset Discovery v1.0"
  : :standard_name_vocabul = "CF Standard Name Table (v25, 05 July 2013)"
  : :Title = "FY4A PGS L2 Blackbody Temperature"
  : :Summary = "Blackbody Temperature"
  : :platform_ID = "FY4A"
  : :instrument_type = "FY4A  Advanced Geosynchronous Radiation Image"
  : :instrument_ID = "AGRI"
  : :processing_level = "L2"
  : :date_created = "2019-07-28T12:40:48Z"
  : :production_site = "NSMC"
  : :production_environment = "UNIX"
  : :Version_Of_Software = "V1.0"
  : :Software_Revision_Date = "2017-09-12"
  : :scene_id = "Full Disk possible values are Full Disk,Southern HEMisphere,the Northern HEMisphere,Regional,China Regional"
  : :spatial_resolution = "4km at nadir."
  : :time_coverage_start = "2019-07-28T12:30:01.802Z"
  : :time_coverage_end = "2019-07-28T12:34:10.952Z"
  : :Data_Quality = 0US
  : :L0QualityFlag = "0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000"
  : :PosQualityFlag = "0 0 0 0 0 0 0 0 0 0 0 0 0 0"
  : :CalQualityFlag = "0 0 0 0 0 0 0 0 0 0 0 0 0 0"
Variations: 18
  float NOMChannel07(x,y);
    NOMChannel07: :long_name = "FY4A PGS L2 Blackbody Temperarure of 3.72um Channel(high)"
    NOMChannel07: :standard_name = "Blackbody Temperarure"
    NOMChannel07: :_Unsigned = "FALSE"
    NOMChannel07: :FillValue = 65534.0f
    NOMChannel07: :valid_range = 0.0f, 400.0f
    NOMChannel07: :scale_factor = 1.0f
    NOMChannel07: :add_offset = 0.0f
    NOMChannel07: :units = "K"
    NOMChannel07: :resolution = "4KM"
    NOMChannel07: :coordinates = "y x"
    NOMChannel07: :ancillary_variables = "DQF"
    NOMChannel07: :Description = "65535.0: space view; 65534.0: fill value"
    NOMChannel07: :_ChunkSizes = 554U, 1374U
  float NOMChannel08(x,y);
    NOMChannel08: :long_name = "FY4A PGS L2 Blackbody Temperarure of 3.72um Channel(low)"
    NOMChannel08: :standard_name = "Blackbody Temperarure"
    NOMChannel08: :_Unsigned = "FALSE"
    NOMChannel08: :FillValue = 65534.0f
    NOMChannel08: :valid_range = 0.0f, 400.0f
    NOMChannel08: :scale_factor = 1.0f
    NOMChannel08: :add_offset = 0.0f
    NOMChannel08: :units = "K"
    NOMChannel08: :resolution = "4KM"
    NOMChannel08: :coordinates = "y x"
    NOMChannel08: :ancillary_variables = "DQF"
    NOMChannel08: :Description = "65535.0: space view; 65534.0: fill value"
    NOMChannel08: :_ChunkSizes = 554U, 1374U
  float NOMChannel09(x,y);
    NOMChannel09: :long_name = "FY4A PGS L2 Blackbody Temperarure of 6.25um Channel"
    NOMChannel09: :standard_name = "Blackbody Temperarure"
    NOMChannel09: :_Unsigned = "FALSE"
    NOMChannel09: :FillValue = 65534.0f
    NOMChannel09: :valid_range = 0.0f, 400.0f
    NOMChannel09: :scale_factor = 1.0f
    NOMChannel09: :add_offset = 0.0f
    NOMChannel09: :units = "K"
    NOMChannel09: :resolution = "4KM"
    NOMChannel09: :coordinates = "y x"
    NOMChannel09: :ancillary_variables = "DQF"
    NOMChannel09: :Description = "65535.0: space view; 65534.0: fill value"
    NOMChannel09: :_ChunkSizes = 554U, 1374U
  float NOMChannel10(x,y);
    NOMChannel10: :long_name = "FY4A PGS L2 Blackbody Temperarure of 7.10um Channel"
    NOMChannel10: :standard_name = "Blackbody Temperarure"
    NOMChannel10: :_Unsigned = "FALSE"
    NOMChannel10: :FillValue = 65534.0f
    NOMChannel10: :valid_range = 0.0f, 400.0f
    NOMChannel10: :scale_factor = 1.0f
    NOMChannel10: :add_offset = 0.0f
    NOMChannel10: :units = "K"
    NOMChannel10: :resolution = "4KM"
    NOMChannel10: :coordinates = "y x"
    NOMChannel10: :ancillary_variables = "DQF"
    NOMChannel10: :Description = "65535.0: space view; 65534.0: fill value"
    NOMChannel10: :_ChunkSizes = 554U, 1374U
  float NOMChannel11(x,y);
    NOMChannel11: :long_name = "FY4A PGS L2 Blackbody Temperarure of 8.50um Channel"
    NOMChannel11: :standard_name = "Blackbody Temperarure"
    NOMChannel11: :_Unsigned = "FALSE"
    NOMChannel11: :FillValue = 65534.0f
    NOMChannel11: :valid_range = 0.0f, 400.0f
    NOMChannel11: :scale_factor = 1.0f
    NOMChannel11: :add_offset = 0.0f
    NOMChannel11: :units = "K"
    NOMChannel11: :resolution = "4KM"
    NOMChannel11: :coordinates = "y x"
    NOMChannel11: :ancillary_variables = "DQF"
    NOMChannel11: :Description = "65535.0: space view; 65534.0: fill value"
    NOMChannel11: :_ChunkSizes = 554U, 1374U
  float NOMChannel12(x,y);
    NOMChannel12: :long_name = "FY4A PGS L2 Blackbody Temperarure of 10.8um Channel"
    NOMChannel12: :standard_name = "Blackbody Temperarure"
    NOMChannel12: :_Unsigned = "FALSE"
    NOMChannel12: :FillValue = 65534.0f
    NOMChannel12: :valid_range = 0.0f, 400.0f
    NOMChannel12: :scale_factor = 1.0f
    NOMChannel12: :add_offset = 0.0f
    NOMChannel12: :units = "K"
    NOMChannel12: :resolution = "4KM"
    NOMChannel12: :coordinates = "y x"
    NOMChannel12: :ancillary_variables = "DQF"
    NOMChannel12: :Description = "65535.0: space view; 65534.0: fill value"
    NOMChannel12: :_ChunkSizes = 554U, 1374U
  float NOMChannel13(x,y);
    NOMChannel13: :long_name = "FY4A PGS L2 Blackbody Temperarure of 12um Channel"
    NOMChannel13: :standard_name = "Blackbody Temperarure"
    NOMChannel13: :_Unsigned = "FALSE"
    NOMChannel13: :FillValue = 65534.0f
    NOMChannel13: :valid_range = 0.0f, 400.0f
    NOMChannel13: :scale_factor = 1.0f
    NOMChannel13: :add_offset = 0.0f
    NOMChannel13: :units = "K"
    NOMChannel13: :resolution = "4KM"
    NOMChannel13: :coordinates = "y x"
    NOMChannel13: :ancillary_variables = "DQF"
    NOMChannel13: :Description = "65535.0: space view; 65534.0: fill value"
    NOMChannel13: :_ChunkSizes = 554U, 1374U
  float NOMChannel14(x,y);
    NOMChannel14: :long_name = "FY4A PGS L2 Blackbody Temperarure of 13.5um Channel"
    NOMChannel14: :standard_name = "Blackbody Temperarure"
    NOMChannel14: :_Unsigned = "FALSE"
    NOMChannel14: :FillValue = 65534.0f
    NOMChannel14: :valid_range = 0.0f, 400.0f
    NOMChannel14: :scale_factor = 1.0f
    NOMChannel14: :add_offset = 0.0f
    NOMChannel14: :units = "K"
    NOMChannel14: :resolution = "4KM"
    NOMChannel14: :coordinates = "y x"
    NOMChannel14: :ancillary_variables = "DQF"
    NOMChannel14: :Description = "65535.0: space view; 65534.0: fill value"
    NOMChannel14: :_ChunkSizes = 554U, 1374U
  float nominal_satellite_subpoint_lat(o);
    nominal_satellite_subpoint_lat: :long_name = "nominal satellite subpoint latitude (platform latitude)"
    nominal_satellite_subpoint_lat: :standard_name = "Latitude"
    nominal_satellite_subpoint_lat: :units = "degrees_north"
    nominal_satellite_subpoint_lat: :_ChunkSizes = 1U
  float nominal_satellite_subpoint_lon(o);
    nominal_satellite_subpoint_lon: :units = "degrees_east"
    nominal_satellite_subpoint_lon: :long_name = "nominal satellite subpoint longitude (platformlongitude)"
    nominal_satellite_subpoint_lon: :standard_name = "Longitude"
    nominal_satellite_subpoint_lon: :_ChunkSizes = 1U
  float nominal_satellite_height(o);
    nominal_satellite_height: :units = "km"
    nominal_satellite_height: :long_name = "nominal satellite height above GRS 80 ellipsoid(platform altitude)"
    nominal_satellite_height: :standard_name = "height_above_reference_ellipsoid"
    nominal_satellite_height: :_ChunkSizes = 1U
  float geospatial_lat_lon_extent(o);
    geospatial_lat_lon_extent: :long_name = "geospatial latitude and longitude references"
    geospatial_lat_lon_extent: :begin_line_number = 179US
    geospatial_lat_lon_extent: :end_line_number = 1286US
    geospatial_lat_lon_extent: :begin_pixel_number = 0US
    geospatial_lat_lon_extent: :end_pixel_number = 2747US
    geospatial_lat_lon_extent: :RegCenterLat = 65535.0f
    geospatial_lat_lon_extent: :RegCenterLon = 65535.0f
    geospatial_lat_lon_extent: :RegLength = 1108.0f
    geospatial_lat_lon_extent: :RegWidth = 2748.0f
    geospatial_lat_lon_extent: :geospatial_lat_units = "degrees_north"
    geospatial_lat_lon_extent: :geospatial_lon_units = "degrees_east"
    geospatial_lat_lon_extent: :_ChunkSizes = 1U
  int OBIType(o);
    OBIType: :OBIType_meanings = "0:Full_disk_observation 1:Southern_hemisphere_observation 2:Northern _hemisphere_observation 3:Regional observation"
    OBIType: :long_name = "Observing Type"
    OBIType: :standard_name = "OBIType"
    OBIType: :OBIType_values = 0, 1, 2, 3
    OBIType: :_ChunkSizes = 1U
  int processing_parm_version_container(o);
    processing_parm_version_container: :long_name = "container for processing parameter package filename and product version"
    processing_parm_version_container: :product_version = "2016-10-10"
    processing_parm_version_container: :processing_parm_version = "2016-10-10"
    processing_parm_version_container: :_ChunkSizes = 1U
  int algorithm_product_version_container(o);
    algorithm_product_version_container: :long_name = "container for algorithm package filename and product version"
    algorithm_product_version_container: :product_version = "2016-10-10"
    algorithm_product_version_container: :algorithm_version = "2016-10-10"
    algorithm_product_version_container: :_ChunkSizes = 1U
  byte DQF(x,y);
    DQF: :long_name = "TBB data quality flags"
    DQF: :standard_name = "status_flag"
    DQF: :_Unsigned = "FALSE"
    DQF: :FillValue = 127B
    DQF: :valid_range = 0B, 3B
    DQF: :units = "NULL"
    DQF: :coordinates = "y x"
    DQF: :flag_values = 0B, 1B, 2B, 3B
    DQF: :flag_meanings = "good_pixel conditionally_usable_pixel out_of_range_pixel no_value_pixel"
    DQF: :number_of_qf_values = 4B
    DQF: :Description = "0:good_pixel; 1:conditionally_usable_pixel; 2:out_of_range_pixel; 3:no_value_pixel"
    DQF: :_ChunkSizes = 1108U, 2748U
  float x(x);
    x: :long_name = "FY4A fixed grid projection x-coordinate"
    x: :_ChunkSizes = 1108U
  float y(y);
    y: :long_name = "FY4A fixed grid projection y-coordinate"
    y: :_ChunkSizes = 2748U
复制成功

可以看出数据的X, Y维长度分别为1108和2748，同样的问题，X, Y维的名称其实是写反了。实际数据X方向的格点数是2748，Y方向的格点数是1108。读取数据变量NOMChannel07看看读出来的二维数组的shape属性就可以验证出来（同样需要注意读取数组时第一维Y维需要反向）：

代码块
Python
自动换行
复制代码
data = f['NOMChannel07'][::-1,:]
data[data>65500] = nan
复制成功

代码块
Python
自动换行
复制代码
>>> data.shape
(1108, 2748)
复制成功

数据的空间分辨率也是4000米，所以X方向的数据坐标和全圆盘数据是一样的。问题是怎么来确定Y方向的数据坐标？可以用geospatial_lat_lon_extent变量中的属性信息，首先获取文件对象中该变量对象llextent，然后从变量对象中读取end_line_number属性值：

代码块
Python
自动换行
复制代码
llextent = f['geospatial_lat_lon_extent']
eline = llextent.attrvalue('end_line_number')[0]
复制成功

代码块
Python
自动换行
复制代码
>>> eline
1286
复制成功

begin_line_number和end_line_number是指区域格点数据Y方向第一行和最后一行数据在全圆盘数据中的位置，注意这个数字是Y方向从上往下来数的。我们知道全圆盘数据Y方向的长度是2748，数据Y坐标是从下往上的顺序，而全圆盘数据的Y坐标我们已经知道如何计算了，那么只需要计算出区域数据最南边一行在全圆盘Y坐标中是第几个就可以了。begin_pixel_number和end_pixel_number是X方向区域数据在全圆盘数据中的起始结束位置。区域数据地理坐标的计算代码如下：

代码块
Python
自动换行
复制代码
ny, nx = data.shape
sx = -(4000 * (nx - 1)) / 2.
sy = sx + (nx - eline - 1) * 4000
y = arange1(sy, ny, 4000)
x = arange1(sx, nx, 4000)
复制成功

其他的代码就很类似了。完整的代码如下：

代码块
Python
自动换行
复制代码
fn = 'D:/Temp/FY/FY4A-_AGRI--_N_REGC_1047E_L2-_TBB-_MULT_NOM_20190728123000_20190728123416_4000M_V0001.NC'
f = addfile(fn)
data = f['NOMChannel07'][::-1,:]
data[data>65500] = nan
lon0 = f['nominal_satellite_subpoint_lon'][:]
height = f['nominal_satellite_height'][:]
llextent = f['geospatial_lat_lon_extent']
eline = llextent.attrvalue('end_line_number')[0]
ny, nx = data.shape
sx = -(4000 * (nx - 1)) / 2.
sy = sx + (nx - eline - 1) * 4000
y = arange1(sy, ny, 4000)
x = arange1(sx, nx, 4000)

#Plot
proj = geolib.Geostationary(central_longitude=lon0, satellite_height=height)
ax = axesm(projection=proj, axison=False, frameon=False)
geoshow('cn_province', edgecolor='gray')
geoshow('country', edgecolor='k')
grid(True, tickvisible=True, tickposition='all')
xticks(arange(60, 151, 30))
yticks(arange(-60, 61, 30))
imshow(x, y, data, transform=proj, fill_color='lightgray')
colorbar(shrink=0.8, xshift=15)
axis()
复制成功

绘图结果如下：

MeteoInfoLab的geolib包中有一些投影函数，可以将Geostationary投影下的地理坐标转换为经纬度，projinfo函数不带任何参数返回等经纬度投影。

代码块
Python
自动换行
复制代码
>>> projinfo()
+title=long/lat:WGS84 +proj=longlat +ellps=WGS84 +datum=WGS84 +units=degrees
复制成功

geolib的project函数进行不同投影下的坐标转换，转为经纬度后在等经纬度投影坐标系下数据就不是规则的矩形网格，可以用pcolor函数绘图（速度较慢）。

代码块
Python
自动换行
复制代码
fn = 'D:/Temp/FY/FY4A-_AGRI--_N_REGC_1047E_L2-_TBB-_MULT_NOM_20190728123000_20190728123416_4000M_V0001.NC'
f = addfile(fn)
data = f['NOMChannel07'][::-1,:]
data[data>65500] = nan
lon0 = f['nominal_satellite_subpoint_lon'][:]
height = f['nominal_satellite_height'][:]
llextent = f['geospatial_lat_lon_extent']
eline = llextent.attrvalue('end_line_number')[0]
ny, nx = data.shape
sx = -(4000 * (nx - 1)) / 2.
sy = sx + (nx - eline - 1) * 4000
y = arange1(sy, ny, 4000)
x = arange1(sx, nx, 4000)

#project data from geostationary to long/lat
from_proj = geolib.Geostationary(central_longitude=lon0, satellite_height=height)
to_proj = geolib.PlateCarree()
xx, yy = meshgrid(x, y)
lon, lat = geolib.project(xx, yy, fromproj=from_proj, toproj=to_proj)
lon[lon<0] = nan

#Plot
ax = axesm()
geoshow('cn_province', edgecolor='gray')
geoshow('country', edgecolor='k')
grid(True)
pcolor(lon, lat, data, zorder=0)
colorbar(shrink=0.8)
复制成功

图形结果：

也可以通过geolib的reproject函数将Geostationary投影下的原始规则矩形网格数据投影到等经纬度，并保持规则矩形：

代码块
Python
自动换行
复制代码
fn = 'D:/Temp/FY/FY4A-_AGRI--_N_REGC_1047E_L2-_TBB-_MULT_NOM_20190728123000_20190728123416_4000M_V0001.NC'
f = addfile(fn)
data = f['NOMChannel07'][::-1,:]
data[data>65500] = nan
lon0 = f['nominal_satellite_subpoint_lon'][:]
height = f['nominal_satellite_height'][:]
llextent = f['geospatial_lat_lon_extent']
eline = llextent.attrvalue('end_line_number')[0]
ny, nx = data.shape
sx = -(4000 * (nx - 1)) / 2.
sy = sx + (nx - eline - 1) * 4000
y = arange1(sy, ny, 4000)
x = arange1(sx, nx, 4000)

#project data from geostationary to long/lat
from_proj = geolib.Geostationary(central_longitude=lon0, satellite_height=height)
to_proj = geolib.PlateCarree()
lon = arange(20, 190.1, 0.1)
lat = arange(0, 60.1, 0.1)
pdata = geolib.reproject(data, x, y, from_proj, lon, lat, to_proj)

#Plot
ax = axesm()
geoshow('cn_province', edgecolor='gray')
geoshow('country', edgecolor='k')
grid()
imshow(lon, lat, pdata, fill_color='lightgray', zorder=0)
colorbar(shrink=0.8)
复制成功

图形结果：

cv30670541