FY2G经纬度

nlj917

我已获取到风云2G降水数据（hdf），2288*1144矩阵，请问如何得到各栅格对应的经纬度坐标？

edwardli

# -*- coding: utf-8 -*-
"""
Created on Wed Sep  1 15:31:46 2021

@author: OLDLee
"""

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from projection import latlon2linecolumn
import cartopy.crs as ccrs

## 打开
sat = xr.open_dataset('data/Satellite/Z_SATE_C_BAWX_20201205125337_P_FY4A-_AGRI--_N_REGC.NC')

## 从人家封装好的方法里摘出来的，（有的可能没用的）属性们
resolution = '4000M'
line_begin = sat.geospatial_lat_lon_extent.attrs['begin_line_number']
line_end = sat.geospatial_lat_lon_extent.attrs['end_line_number']
column_begin = sat.geospatial_lat_lon_extent.attrs['begin_pixel_number']
column_end = sat.geospatial_lat_lon_extent.attrs['end_pixel_number']

## 这个geo_desc很重要，决定了最终导出哪个经纬度范围、分辨率的数据
geo_desc = [0, 60, 50, 150, 0.1]
lat_S, lat_N, lon_W, lon_E, step = [1000 * x for x in geo_desc]
lat = np.arange(lat_N, lat_S-1, -step) / 1000
lon = np.arange(lon_W, lon_E+1, step) / 1000
lon_mesh, lat_mesh = np.meshgrid(lon, lat)
## 求geo_desc对应的标称全圆盘行列号。。。
## ！！！！这个projection.latlon2linecolumn方法就是按照NSMC上的行列号查找表/方法写的，我大致比对过
line, column = latlon2linecolumn(lat_mesh, lon_mesh, resolution)
line = xr.DataArray(line, coords=(('lat', lat), ('lon', lon)), name='line')
column = xr.DataArray(column, coords=(('lat', lat), ('lon', lon)), name='column')


## “构造”DataArray
dn_values = sat['Precipitation']
dn_values = dn_values.rename({dn_values.dims[0]: 'line', dn_values.dims[1]: 'column'})
dn_values = dn_values.assign_coords(line=range(line_begin, line_end+1),
                                    column=range(column_begin, column_end+1))

if geo_desc:
    # 若geo_desc已指定，则插值到对应网格
    dn_values = dn_values.interp(line=line, column=column, method='nearest')
    del dn_values.coords['line'], dn_values.coords['column']
else:
    # 若geo_desc为None，则保留原始NOM网格
    pass

##  后面比较乱了，也没剔除65534空白值…………凑合看吧

# fig =plt.figure(figsize=(9,6),dpi=200)
# ax = plt.subplot(111,projection=ccrs.PlateCarree())
# plt.contourf(column.lon,line.lat,dn_values,cmap='jet',levels=np.arange(0,30,1),transform=ccrs.PlateCarree())
# plt.show()
da = dn_values
p = da.plot(
    subplot_kws=dict(projection=ccrs.PlateCarree(), facecolor="gray"),
    transform=ccrs.PlateCarree(),
    levels=[0, 0.1],
    # size=3
)

p.axes.set_global()
p.axes.coastlines()
p.axes.gridlines(draw_labels=True)

plt.show()

二爷名声在外

hdf描述里面有吧，我记得我读过……

nlj917

edwardli 发表于 2022-3-13 10:31
# -*- coding: utf-8 -*-
"""
Created on Wed Sep  1 15:31:46 2021

好的，谢谢

nlj917

edwardli 发表于 2022-3-13 10:31
# -*- coding: utf-8 -*-
"""
Created on Wed Sep  1 15:31:46 2021

您好，请问您那边有没有FY产品各网格的经纬度？我在进行评估，需要经纬度栅格坐标，可是一直获取不到