FY-4A/LPW产品4km经纬度查找表生成代码-风云四号

实在是太忙，没时间整理，没办法，抽时间整理一下吧。
FY-4A/LPW水汽产品4km经纬度查找表生成代码：

# -*- coding: utf-8 -*-
"""
Created on Tue Apr 23 11:39:02 2019

@author: Administrator
"""


import os
import numpy as np 
import math
import netCDF4 as nc
import h5py 
from osgeo import gdal
import operator


#读取FY4A-_AGRI--_N_DISK_1047E_L1-_GEO-_MULT_NOM文件，空间分辨率4km
#此程序仅限于此类文件

#该文件中仅用于获得FY4A-_AGRI--_N_REGC_1047E_L1_***_4000M_***.HDF文件
#中数据的行列号转换成经纬度

def Get_Line_Col_in_Disk(File_Full_Path):
    #print('#########################################')
    #判读输入文件是否为FY4A-_AGRI--_N_REGC_1047E_L1_***_4000M_***.HDF文件
    #FY4A-_AGRI--_N_REGC_1047E_L2-_LPW-_MULT_NOM_20210303074918_20210303075335_4000M_V0001.NC



    #以只读的方式打开hdf5文件
    InFile_REGC_GEO_4KM = nc.Dataset(File_Full_Path,'r')
    #data.variables
    #外围无效列标
    TPW_data = InFile_REGC_GEO_4KM.variables['TPW'][:].data
    Valid_Value_Index = np.where(TPW_data != 65535)
    
    line_begin = InFile_REGC_GEO_4KM.variables['geospatial_lat_lon_extent'].begin_line_number
    line_end = InFile_REGC_GEO_4KM.variables['geospatial_lat_lon_extent'].end_line_number
    column_begin = InFile_REGC_GEO_4KM.variables['geospatial_lat_lon_extent'].begin_pixel_number
    column_end = InFile_REGC_GEO_4KM.variables['geospatial_lat_lon_extent'].end_pixel_number
    x = np.array(list(range(column_begin,column_end+1)))
    y = np.array(list(range(line_begin,line_end+1)))
    X,Y = np.meshgrid(x,y)

    
    InFile_REGC_GEO_4KM.close() 
    return Y,X,Valid_Value_Index
    
def Convert_Line_Col_To_Lon_Lat(Data_Line_Num_in_Disk,Data_Col_Num_in_Disk,Valid_Value_Index):
    #print('##########################################')  
    ##根据中国区域在全圆盘中的行列号计算经纬度
    #注意：各个文件中在全圆盘的行列号会有所不同
    #中国区域内的有效行列号
           

    Line = Data_Line_Num_in_Disk[Valid_Value_Index]     
    Column = Data_Col_Num_in_Disk[Valid_Value_Index] 
     
   
    #地球的长半轴，单位km
    ea = 6378.137
    #地球的短半轴，单位km
    eb = 6356.7523
    #地心到卫星质心的距离，单位km
    h = 42164
    #卫星星下点的经度
    Lon_Sate = 104.7
    #COFF全圆盘列偏移，4km对应的数值
    COFF = 1373.5
    #CFAC全圆盘列比例因子，4km对应的数值
    CFAC = 10233137
    #LOFF全圆盘行偏移，4km对应的数值
    LOFF = 1373.5
    #LFAC全圆盘列比例因子，4km对应的数值
    LFAC = 10233137
    
    tmp = np.pi/(180*math.pow(2,-16))
    #print('tmp =',tmp)
    
    x = tmp * (Column - COFF)/CFAC
    #print('x =',x)
    
    y = tmp * (Line - LOFF)/LFAC
    #print('y =',y)

    sinx = np.sin(x)
    cosx = np.cos(x)
    siny = np.sin(y)
    cosy = np.cos(y)
    
    cosxy = cosx * cosy
    tmp2 = cosy*cosy + siny*siny * ea*ea /(eb*eb) 
    #print('tmp2=',tmp)
    sd = np.power((h*cosxy)*(h*cosxy) - tmp2*(h*h-ea*ea),0.5)
    sn = ( h*cosxy - sd ) / tmp2
    s1 = h - sn * cosxy
    s2 = sn * sinx * cosy
    s3 = -sn * siny
    sxy = np.power(s1*s1 + s2*s2,0.5)
    
    tmp3 = 180 / np.pi
    lon = tmp3 * np.arctan(s2/s1) + Lon_Sate
    lat = tmp3 * np.arctan((ea*ea*s3)/(eb*eb*sxy))
    #print('lon =',lon)
    #print('lat =',lat)
    print(lon)
    Data_Lon = np.zeros(Data_Line_Num_in_Disk.shape, dtype = np.float64,order = 'C' )
    Data_Lon[:,:] = 999    
    Data_Lon[Valid_Value_Index] = lon
    
    #print(Data_Lon)
    Data_Lat = np.zeros(Data_Line_Num_in_Disk.shape, dtype = np.float64,order = 'C' )
    Data_Lat[:,:] = 999
    Data_Lat[Valid_Value_Index] = lat 
    
    return Data_Lon,Data_Lat
    
def writeTiff(im_data,im_width,im_height,im_bands,path):#im_geotrans,im_proj,path):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float64

    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    elif len(im_data.shape) == 2:
        im_data = np.array([im_data])
    else:
        im_bands, (im_height, im_width) = 1,im_data.shape
        #创建文件
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(path, im_width, im_height, im_bands, datatype)
    #if(dataset!= None):
        #dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数
        #dataset.SetProjection(im_proj) #写入投影
    for i in range(im_bands):
        dataset.GetRasterBand(i+1).WriteArray(im_data[i])
    del dataset
     

#主程序
if __name__ =='__main__': 
    
    GEO_File_Full_Path = 'E:\FY4A-_AGRI--_N_REGC_1047E_L2-_LPW-_MULT_NOM_20210303074918_20210303075335_4000M_V0001.NC'   
    GEO_file_name = os.path.basename(GEO_File_Full_Path)
    first = GEO_file_name[44:57]
    #print(first)
    second = GEO_file_name[73:79]
    #print(second)
    out_name = first + second + '经纬度查找表'
  
    Data_line_col = Get_Line_Col_in_Disk(GEO_File_Full_Path)      
    #print(Data_line_col)
    Lon_Lat = Convert_Line_Col_To_Lon_Lat(Data_line_col[0],Data_line_col[1],Data_line_col[2])

       
    #print(type(Lon_Lat))
    if Lon_Lat is None :
        print ('wrong file')
    #elif (id(GEO_file_name[44:73]) is not id(Ocean_file_name[44:73])):
        #print(type(GEO_file_name[44:73]))
        #print('Geo file does not match Ocean file')
    else:
        
        Lon = Lon_Lat[0] 
        Lat = Lon_Lat[1]        
    
    im_data = np.array([Lon,Lat])
    im_bands, im_height, im_width = im_data.shape
    
    #保存tif文件函数
    path = 'E:\\' + out_name + '.tif'
    
    writeTiff(im_data, im_width, im_height,im_bands,path)