609 lines
26 KiB
Python
609 lines
26 KiB
Python
"""
|
||
@Project ：microproduct
|
||
@File ：AtmosphericDelayAlg.py
|
||
@Function ：主程序，执行入口
|
||
@Author ：LMM
|
||
@Date ：2021/09/3 14:39
|
||
@Version ：1.0.0
|
||
"""
|
||
|
||
import numpy as np
|
||
import netCDF4 as Nc
|
||
from osgeo import gdal
|
||
import math
|
||
from PIL import Image
|
||
import scipy.spatial
|
||
import scipy.spatial.transform._rotation_groups
|
||
# from scipy.interpolate import interp1d
|
||
# from scipy import interpolate, integrate
|
||
from scipy import interpolate, integrate
|
||
# from pykrige.ok import OrdinaryKriging
|
||
# import pykrige
|
||
from pykrige import OrdinaryKriging
|
||
from xml.etree.ElementTree import ElementTree
|
||
from tool.algorithm.image.ImageHandle import ImageHandler
|
||
from AtmosphericDelayAuxData import ReadImage
|
||
import logging
|
||
from tool.algorithm.algtools.logHandler import LogHandler
|
||
import gc
|
||
LogHandler.init_log_handler(r"run_log\AtmosphericDelay")
|
||
logger = logging.getLogger("mylog")
|
||
|
||
Rd = 287.05 # [j/kg/K] specific gas constants dry air
|
||
Rv = 461.495 # [j/kg/K] water vapour
|
||
k1 = 0.776*100 # [K/hPa]
|
||
k2 = 0.716*100 # [K/hPa]
|
||
k22 = 0.233*100 # [K/hPa]
|
||
k3 = 3.75*1000*100 # [K^2/hPa]
|
||
zref = 15000 # zref for integration- tropopause
|
||
zincr = 15 # z increment spacing for vertical interpolation before integration
|
||
vertres = 5 # vertical resolution for comparing dem to vert profiles of delay
|
||
|
||
class AtmosphericDelay:
|
||
"""
|
||
function： 计算天顶距延迟值ZTD
|
||
"""
|
||
def __init__(self):
|
||
pass
|
||
|
||
def aps_wrf_sar(self, file, temp, re_hum, geo_height, dem_file, p_level=None):
|
||
"""
|
||
计算干湿图
|
||
para: file：NC_geo_h.tif :参考数据，用来获取气象数据的经纬度、分辨率等信息
|
||
para: temp：NC_temp.tif
|
||
para: re_hum：NC_re_hum.tif
|
||
para: geo_height：NC_geo_h.tif
|
||
para: p_level：NC_level.tif
|
||
para: dem_file：dem影像数据
|
||
"""
|
||
dem, max_dem = ReadImage().read_dem(dem_file)
|
||
lon, lat, lon_min, lat_max, lon_max, lat_min, n_lat, n_lon, nncols, nnrows, lon_res, lat_res,\
|
||
lats, lons = ReadImage().read_lon_lat(file)
|
||
lon_dem, lat_dem, lon_min_dem, lat_max_dem, lon_max_dem, lat_min_dem, n_lat_dem, n_lon_dem, nncols_dem, nnrows_dem, lon_res_dem, lat_res_dem, \
|
||
lats_dem, lons_dem = ReadImage().read_lon_lat(dem_file)
|
||
min_lon_lat_dem = min(n_lon_dem, n_lat_dem)
|
||
|
||
logger.info('ERA data ready !')
|
||
# 1、 计算e
|
||
svpw = 6.1121 * np.exp(((17.502 * (temp - 273.16)) / (240.97 + temp - 273.16)))
|
||
svpi = 6.1121 * np.exp((22.587 * (temp - 273.16)) / (273.86 + temp - 273.16))
|
||
temp_bound1 = 273.16 # 0
|
||
temp_bound2 = 250.16 # -23
|
||
wgt = (temp - temp_bound2)/(temp_bound1 - temp_bound2)
|
||
svp = svpi + (svpw - svpi) * (wgt ** 2)
|
||
for k in range(0, 37):
|
||
for i in range(0, n_lat):
|
||
for j in range(0, n_lon):
|
||
if temp[i, j, k].all() > 273.16:
|
||
svp[i, j, k] = svpw[i, j, k]
|
||
elif temp[i, j, k].all() < 250.16:
|
||
svp[i, j, k] = svpi[i, j, k]
|
||
else:
|
||
svp[i, j, k] = svp[i, j, k]
|
||
|
||
# svp = 6.11 * np.exp(2500000 / Rv * ((temp - 273.16) / (temp * 273.16)))
|
||
e = re_hum*svp/100
|
||
# 2、 计算pressure
|
||
pressure = p_level # {ndarray:(5,5,37)}
|
||
# 3、 计算get_z
|
||
g0 = 9.80665
|
||
# 重力加速度随纬度和高度而变化
|
||
# geo_height = geo_height / g0
|
||
# g = 9.80616 * (1 - 0.002637 * np.cos(2 * lats*(math.pi/180)) + 0.0000059 * (np.cos(2 * lats*(math.pi/180)))**2)
|
||
# r_max = 6378137
|
||
# r_min = 6356752
|
||
# r_e = np.sqrt(1/(((np.cos(lats*(math.pi/180))**2)/(r_max**2)) + ((np.sin(lats*(math.pi/180))**2)/(r_min**2))))
|
||
# get_z = (geo_height * r_e) / (g / g0 * r_e - geo_height)
|
||
get_z = geo_height/g0
|
||
|
||
cd_slices = max_dem // vertres + 1
|
||
min_lon_lat = min(n_lon, n_lat)
|
||
max_lon_lat = max(n_lon, n_lat)
|
||
cd_stack_dry = np.zeros((min_lon_lat, cd_slices), dtype=float)
|
||
cd_stack_wet = np.zeros((min_lon_lat, cd_slices), dtype=float)
|
||
# 4、 计算0-15000插值
|
||
x_i = np.array(range(0, zref, zincr))
|
||
logger.info('begin calc wet and dry of diff levels !')
|
||
for m in range(0, min_lon_lat):
|
||
yn = m
|
||
xn = m
|
||
geopo_hgt = get_z[yn, xn, :] # 位势高度 [37,]
|
||
# y_p = pressure[yn, xn, :] # 压力水平level [37,]
|
||
y_p = np.array([1, 2, 3, 5, 7, 10, 20, 30, 50, 70, 100, 125, 150, 175, 200, 225, 250, 300, 350, 400,
|
||
450, 500, 550, 600, 650, 700, 750, 775, 800, 825, 850, 875, 900, 925, 950, 975, 1000])
|
||
y_e = e[yn, xn, :] # 相对湿度 [37,]
|
||
y_t = temp[yn, xn, :] # 温度 [37,]
|
||
|
||
y_e_i = np.zeros((1, 1000), dtype=float)
|
||
y_p_i = np.zeros((1, 1000), dtype=float)
|
||
y_t_i = np.zeros((1, 1000), dtype=float)
|
||
# 判断位势能是否在0-15000之间，避免超出范围的位势能出现插值异常
|
||
if np.min(x_i) > np.min(geopo_hgt) and np.max(x_i) < np.max(geopo_hgt):
|
||
tre_yei = interpolate.interp1d(geopo_hgt, y_e)
|
||
y_e_i = tre_yei(x_i)
|
||
tre_yt = interpolate.interp1d(geopo_hgt, y_t)
|
||
y_t_i = tre_yt(x_i)
|
||
tre_yp = interpolate.interp1d(geopo_hgt, y_p)
|
||
y_p_i = tre_yp(x_i)
|
||
if np.min(x_i) < np.min(geopo_hgt) and np.max(x_i) < np.max(geopo_hgt):
|
||
for f in range(0, 999):
|
||
if x_i[f+1] > np.min(geopo_hgt) > x_i[f]:
|
||
tre_yei = interpolate.interp1d(geopo_hgt, y_e)
|
||
y_e_i[0, f+1:] = tre_yei(x_i[f+1:])
|
||
y_e_i[0, 0:f + 1] = y_e[np.argmin(geopo_hgt)]
|
||
|
||
tre_yt = interpolate.interp1d(geopo_hgt, y_t)
|
||
y_t_i[0, f+1:] = tre_yt(x_i[f+1:])
|
||
y_t_i[0, 0:f + 1] = y_t[np.argmin(geopo_hgt)]
|
||
|
||
tre_yp = interpolate.interp1d(geopo_hgt, y_p)
|
||
y_p_i[0, f+1:] = tre_yp(x_i[f+1:])
|
||
y_p_i[0, 0:f + 1] = y_p[np.argmin(geopo_hgt)] # {ndarray:(1:1000)}
|
||
break
|
||
if np.min(x_i) > np.min(geopo_hgt) and np.max(x_i) > np.max(geopo_hgt):
|
||
for g in range(0, 999):
|
||
if x_i[g + 1] > np.max(geopo_hgt) > x_i[g]:
|
||
tre_yei = interpolate.interp1d(geopo_hgt, y_e)
|
||
y_e_i[0, 0:g+1] = tre_yei(x_i[0:g+1])
|
||
y_e_i[0, g+1:] = y_e[np.argmax(geopo_hgt)]
|
||
|
||
tre_yt = interpolate.interp1d(geopo_hgt, y_t)
|
||
y_t_i[0, 0:g+1] = tre_yt(x_i[0:g+1])
|
||
y_t_i[0, g+1:] = y_t[np.argmax(geopo_hgt)]
|
||
|
||
tre_yp = interpolate.interp1d(geopo_hgt, y_p)
|
||
y_p_i[0, 0:g+1] = tre_yp(x_i[0:g+1])
|
||
y_p_i[0, g+1:] = y_p[np.argmax(geopo_hgt)]
|
||
if np.min(x_i) < np.min(geopo_hgt) and np.max(x_i) > np.max(geopo_hgt):
|
||
for g in range(0, 999):
|
||
if x_i[g + 1] > np.min(geopo_hgt) > x_i[g]:
|
||
m = g+1
|
||
y_e_i[0, 0:m] = y_e[np.argmin(geopo_hgt)]
|
||
y_t_i[0, 0:m] = y_t[np.argmin(geopo_hgt)]
|
||
y_p_i[0, 0:m] = y_p[np.argmin(geopo_hgt)]
|
||
if x_i[g + 1] > np.max(geopo_hgt) > x_i[g]:
|
||
n = g+1
|
||
y_e_i[0, n:] = y_e[np.argmax(geopo_hgt)]
|
||
y_t_i[0, n:] = y_t[np.argmax(geopo_hgt)]
|
||
y_p_i[0, n:] = y_p[np.argmax(geopo_hgt)]
|
||
|
||
tre_yei = interpolate.interp1d(geopo_hgt, y_e)
|
||
y_e_i[0, m:n] = tre_yei(x_i[m:n])
|
||
|
||
tre_yt = interpolate.interp1d(geopo_hgt, y_t)
|
||
y_t_i[0, m:n] = tre_yt(x_i[m:n])
|
||
|
||
tre_yp = interpolate.interp1d(geopo_hgt, y_p)
|
||
y_p_i[0, m:n] = tre_yp(x_i[m:n])
|
||
|
||
# 1、 计算湿延迟
|
||
# tmp_01 = ((k2 - (Rd * k1 / Rv)) * YeI / YTI + k3 * YeI / (YTI ** 2))
|
||
tmp_01 = k22 * y_e_i / y_t_i+k3 * y_e_i / (y_t_i ** 2)
|
||
lw = np.zeros((1, 1000), dtype=float)
|
||
lw1 = (10 ** (-6)) * integrate.cumtrapz(tmp_01, x_i) # (1,999)
|
||
lw1 = np.fliplr(lw1) # 左右翻转(1,999)
|
||
lw[0, 0:999] = lw1[0, :]
|
||
|
||
# 2、 计算干延迟
|
||
tmp_02 = k1 * y_p_i / y_t_i
|
||
|
||
ld = np.zeros((1, 1000), dtype=float)
|
||
ld_1 = (10 ** (-6)) * integrate.cumtrapz(tmp_02, x_i) # （1,999）
|
||
ld_1 = np.fliplr(ld_1) # （1,999）
|
||
ld[0, 0:999] = ld_1[0, :]
|
||
|
||
# ld = 2.2768 * YPI/(1-0.00266*np.cos(lats[0, m, 0]*(math.pi/180))-0.00028*XI) # ld计算方法2 横着（1000，）
|
||
# ld =(10 ** (-6)) *(k1*Rd/g0)*(YPI-YPI[0,999]) # ld计算方法3
|
||
# 3、在(0,maxdem)上做插值
|
||
cdi = np.array(range(0, max_dem, vertres))
|
||
xi_range = x_i
|
||
|
||
ld_range = ld[0, :]
|
||
lw_range = lw[0, :]
|
||
# ld_i = interpolate.interp1d(xi_range, ld_range, kind="cubic")
|
||
ld_i = interpolate.interp1d(xi_range, ld_range)
|
||
# lw_i = interpolate.interp1d(xi_range, lw_range, kind="cubic")
|
||
lw_i = interpolate.interp1d(xi_range, lw_range)
|
||
ldi = np.zeros((max_dem//vertres+1, 1), dtype=float)
|
||
lwi = np.zeros((max_dem//vertres+1, 1), dtype=float)
|
||
for i in range(0, max_dem // vertres):
|
||
ldi[i, 0] = ld_i(cdi[i])
|
||
lwi[i, 0] = lw_i(cdi[i])
|
||
cd_stack_dry[m, :] = ldi[:, 0].T
|
||
cd_stack_wet[m, :] = lwi[:, 0].T
|
||
logger.info('finish calc wet and dry delay of diff levels !')
|
||
|
||
cdstack_interp_dry_wet = np.zeros((min_lon_lat_dem, min_lon_lat_dem, cd_stack_dry.shape[1]), dtype=float) # 修改插值网格分辨率
|
||
# cdstack_interp_dry_wet = np.zeros((min_lon_lat, min_lon_lat, cd_stack_dry.shape[1]), dtype=float)
|
||
|
||
# np.zeros((max_lon_lat, 1), dtype=float)
|
||
# 创建data01、data02
|
||
data01 = np.zeros((min_lon_lat, 1), dtype=float)
|
||
data02 = np.zeros((min_lon_lat, 1), dtype=float)
|
||
new_lon = np.zeros((n_lon, 1), dtype=float)
|
||
new_lat = np.zeros((n_lat, 1), dtype=float)
|
||
new_lon[:, 0] = lon
|
||
new_lat[:, 0] = lat
|
||
data01[0:min_lon_lat, 0] = new_lon[0:min_lon_lat, 0]
|
||
data02[0:min_lon_lat, 0] = new_lat[0:min_lon_lat, 0]
|
||
data01 = data01.reshape(data01.shape[0], 1) # 水平转垂直
|
||
data02 = data02.reshape(data02.shape[0], 1) # 水平转垂直
|
||
|
||
lat_max = lat_min+(min_lon_lat-1) * (-lat_res)
|
||
lon_max = lon_min+(min_lon_lat-1) * lon_res
|
||
logger.info('start interpolation dry and wet delay...')
|
||
# 在网格上进行克里金插值
|
||
# 1. 先插值到 气象数据给的分辨率
|
||
try:
|
||
for i in range(0, cd_stack_dry.shape[1] - 1):
|
||
data03 = cd_stack_dry[:, i] + cd_stack_wet[:, i]
|
||
cdstack_interp_dry_wet[:, :, i] = self.kringing(data01, data02, data03, lon_min, lat_max, lon_max,
|
||
lat_min, lon_res_dem, -lat_res_dem, min_lon_lat_dem)
|
||
except:
|
||
logger.info('calc kringing error !')
|
||
logger.info('start calc rbf ... !')
|
||
for i in range(0, cd_stack_dry.shape[1] - 1):
|
||
data03 = cd_stack_dry[:, i] + cd_stack_wet[:, i]
|
||
cdstack_interp_dry_wet[:, :, i] = self.calc_rbf(data01, data02, data03, lon_min, lat_max, lon_max,
|
||
lat_min, lon_res_dem, -lat_res_dem, min_lon_lat_dem)
|
||
logger.info('calc kringing interv finish !')
|
||
return cdstack_interp_dry_wet
|
||
|
||
# try:
|
||
# for i in range(0, cd_stack_dry.shape[1]-1):
|
||
# data03 = cd_stack_dry[:, i]+cd_stack_wet[:, i]
|
||
# cdstack_interp_dry_wet[:, :, i] = self.kringing(data01, data02, data03, lon_min, lat_max, lon_max,
|
||
# lat_min, lon_res, -lat_res, min_lon_lat)
|
||
# except:
|
||
# logger.info('calc kringing error !')
|
||
# logger.info('start calc rbf ... !')
|
||
# for i in range(0, cd_stack_dry.shape[1]-1):
|
||
# data03 = cd_stack_dry[:, i]+cd_stack_wet[:, i]
|
||
# cdstack_interp_dry_wet[:, :, i] = self.calc_rbf(data01, data02, data03, lon_min, lat_max, lon_max,
|
||
# lat_min, lon_res, -lat_res, min_lon_lat)
|
||
# logger.info('calc kringing interv finish !')
|
||
|
||
|
||
@staticmethod
|
||
def calc_ztd(dry_wet_file_name, dem_file):
|
||
"""
|
||
计算ztd
|
||
wet_file_name:干湿校正文件路径
|
||
dem_file: dem文件路径
|
||
"""
|
||
cdstack_interp_dry_wet = ImageHandler().get_all_band_array(dry_wet_file_name)
|
||
dem = ImageHandler().get_all_band_array(dem_file)
|
||
nncols = ImageHandler().get_img_width(dem_file)
|
||
nnrows = ImageHandler().get_img_height(dem_file)
|
||
|
||
wet_dry_wid = ImageHandler().get_img_width(dry_wet_file_name)
|
||
dem_wid = ImageHandler().get_img_width(dem_file)
|
||
|
||
wet_dry_hgt = ImageHandler().get_img_height(dry_wet_file_name)
|
||
dem_hgt = ImageHandler().get_img_height(dem_file)
|
||
|
||
wetcorrection = np.zeros((nnrows, nncols), dtype=float)
|
||
rounddem = np.round(dem/vertres)
|
||
rounddem[np.where(dem < 0)] = 0
|
||
|
||
ratio_wid = (dem_wid/wet_dry_wid)
|
||
ratio_hgt = (dem_hgt/wet_dry_hgt)
|
||
|
||
for i in range(1, wet_dry_hgt+1):
|
||
for j in range(1, wet_dry_wid+1):
|
||
|
||
m1 = round((i - 1) * ratio_hgt)
|
||
m2 = round(i * ratio_hgt)
|
||
n1 = round((j - 1) * ratio_wid)
|
||
n2 = round(j * ratio_wid)
|
||
|
||
t1_array02 = np.zeros((m2-m1, n2-n1), dtype=float)
|
||
|
||
for h in range(0, cdstack_interp_dry_wet.shape[2]):
|
||
t1_array = rounddem[m1:m2, n1:n2]
|
||
t1_array01 = np.where(t1_array == h, 1, 0)
|
||
t1_array01 = np.where(t1_array01 > 0, cdstack_interp_dry_wet[i-1, j-1, h], 0)
|
||
t1_array02 = t1_array02+t1_array01
|
||
wet_block_array = t1_array02
|
||
wetcorrection[m1:m2, n1:n2] = wet_block_array
|
||
|
||
ztd = wetcorrection
|
||
logger.info('calc ztd finish !')
|
||
return ztd
|
||
|
||
@staticmethod
|
||
def read_nc(infile):
|
||
"""
|
||
读取气象数据
|
||
infile：气象数据路径
|
||
"""
|
||
nc_data_obj = Nc.Dataset(infile)
|
||
keys = nc_data_obj.variables.keys()
|
||
|
||
if 'longitude' in keys:
|
||
lons = nc_data_obj.variables['longitude'][:]
|
||
lats = nc_data_obj.variables['latitude'][:]
|
||
|
||
elif 'lon' in keys:
|
||
lons = nc_data_obj.variables['lon'][:]
|
||
lats = nc_data_obj.variables['lat'][:]
|
||
else:
|
||
logger.info("Specify variable names in the .nc file !")
|
||
temp = np.array(nc_data_obj.variables["t"]) # temp in K 温度
|
||
# temp =(temp-255.6357)/0.00153279
|
||
|
||
hum = np.array(nc_data_obj.variables["r"]) # relative humidity in percent 湿度
|
||
# hum = (hum-41.8782447)/0.0012780
|
||
|
||
geo_h = np.array(nc_data_obj.variables["z"]) # Geopotential Height in m 重力势
|
||
# geo_h =(H-237742.3474)/7.243160
|
||
|
||
p_levs = np.array(nc_data_obj.variables["level"])
|
||
|
||
nc_data_obj.close() # 关闭文件
|
||
return lons, lats, temp, hum, geo_h, p_levs
|
||
|
||
@staticmethod
|
||
def read_dem(demfile):
|
||
"""
|
||
读取dem数据
|
||
demfile：dem数据路径
|
||
"""
|
||
dem = np.array(Image.open(demfile))
|
||
maxdem = math.ceil(np.max(dem)/100)*100+50 # 四舍五入向上取值
|
||
return dem, maxdem
|
||
|
||
@staticmethod
|
||
def copy_channel(num_lat, num_lon, air_temp, rel_hum, h, lons, lats, plevs):
|
||
"""
|
||
改变数据维度
|
||
"""
|
||
temp_1 = np.zeros((num_lat, num_lon, 37), dtype=float)
|
||
hum_1 = np.zeros((num_lat, num_lon, 37), dtype=float)
|
||
h_1 = np.zeros((num_lat, num_lon, 37), dtype=float)
|
||
lat_1 = np.zeros((1, num_lat, 37), dtype=float)
|
||
lon_1 = np.zeros((1, num_lon, 37), dtype=float)
|
||
p_level_1 = np.zeros((num_lat, num_lon, 37), dtype=float)
|
||
|
||
for i in range(0, 37):
|
||
a = air_temp[0, i, :, :]
|
||
temp_1[0:num_lat, 0:num_lon, i] = a
|
||
|
||
b = rel_hum[0, i, :, :]
|
||
hum_1[0:num_lat, 0:num_lon, i] = b
|
||
|
||
c = h[0, i, :, :]
|
||
h_1[:, :, i] = c
|
||
del air_temp, rel_hum, h
|
||
air_temp, rel_hum, h = temp_1, hum_1, h_1
|
||
|
||
for i in range(0, 37):
|
||
d = lons
|
||
lon_1[0, :, i] = d
|
||
|
||
e = lats
|
||
lat_1[0, :, i] = e
|
||
del lons, lats
|
||
lons, lats = lon_1, lat_1
|
||
|
||
# Plevs = Plevs/100
|
||
for i in range(0, 37):
|
||
p_level_1[:, :, i] = p_level_1[:, :, i] + plevs[i]
|
||
# print("输出Plevs1", Plevs1[:, :, i])
|
||
p_level = p_level_1
|
||
return air_temp, rel_hum, h, lons, lats, p_level
|
||
|
||
@staticmethod
|
||
def calc_rbf(data01, data002, data003, lon_min, lat_max, lon_max, lat_min, lon_res, lat_res, min_lon_lat):
|
||
"""
|
||
rbf插值
|
||
"""
|
||
# gridx = np.arange(lon_min, lon_max, lon_res) # 三个参数的意思：范围0.0 - 0.6 ，每隔0.1划分一个网格
|
||
# gridy = np.arange(lat_min, lat_max, lat_res)
|
||
|
||
gridx = np.linspace(lon_min, lon_max, min_lon_lat) # 三个参数的意思：范围(lon_min, lon_max) ,每隔Lon_Res划分一个网格
|
||
gridy = np.linspace(lat_min, lat_max, min_lon_lat)
|
||
|
||
func = interpolate.Rbf(data01, data002, data003, function='multiquadric')
|
||
znew = func(gridx, gridy)
|
||
return znew
|
||
|
||
@staticmethod
|
||
def kringing(data01, data02, data03, lon_min, lat_max, lon_max, lat_min, lon_res, lat_res, min_lon_lat):
|
||
"""
|
||
克里金插值方法
|
||
"""
|
||
# gridx = np.arange(lon_min, lon_max+lon_res, lon_res) # 三个参数的意思：范围(lon_min, lon_max) ,每隔Lon_Res划分一个网格
|
||
# gridy = np.arange(lat_min, lat_max+lat_res, lat_res)
|
||
|
||
# ok3d = ok.OrdinaryKriging(data01, data02, data03, variogram_model="linear") # 模型
|
||
# ok3d = ok.OrdinaryKriging(data01, data02, data03) # 模型
|
||
gridx = np.linspace(lon_min, lon_max, min_lon_lat) # 三个参数的意思：范围(lon_min, lon_max) ,每隔Lon_Res划分一个网格
|
||
gridy = np.linspace(lat_min, lat_max, min_lon_lat)
|
||
|
||
# ok3d = pykrige.OrdinaryKriging(data01, data02, data03) # 模型
|
||
ok3d = OrdinaryKriging(data01, data02, data03)
|
||
|
||
k3d1, ss3d = ok3d.execute("grid", gridx, gridy) # k3d1是结果，给出了每个网格点处对应的值
|
||
# print(np.round(k3d1, 6)) # 小数点后6位
|
||
return k3d1
|
||
|
||
@staticmethod
|
||
def write_dry_wet_tif(ref_image_dem, ref_image_nc, in_array, out_path):
|
||
"""
|
||
存为tif数据
|
||
ref_image：参考影像路径
|
||
in_array：输入数组
|
||
outpath：输出路径
|
||
"""
|
||
dataset_dem = gdal.Open(ref_image_dem) # 获取dem宽高及分辨率
|
||
x_size = dataset_dem.RasterXSize
|
||
y_size = dataset_dem.RasterYSize
|
||
trans_dem = dataset_dem.GetGeoTransform()
|
||
|
||
dataset_nc = gdal.Open(ref_image_nc) # 获取气象数据坐标信息
|
||
trans_nc = dataset_nc.GetGeoTransform()
|
||
|
||
tran_geo = [trans_nc[0], trans_dem[1], 0, trans_nc[3], 0, trans_dem[5]]
|
||
if len(in_array.shape) == 3:
|
||
im_height, im_width, im_bands = in_array.shape
|
||
else:
|
||
im_bands, (im_height, im_width) = 1, in_array.shape
|
||
|
||
pro = dataset_nc.GetProjection()
|
||
|
||
gdal.AllRegister()
|
||
drive = gdal.GetDriverByName("GTiff")
|
||
newfile = drive.Create(out_path, y_size, y_size, im_bands, gdal.GDT_Float32)
|
||
newfile.SetProjection(pro)
|
||
newfile.SetGeoTransform(tran_geo)
|
||
if im_bands == 1:
|
||
newfile.GetRasterBand(1).WriteArray(in_array) # 写入数组数据
|
||
else:
|
||
for i in range(0, im_bands):
|
||
newfile.GetRasterBand(i + 1).WriteArray(in_array[:, :, i])
|
||
# dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
|
||
del newfile
|
||
|
||
@staticmethod
|
||
def write_tif(ref_image, in_array, out_path):
|
||
"""
|
||
存为tif数据
|
||
ref_image：参考影像路径
|
||
in_array：输入数组
|
||
outpath：输出路径
|
||
"""
|
||
dataset = gdal.Open(ref_image)
|
||
x_size = dataset.RasterXSize
|
||
y_size = dataset.RasterYSize
|
||
|
||
if len(in_array.shape) == 3:
|
||
im_height, im_width, im_bands = in_array.shape
|
||
else:
|
||
im_bands, (im_height, im_width) = 1, in_array.shape
|
||
|
||
pro = dataset.GetProjection()
|
||
geo = dataset.GetGeoTransform()
|
||
gdal.AllRegister()
|
||
drive = gdal.GetDriverByName("GTiff")
|
||
newfile = drive.Create(out_path, x_size, y_size, im_bands, gdal.GDT_Float32)
|
||
newfile.SetProjection(pro)
|
||
newfile.SetGeoTransform(geo)
|
||
if im_bands == 1:
|
||
newfile.GetRasterBand(1).WriteArray(in_array) # 写入数组数据
|
||
else:
|
||
for i in range(0, im_bands):
|
||
newfile.GetRasterBand(i + 1).WriteArray(in_array[:, :, i])
|
||
# dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
|
||
del newfile
|
||
|
||
@staticmethod
|
||
def write_ztd_tif(ref_image_dem, ref_image_nc, in_array, out_path):
|
||
"""
|
||
存为tif数据
|
||
ref_image：参考影像路径
|
||
in_array：输入数组
|
||
outpath：输出路径
|
||
"""
|
||
dataset_dem = gdal.Open(ref_image_dem)
|
||
x_size = dataset_dem.RasterXSize
|
||
y_size = dataset_dem.RasterYSize
|
||
trans_dem = dataset_dem.GetGeoTransform()
|
||
|
||
dataset_nc = gdal.Open(ref_image_nc) # 获取气象数据坐标信息
|
||
trans_nc = dataset_nc.GetGeoTransform()
|
||
|
||
tran_geo = [trans_nc[0], trans_dem[1], 0, trans_nc[3], 0, trans_dem[5]]
|
||
|
||
if len(in_array.shape) == 3:
|
||
im_height, im_width, im_bands = in_array.shape
|
||
else:
|
||
im_bands, (im_height, im_width) = 1, in_array.shape
|
||
|
||
pro = dataset_nc.GetProjection()
|
||
|
||
gdal.AllRegister()
|
||
drive = gdal.GetDriverByName("GTiff")
|
||
newfile = drive.Create(out_path, x_size, y_size, im_bands, gdal.GDT_Float32)
|
||
newfile.SetProjection(pro)
|
||
newfile.SetGeoTransform(tran_geo)
|
||
if im_bands == 1:
|
||
newfile.GetRasterBand(1).WriteArray(in_array) # 写入数组数据
|
||
else:
|
||
for i in range(0, im_bands):
|
||
newfile.GetRasterBand(i + 1).WriteArray(in_array[:, :, i])
|
||
# dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
|
||
del newfile
|
||
|
||
@staticmethod
|
||
def cal_phase_bias(local_file, ztd_file, waveLength):
|
||
"""
|
||
function: 计算延迟转换,输出像素-DLOS表
|
||
param: local_file : 输入-局部入射角文件路径
|
||
param: ztd_file : 输入-大气延迟表 ZTD路径
|
||
param: waveLength : 影像对应的波长
|
||
"""
|
||
ztd_array = ImageHandler.get_band_array(ztd_file, 1) # 获取大气延迟表ztd
|
||
local_angle_array = ImageHandler.get_band_array(local_file, 1) # 获取局部入射角
|
||
dlos = (1 / np.cos(local_angle_array)) * ztd_array # 计算公式
|
||
phase_bias = 4 * math.pi /waveLength * dlos # 计算主影像相位误差
|
||
return dlos#phase_bias
|
||
|
||
def complex_phase_removal(self, PhaseError_m, PhaseError_a, intervene_sar,waveLength):
|
||
"""
|
||
function: isce输出的复数干涉图 除以e^(0+wj),得到去除延迟相位的复数干涉图
|
||
param: PhaseError_m: 主影像相位误差图
|
||
param: PhaseError_a：辅影像相位误差图
|
||
param: intervene_sar：原始干涉图
|
||
"""
|
||
phase_bias_m = ImageHandler.get_band_array(PhaseError_m, 1) # 计算主影像相位误差
|
||
phase_bias_a = ImageHandler.get_band_array(PhaseError_a, 1) # 计算辅影像相位误差
|
||
|
||
bias =(4 * math.pi /waveLength)*(phase_bias_m-phase_bias_a) # 延迟相位
|
||
sar_phase = ImageHandler.get_band_array(intervene_sar, 1)
|
||
complex64_sar_phase = sar_phase.astype(np.complex64)
|
||
complex_phase = bias * (1j)
|
||
fl64_com = complex_phase.astype(np.complex64) # 另存为complex64
|
||
bias_result = complex64_sar_phase/np.exp(fl64_com)
|
||
bias_result[np.isnan(bias_result)] = 0 + 0j
|
||
return bias_result
|
||
|
||
@staticmethod
|
||
def get_radarwavelength(IW1):
|
||
"""
|
||
读取波长
|
||
:param IW1:
|
||
:return:
|
||
"""
|
||
tree = ElementTree()
|
||
tree.parse(IW1)
|
||
root = tree.getroot()
|
||
RadarCenterFrequency = root.find('component').find('component').find('component')
|
||
element_trees = list(RadarCenterFrequency)
|
||
value = 0
|
||
for element in element_trees:
|
||
if len(element.attrib) == 1:
|
||
if element.attrib["name"]=="radarwavelength":
|
||
value=float(element.find('value').text)
|
||
return value
|
||
|
||
@staticmethod
|
||
def get_polarization(IW1):
|
||
"""
|
||
读取极化方式
|
||
:param IW1:
|
||
:return:
|
||
"""
|
||
tree = ElementTree()
|
||
tree.parse(IW1)
|
||
root = tree.getroot()
|
||
RadarCenterFrequency = root.find('component').find('component').find('component')
|
||
element_trees = list(RadarCenterFrequency)
|
||
value = 0
|
||
for element in element_trees:
|
||
if len(element.attrib) == 1:
|
||
if element.attrib["name"]=="polarization":
|
||
value=element.find('value').text
|
||
return value
|
||
|