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'][:].dataValid_Value_Index = np.where(TPW_data != 65535)line_begin = InFile_REGC_GEO_4KM.variables['geospatial_lat_lon_extent'].begin_line_numberline_end = InFile_REGC_GEO_4KM.variables['geospatial_lat_lon_extent'].end_line_numbercolumn_begin = InFile_REGC_GEO_4KM.variables['geospatial_lat_lon_extent'].begin_pixel_numbercolumn_end = InFile_REGC_GEO_4KM.variables['geospatial_lat_lon_extent'].end_pixel_numberx = 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_Indexdef 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] #地球的长半轴，单位kmea = 6378.137#地球的短半轴，单位kmeb = 6356.7523#地心到卫星质心的距离，单位kmh = 42164#卫星星下点的经度Lon_Sate = 104.7#COFF全圆盘列偏移，4km对应的数值COFF = 1373.5#CFAC全圆盘列比例因子，4km对应的数值CFAC = 10233137#LOFF全圆盘行偏移，4km对应的数值LOFF = 1373.5#LFAC全圆盘列比例因子，4km对应的数值LFAC = 10233137tmp = 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 * cosytmp2 = 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 ) / tmp2s1 = h - sn * cosxys2 = sn * sinx * cosys3 = -sn * sinysxy = np.power(s1*s1 + s2*s2,0.5)tmp3 = 180 / np.pilon = tmp3 * np.arctan(s2/s1) + Lon_Satelat = 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[:,:] = 999Data_Lat[Valid_Value_Index] = lat return Data_Lon,Data_Latdef 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_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float64if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelif 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)

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

风云4A卫星行列号和经纬度查找表文件下载及解析
FY-4A是中国自主研发的静止气象卫星,其在太阳同步轨道上运行,每天可以全天候.全球范围内获取高空大气和地表的红外.可见光.微波等多种信息,对天气.气候研究及预报有着重要的作用. 一.风云卫星官网 ...
idea access数据库连接_idea代码神器：根据表生成代码
Easycode是idea的一个插件,可以直接对数据的表生成entity,controller,service,dao,mapper,无需任何编码,简单而强大. 1.安装(EasyCode) 我这里的 ...
sqlyog怎么查找表_VBA代码解决方案第58讲：在VBA中查找指定工作表的实用方法
大家好,我们今日继续讲解VBA代码解决方案的第58讲内容:在VBA中查找指定工作表的方法.在上一个例子中,我们通过一个自定义函数解决了删除工作表的方法.其实实现目的的方法有很多种,不必要必须有某种办法 ...
sin查找表 matlab,利用Xilinx中的ROM构造查找表来计算sin和cos的方法探讨
1.使用matlab制作.coe文件查找表的构造构造256点的正余弦表 exp(-j*2*pi*(0:255)/256),分别得到 cos和sin的查找表 matlab代码: 求sin fid = ...
SQL Server：查找表的生成或顺序
目录介绍背景临时表查找表关系详细信息查找表生成详细信息查找第一个生成表结果表生成表关系不需要的场景局限性 SqlServer_TableGeneration.zip - 2.6 ...
VS+openCV 用直方图统计像素（上）计算图像直方图、利用查找表修改图像外观
一.计算图像直方图图像由各种数值的像素构成.例如在单通道灰度图像中,每个像素都有一个 0(黑色)~255(白色)的整数.对于每个灰度,都有不同数量的像素分布在图像内,具体取决于图片内容. 直方图是一 ...
mybatis基础学习4-插件生成器(根据数据库的表生成文件)
1:安装(根据数据库的表生成文件) 2:在所建项目单击右键输入mybatis如下图 *建项目文件时不用建包和类,插件可以根据数据表自动生成,在配置文件(generatorConfig.xml)里写即可 ...
修改6S Fortran77 代码，建立查找表
逐像元大气校正,常预先计算查找表(LUT,LookUp Tabel),6S大气辐射传输模式也可以用来计算LUT.但6S源程序输出信息多,且浮点数输出精度低,不利于提取关键信息生成LUT,本文描述了怎样 ...
（18）FPGA面试题查找表的原理与结构
1.1 FPGA面试题查找表的原理与结构 1.1.1 本节目录 1)本节目录: 2)本节引言: 3)FPGA简介: 4)FPGA面试题查找表的原理与结构: 5)结束语. 1.1.2 本节引言 &quo ...

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

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

最新文章

热门文章