321 lines
13 KiB
Python
321 lines
13 KiB
Python
"""
|
||
@Project ：microproduct
|
||
@File ：AtmosphericAuxData.PY
|
||
@Function ：辅助文件读写
|
||
@Author ：LMM
|
||
@Date ：2021/09/3 14:39
|
||
@Version ：1.0.0
|
||
"""
|
||
import os
|
||
import shutil
|
||
import cftime
|
||
import numpy as np
|
||
import netCDF4 as Nc
|
||
import math
|
||
from tool.algorithm.image.ImageHandle import ImageHandler
|
||
from osgeo import gdal, osr
|
||
import gc
|
||
import logging
|
||
|
||
logger = logging.getLogger("mylog")
|
||
|
||
|
||
class NcHandle:
|
||
"""
|
||
气象数据的检查
|
||
气象数据的读取、降维、转tif
|
||
"""
|
||
def __int__(self):
|
||
pass
|
||
|
||
@staticmethod
|
||
def check_nc(nc_infile):
|
||
"""
|
||
检查气象数据
|
||
nc_infile： 气象数据路径
|
||
"""
|
||
nc_data_obj = Nc.Dataset(nc_infile)
|
||
keys = nc_data_obj.variables.keys() # 查看气象数据文件包含的变量列表
|
||
inform = nc_data_obj.variables
|
||
|
||
variables = ['longitude', 'latitude', 'level', 'time', 'z', 'r', 't']
|
||
|
||
if len(keys) < len(variables): # 检查缺少了哪种气象数据
|
||
logger.error("下载气象数据有误，缺少参数")
|
||
return False
|
||
else:
|
||
for var in variables:
|
||
if var not in keys:
|
||
logger.error("缺少气象数据：", inform[var].long_name)
|
||
return False
|
||
logger.info("variables check true !")
|
||
|
||
all_vars_units = [] # 检查气象数据单位是否正确
|
||
for key in keys:
|
||
all_vars_units.append(nc_data_obj.variables[key].units)
|
||
if key == "level":
|
||
if nc_data_obj.variables[key].units != "millibars":
|
||
logger.error("units errors", inform[key].long_name)
|
||
if key == "t":
|
||
if nc_data_obj.variables[key].units != "K":
|
||
logger.error("units errors", inform[key].long_name)
|
||
if key == "r":
|
||
if nc_data_obj.variables[key].units != "%":
|
||
logger.error("units errors", inform[key].long_name)
|
||
if key == "z":
|
||
if nc_data_obj.variables[key].units != "m**2 s**-2":
|
||
logger.error("units errors", inform[key].long_name)
|
||
logger.info("The units of variables check true !")
|
||
|
||
all_vars_shape = [] # 检查气象数据维数是否正确
|
||
for key in keys:
|
||
all_vars_shape.append(inform[key].shape)
|
||
|
||
if inform["level"].shape[0] != 37: # 检查level是否为37压力面
|
||
logger.error("shape errors", inform[key].long_name)
|
||
if inform["time"].shape[0] != 1: # 检查time是否只有一个时间点
|
||
logger.error("shape errors", inform[key].long_name)
|
||
if (inform["t"].shape[0]) != 1 or (inform["t"].shape[1]) != 37:
|
||
logger.error("shape errors", inform[key].long_name)
|
||
if inform["r"].shape[0] != 1 or inform["r"].shape[1] != 37:
|
||
logger.error("shape errors", inform[key].long_name)
|
||
if inform["z"].shape[0] != 1 or inform["z"].shape[1] != 37:
|
||
logger.error("shape errors", inform[key].long_name)
|
||
logger.info("The shape of variables check true !")
|
||
|
||
@staticmethod
|
||
def read_nc(nc_file):
|
||
"""
|
||
function: 从气象数据中读取数组
|
||
|
||
"""
|
||
nc_data_obj = Nc.Dataset(nc_file)
|
||
inform = nc_data_obj.variables
|
||
|
||
lon_array = inform['longitude'][:] # lon数组
|
||
lat_array = inform['latitude'][:] # lat数组
|
||
n_lat = len(lat_array) # 经线列数量
|
||
n_lon = len(lon_array) # 纬线行数量
|
||
lon_min, lat_max, lon_max, lat_min = [lon_array.min(), lat_array.max(), lon_array.max(), lat_array.min()]
|
||
lon_res = (lon_max - lon_min) / (float(n_lon) - 1) # 计算经度方向分辨率
|
||
lat_res = (lat_max - lat_min) / (float(n_lat) - 1) # 计算纬度方向分辨率
|
||
|
||
pre_level = np.array(inform["level"]) # pressure level in m 压力面
|
||
|
||
temp_array = np.array(inform["t"]) # temp--k (1,37,25,25)
|
||
rel_hum_array = np.array(inform["r"]) # relative_humidity--% (1,37,25,25)
|
||
geo_h_array = np.array(inform["z"]) # Geopotential--m**2 s**-2 (1,37,25,25)
|
||
|
||
return lon_array, lat_array, pre_level, temp_array, rel_hum_array, geo_h_array, n_lat, n_lon, lon_res, lat_res
|
||
|
||
def redu_dim(self, nc_file):
|
||
"""
|
||
function : 降低数组维度
|
||
output:
|
||
lon_array: (25)-->(1,25,37)
|
||
lat_array: (25)-->(1,25,37)
|
||
pre_level: (37)-->(25,25,37)
|
||
temp_array: (1,37,25,25)-->(25,25,37)
|
||
rel_hum_array: (1,37,25,25)-->(25,25,37)
|
||
geo_h_array: (1,37,25,25)-->(25,25,37)
|
||
|
||
"""
|
||
lon_array, lat_array, pre_level, temp_array, rel_hum_array, geo_h_array, n_lat,\
|
||
n_lon, lon_res, lat_res = self.read_nc(nc_file)
|
||
|
||
lat_01 = np.zeros((1, n_lat, 37), dtype=float)
|
||
lon_01 = np.zeros((1, n_lon, 37), dtype=float)
|
||
level_01 = np.zeros((n_lat, n_lon, 37), dtype=float)
|
||
temp_01 = np.zeros((n_lat, n_lon, 37), dtype=float)
|
||
rel_hum_01 = np.zeros((n_lat, n_lon, 37), dtype=float)
|
||
geo_h_01 = np.zeros((n_lat, n_lon, 37), dtype=float)
|
||
|
||
for i in range(0, 37):
|
||
lon_01[0, :, i] = lon_array
|
||
lat_01[0, :, i] = lat_array
|
||
# pre_level = pre_level/100
|
||
level_01[:, :, i] = level_01[:, :, i] + pre_level[i]
|
||
temp_01[0:n_lat, 0:n_lon, i] = temp_array[0, i, :, :]
|
||
rel_hum_01[0:n_lat, 0:n_lon, i] = rel_hum_array[0, i, :, :]
|
||
geo_h_01[:, :, i] = geo_h_array[0, i, :, :]
|
||
|
||
del lon_array, lat_array, pre_level, temp_array, rel_hum_array, geo_h_array
|
||
lon_array, lat_array, pre_level, temp_array, rel_hum_array, geo_h_array =\
|
||
lon_01, lat_01, level_01, temp_01, rel_hum_01, geo_h_01
|
||
return lon_array, lat_array, pre_level, temp_array, rel_hum_array, geo_h_array, lon_res, lat_res
|
||
|
||
def tran_tif(self, nc_file, out_path, str_name):
|
||
"""
|
||
function: 气象数据矩阵转成tiff数据
|
||
param: nc_infile 气象数据文件
|
||
param: out_path 输出的文件夹路径
|
||
param: str 气象数据名字
|
||
"""
|
||
|
||
lon_array, lat_array, pre_level, temp_array, rel_hum_array, geo_h_array,\
|
||
lon_res, lat_res = self.redu_dim(nc_file)
|
||
|
||
# im_bands, n_lat, n_lon = temp_array.shape
|
||
dst_srs = osr.SpatialReference()
|
||
dst_srs.ImportFromEPSG(4326) # 定义一个WGS84坐标系
|
||
im_proj = dst_srs.ExportToWkt()
|
||
im_geo = [lon_array[0, 0, 0], lon_res, 0, lat_array[0, 0, 0], 0, -lat_res]
|
||
# first_name = (nc_file.split("\\", -1)[-1]).split(".nc")[0]
|
||
|
||
# out_level_path = out_path+str_name+"_"+"level.tif"
|
||
# ImageHandler.write_era_into_img(out_level_path, im_proj, im_geo, pre_level)
|
||
out_temp_path = out_path+str_name+"_"+"temp.tif"
|
||
ImageHandler.write_era_into_img(out_temp_path, im_proj, im_geo, temp_array)
|
||
out_hum_path = out_path+str_name+"_"+"re_hum.tif"
|
||
ImageHandler.write_era_into_img(out_hum_path, im_proj, im_geo, rel_hum_array)
|
||
out_geo_h_path = out_path+str_name+"_"+"geo_h.tif"
|
||
ImageHandler.write_era_into_img(out_geo_h_path, im_proj, im_geo, geo_h_array)
|
||
|
||
def copy_tif(self, nc_tif_file, out_path, key_name):
|
||
"""
|
||
tif格式的气象数据：1.识别文件夹中的气象参数；2.复制气象参数至临时文件夹
|
||
param: nc_tif_file:存放气象数据的文件夹路径
|
||
param: out_path：气象数据的复制路径
|
||
param: key_name：MasterNC/AuxiliaryNC
|
||
"""
|
||
filename = os.listdir(nc_tif_file)
|
||
for name in filename:
|
||
if "geo" in name:
|
||
geo_path = os.path.join(nc_tif_file, name)
|
||
out_path1 = os.path.join(out_path, key_name+"_geo_h.tif")
|
||
shutil.copy(geo_path, out_path1)
|
||
elif "hum" in name:
|
||
hum_path = os.path.join(nc_tif_file, name)
|
||
out_path2 = os.path.join(out_path, key_name+"_re_hum.tif")
|
||
shutil.copy(hum_path, out_path2)
|
||
elif "temp" in name:
|
||
temp_path = os.path.join(nc_tif_file, name)
|
||
out_path3 = os.path.join(out_path, key_name+"_temp.tif")
|
||
shutil.copy(temp_path, out_path3)
|
||
|
||
|
||
class ReadImage:
|
||
"""
|
||
读取dem影像, 并输出dem矩阵
|
||
读取主辅影像，并计算干涉相位图
|
||
"""
|
||
def __int__(self):
|
||
pass
|
||
|
||
@staticmethod
|
||
def read_dem(filename):
|
||
"""
|
||
function: 从input目录下读取dem,返回DEM对应的数值数组
|
||
param: fileName-- dem文件路径
|
||
"""
|
||
dem = ImageHandler().get_band_array(filename, 1)
|
||
# dem = np.array(Image.open(filename))
|
||
max_dem = math.ceil(np.max(dem) / 100) * 100 + 50 # 四舍五入向上取值
|
||
return dem, max_dem
|
||
|
||
@staticmethod
|
||
def read_lon_lat(filename):
|
||
"""
|
||
读取影像的经度，纬度
|
||
"""
|
||
lon = []
|
||
lat = []
|
||
im_scope = ImageHandler().get_scope(filename)
|
||
width, height = ImageHandler().get_img_width(filename), ImageHandler().get_img_height(filename) # 获取数据的宽度、高度
|
||
point_upleft, point_downleft = im_scope[1], im_scope[3] # 获取左上数据的[经度、纬度]，右下数据的[经度、纬度]
|
||
geotran = ImageHandler().get_geotransform(filename)
|
||
pix_width, pix_hight = geotran[1], geotran[5] # 像元的宽度，像元的高度 # 获取像元的宽度和像元的高度
|
||
for i in range(0, width):
|
||
longitude = point_upleft[0] + i * pix_width
|
||
lon.append(longitude)
|
||
for j in range(0, height):
|
||
latitude = point_upleft[1] + j * pix_hight
|
||
lat.append(latitude)
|
||
lon = np.asarray(lon)
|
||
lat = np.asarray(lat)
|
||
lon_min, lat_max, lon_max, lat_min = [lon.min(), lat.max(), lon.max(), lat.min()]
|
||
n_lat = len(lat) # 经线列数量
|
||
n_lon = len(lon) # 纬线行数量
|
||
nncols = n_lat-1
|
||
nnrows = n_lon-1
|
||
|
||
lat_37 = np.zeros((1, n_lat, 37), dtype=float)
|
||
lon_37 = np.zeros((1, n_lon, 37), dtype=float)
|
||
for i in range(0, 37):
|
||
d = lon
|
||
lon_37[0, :, i] = d
|
||
|
||
e = lat
|
||
lat_37[0, :, i] = e
|
||
return lon, lat, lon_min, lat_max, lon_max, lat_min, n_lat, n_lon,\
|
||
nncols, nnrows, pix_width, pix_hight, lat_37, lon_37
|
||
|
||
@staticmethod
|
||
def read_phase(file1, file2, out_path):
|
||
"""
|
||
function： 计算干涉相位
|
||
"""
|
||
f01_dataset = gdal.Open(file1) # 获取主影像波段
|
||
f01_width = f01_dataset.RasterXSize
|
||
f01_height = f01_dataset.RasterYSize
|
||
f01_bands = f01_dataset.RasterCount
|
||
img_pro = f01_dataset.GetProjection()
|
||
img_tra = f01_dataset.GetGeoTransform()
|
||
f01_bands1 = f01_dataset.GetRasterBand(1)
|
||
f01_bands2 = f01_dataset.GetRasterBand(2)
|
||
m_r_1 = f01_bands1.ReadAsArray(0, 0, f01_width, f01_height)
|
||
m_i_1 = f01_bands2.ReadAsArray(0, 0, f01_width, f01_height)
|
||
|
||
f02_dataset = gdal.Open(file2) # 获取辅影像波段
|
||
f02_width = f02_dataset.RasterXSize
|
||
f02_height = f02_dataset.RasterYSize
|
||
f02_bands = f02_dataset.RasterCount
|
||
f02_bands1 = f02_dataset.GetRasterBand(1)
|
||
f02_bands2 = f02_dataset.GetRasterBand(2)
|
||
a_r_1 = f02_bands1.ReadAsArray(0, 0, f02_width, f02_height)
|
||
a_i_1 = f02_bands2.ReadAsArray(0, 0, f02_width, f02_height)
|
||
|
||
real_int_ms = m_r_1 * a_r_1 + m_i_1 * a_i_1
|
||
imag_int_ms = m_i_1 * a_r_1 - m_r_1 * a_i_1
|
||
|
||
a = real_int_ms
|
||
b = imag_int_ms
|
||
|
||
fan1 = np.where(a <= 0)
|
||
a = np.where(a == 0, -9999, a)
|
||
phase1 = np.arctan(b / a)
|
||
phase1[fan1] = 0
|
||
del fan1
|
||
gc.collect()
|
||
|
||
fan2 = np.where(a == 0, 1, 0)
|
||
fan3 = np.where(b > 0, 1, 0)
|
||
phase2 = fan2 * fan3 * math.pi / 2
|
||
del fan2, fan3
|
||
gc.collect()
|
||
|
||
fan4 = np.where(a == 0, 1, 0)
|
||
fan5 = np.where(b < 0, 1, 0)
|
||
phase3 = fan4 * fan5 * (-math.pi / 2)
|
||
del fan4, fan5
|
||
gc.collect()
|
||
|
||
fan6 = np.where(a < 0, 1, 0)
|
||
fan7 = np.where(b >= 0, 1, 0)
|
||
a = np.where(a == 0, -9999, a)
|
||
phase4 = fan6 * fan7 * (np.arctan(b / a) + math.pi)
|
||
del fan6, fan7
|
||
gc.collect()
|
||
|
||
fan8 = np.where(a < 0, 1, 0)
|
||
fan9 = np.where(b < 0, 1, 0)
|
||
a = np.where(a == 0, -9999, a)
|
||
phase5 = fan8 * fan9 * (np.arctan(b / a) - math.pi)
|
||
del fan8, fan9
|
||
gc.collect()
|
||
|
||
phase = phase1 + phase2 + phase3 + phase4 + phase5
|
||
logger.info("Calculate phase complete !")
|
||
ImageHandler.write_era_into_img(out_path, img_pro, img_tra, phase)
|
||
|