KJ135
/
microproduct
mirror of http://172.16.0.12:5000/WBJY/microproduct.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
microproduct/backScattering/BackScatteringAlg.py

502 lines
21 KiB
Python
Raw Normal View History

入所测试定稿上传
2023-08-28 10:17:29 +00:00
# -*- coding: UTF-8 -*-
"""
@Project ：microproduct
@File ：BackscatteringAlg.py
@Author ：KHZ
@Date ：2021/9/2 11:14
@Version ：1.0.0
修改历史：
[修改序列] [修改日期] [修改者] [修改内容]
1 2022-6-29 石海军 1.int16矩阵转float32；2.对数形式（单位dB）转指数形式（无单位）
"""
import os
import numpy as np
from tool.algorithm.algtools.MetaDataHandler import MetaDataHandler
from tool.algorithm.image.ImageHandle import ImageHandler
from osgeo import gdal, gdalconst
from concurrent.futures import ThreadPoolExecutor, as_completed
# env_str = os.getcwd()
# os.environ['PROJ_LIB'] = env_str
class ScatteringAlg:
def __init__(self):
pass
@staticmethod
def sar_backscattering_coef(in_sar_tif, meta_file_path, out_sar_tif, replece_VV = False, is_DB = True):
# 读取原始SAR影像
proj, geotrans, in_data = ImageHandler.read_img(in_sar_tif)
# 计算强度信息
I = np.array(in_data[0], dtype="float32")
Q = np.array(in_data[1], dtype="float32")
where_9999_0 = np.where(I == -9999)
where_9999_1 = np.where(Q == -9999)
I[where_9999_0] = 1.0
Q[where_9999_1] = 1.0
I2 = np.square(I)
Q2 = np.square(Q)
intensity_arr = I2 + Q2
# 获取极化类型
if 'HH' in os.path.basename(in_sar_tif):
polarization = 'HH'
elif 'HV' in os.path.basename(in_sar_tif):
polarization = 'HV'
elif 'VH' in os.path.basename(in_sar_tif):
polarization = 'VH'
elif 'VV' in os.path.basename(in_sar_tif):
polarization = 'VV'
if replece_VV:
polarization = 'HV' #土壤水分算法中可能会用HV替换VV
else:
raise Exception ('there are not HH、HV、VH、VV in path:',in_sar_tif)
# 获取参数
QualifyValue = MetaDataHandler.get_QualifyValue(meta_file_path, polarization)
Kdb = MetaDataHandler.get_Kdb(meta_file_path, polarization)
# 10 * (alog10((b1 * b1 + b2 * b2) * ((3690.385986 / 32767) * (3690.385986 / 32767)))) - 32.565000
# 计算后向散射系数
#对数形式
coef_arr = 10 * (np.log10(intensity_arr * ((QualifyValue/32767)**2))) - Kdb
coef_arr[np.isnan(coef_arr)] = -9999
coef_arr[np.isinf(coef_arr)] = -9999
coef_arr[where_9999_0] = -9999
coef_arr[where_9999_1] = -9999
# 输出的SAR后向散射系数产品
ImageHandler.write_img(out_sar_tif, proj, geotrans, coef_arr, -9999)
return True
###############
# RPC 模块
###############
def parse_rpc_file(rpc_file):
# rpc_file:.rpc文件的绝对路径
# rpc_dict：符号RPC域下的16个关键字的字典
# 参考网址：http://geotiff.maptools.org/rpc_prop.html；
# https://www.osgeo.cn/gdal/development/rfc/rfc22_rpc.html
rpc_dict = {}
with open(rpc_file) as f:
text = f.read()
# .rpc文件中的RPC关键词
words = ['errBias', 'errRand', 'lineOffset', 'sampOffset', 'latOffset',
'longOffset', 'heightOffset', 'lineScale', 'sampScale', 'latScale',
'longScale', 'heightScale', 'lineNumCoef', 'lineDenCoef','sampNumCoef', 'sampDenCoef',]
# GDAL库对应的RPC关键词
keys = ['ERR_BIAS', 'ERR_RAND', 'LINE_OFF', 'SAMP_OFF', 'LAT_OFF', 'LONG_OFF',
'HEIGHT_OFF', 'LINE_SCALE', 'SAMP_SCALE', 'LAT_SCALE',
'LONG_SCALE', 'HEIGHT_SCALE', 'LINE_NUM_COEFF', 'LINE_DEN_COEFF',
'SAMP_NUM_COEFF', 'SAMP_DEN_COEFF']
for old, new in zip(words, keys):
text = text.replace(old, new)
# 以‘;\n’作为分隔符
text_list = text.split(';\n')
# 删掉无用的行
text_list = text_list[3:-2]
#
text_list[0] = text_list[0].split('\n')[1]
# 去掉制表符、换行符、空格
text_list = [item.strip('\t').replace('\n', '').replace(' ', '') for item in text_list]
for item in text_list:
# 去掉‘=’
key, value = item.split('=')
# 去掉多余的括号‘(’，‘)’
if '(' in value:
value = value.replace('(', '').replace(')', '')
rpc_dict[key] = value
for key in keys[:12]:
# 为正数添加符号‘+’
if not rpc_dict[key].startswith('-'):
rpc_dict[key] = '+' + rpc_dict[key]
# 为归一化项和误差标志添加单位
if key in ['LAT_OFF', 'LONG_OFF', 'LAT_SCALE', 'LONG_SCALE']:
rpc_dict[key] = rpc_dict[key] + ' degrees'
if key in ['LINE_OFF', 'SAMP_OFF', 'LINE_SCALE', 'SAMP_SCALE']:
rpc_dict[key] = rpc_dict[key] + ' pixels'
if key in ['ERR_BIAS', 'ERR_RAND', 'HEIGHT_OFF', 'HEIGHT_SCALE']:
rpc_dict[key] = rpc_dict[key] + ' meters'
# 处理有理函数项
for key in keys[-4:]:
values = []
for item in rpc_dict[key].split(','):
#print(item)
if not item.startswith('-'):
values.append('+'+item)
else:
values.append(item)
rpc_dict[key] = ' '.join(values)
return rpc_dict
def write_rpc_to_tiff(inputpath,rpc_file,ap = True,outpath = None):
rpc_dict = parse_rpc_file(rpc_file)
if ap:
# 可修改读取
dataset = gdal.Open(inputpath,gdal.GA_Update)
# 向tif影像写入rpc域信息
# 注意，这里虽然写入了RPC域信息，但实际影像还没有进行实际的RPC校正
# 尽管有些RS/GIS能加载RPC域信息，并进行动态校正
for k in rpc_dict.keys():
dataset.SetMetadataItem(k, rpc_dict[k], 'RPC')
dataset.FlushCache()
del dataset
else:
dataset = gdal.Open(inputpath,gdal.GA_Update)
tiff_driver= gdal.GetDriverByName('Gtiff')
out_ds = tiff_driver.CreateCopy(outpath, dataset, 1)
for k in rpc_dict.keys():
out_ds.SetMetadataItem(k, rpc_dict[k], 'RPC')
out_ds.FlushCache()
def rpc_correction(inputpath,rpc_file,corrtiff,dem_tif_file = None):
## 设置rpc校正的参数
# 原图像和输出影像缺失值设置为0，输出影像坐标系为WGS84(EPSG:4326), 重采样方法为双线性插值（bilinear，还有最邻近‘near’、三次卷积‘cubic’等可选)
# 注意DEM的覆盖范围要比原影像的范围大，此外，DEM不能有缺失值，有缺失值会报错
# 通常DEM在水域是没有值的（即缺失值的情况），因此需要将其填充设置为0，否则在RPC校正时会报错
# 这里使用的DEM是填充0值后的SRTM V4.1 3秒弧度的DEM(分辨率为90m)
# RPC_DEMINTERPOLATION=bilinear 表示对DEM重采样使用双线性插值算法
# 如果要修改输出的坐标系，则要修改dstSRS参数值，使用该坐标系统的EPSG代码
# 可以在网址https://spatialreference.org/ref/epsg/32650/ 查询得到EPSG代码
write_rpc_to_tiff(inputpath,rpc_file,ap = True,outpath = None)
if dem_tif_file is None:
wo = gdal.WarpOptions(srcNodata=0, dstNodata=0, dstSRS='EPSG:4326', resampleAlg='bilinear',
format='Gtiff',rpc=True, warpOptions=["INIT_DEST=NO_DATA"])
wr = gdal.Warp(corrtiff, inputpath, options=wo)
print("RPC_GEOcorr>>>")
else:
wo = gdal.WarpOptions(srcNodata=0, dstNodata=0, dstSRS='EPSG:4326', resampleAlg='bilinear', format='Gtiff',rpc=True, warpOptions=["INIT_DEST=NO_DATA"],
transformerOptions=["RPC_DEM=%s"%(dem_tif_file), "RPC_DEMINTERPOLATION=bilinear"])
wr = gdal.Warp(corrtiff, inputpath, options=wo)
print("RPC_GEOcorr>>>")
## 对于全海域的影像或者不使用DEM校正的话，可以将transformerOptions有关的RPC DEM关键字删掉
## 即将上面gdal.WarpOptions注释掉，将下面的语句取消注释，无DEM时，影像范围的高程默认全为0
del wr
##########################
# 输出RPC 行列号到 经纬度的变换
##########################
#最大迭代次数超过则报错
class MaxLocalizationIterationsError(Exception):
"""
Custom rpcm Exception.
"""
pass
def apply_poly(poly, x, y, z):
"""
Evaluates a 3-variables polynom of degree 3 on a triplet of numbers.
将三次多项式的统一模式构建为一个单独的函数
Args:
poly: list of the 20 coefficients of the 3-variate degree 3 polynom,
ordered following the RPC convention.
x, y, z: triplet of floats. They may be numpy arrays of same length.
Returns:
the value(s) of the polynom on the input point(s).
"""
out = 0
out += poly[0]
out += poly[1]*y + poly[2]*x + poly[3]*z
out += poly[4]*y*x + poly[5]*y*z +poly[6]*x*z
out += poly[7]*y*y + poly[8]*x*x + poly[9]*z*z
out += poly[10]*x*y*z
out += poly[11]*y*y*y
out += poly[12]*y*x*x + poly[13]*y*z*z + poly[14]*y*y*x
out += poly[15]*x*x*x
out += poly[16]*x*z*z + poly[17]*y*y*z + poly[18]*x*x*z
out += poly[19]*z*z*z
return out
def apply_rfm(num, den, x, y, z):
"""
Evaluates a Rational Function Model (rfm), on a triplet of numbers.
执行20个参数的分子和20个参数的除法
Args:
num: list of the 20 coefficients of the numerator
den: list of the 20 coefficients of the denominator
All these coefficients are ordered following the RPC convention.
x, y, z: triplet of floats. They may be numpy arrays of same length.
Returns:
the value(s) of the rfm on the input point(s).
"""
return apply_poly(num, x, y, z) / apply_poly(den, x, y, z)
def rpc_from_geotiff(geotiff_path):
"""
Read the RPC coefficients from a GeoTIFF file and return an RPCModel object.
该函数返回影像的Gdal格式的影像和RPCmodel
Args:
geotiff_path (str): path or url to a GeoTIFF file
Returns:
instance of the rpc_model.RPCModel class
"""
# with rasterio.open(geotiff_path, 'r') as src:
#
dataset = gdal.Open(geotiff_path, gdal.GA_ReadOnly)
rpc_dict = dataset.GetMetadata("RPC")
# 同时返回影像与rpc
return dataset, RPCModel(rpc_dict,'geotiff')
def read_rpc_file(rpc_file):
"""
Read RPC from a RPC_txt file and return a RPCmodel
从TXT中直接单独读取RPC模型
Args:
rpc_file: RPC sidecar file path
Returns:
dictionary read from the RPC file, or an empty dict if fail
"""
rpc = parse_rpc_file(rpc_file)
return RPCModel(rpc)
class RPCModel:
def __init__(self, d, dict_format="geotiff"):
"""
Args:
d (dict): dictionary read from a geotiff file with
rasterio.open('/path/to/file.tiff', 'r').tags(ns='RPC'),
or from the .__dict__ of an RPCModel object.
dict_format (str): format of the dictionary passed in `d`.
Either "geotiff" if read from the tags of a geotiff file,
or "rpcm" if read from the .__dict__ of an RPCModel object.
"""
if dict_format == "geotiff":
self.row_offset = float(d['LINE_OFF'][0:d['LINE_OFF'].rfind(' ')])
self.col_offset = float(d['SAMP_OFF'][0:d['SAMP_OFF'].rfind(' ')])
self.lat_offset = float(d['LAT_OFF'][0:d['LAT_OFF'].rfind(' ')])
self.lon_offset = float(d['LONG_OFF'][0:d['LONG_OFF'].rfind(' ')])
self.alt_offset = float(d['HEIGHT_OFF'][0:d['HEIGHT_OFF'].rfind(' ')])
self.row_scale = float(d['LINE_SCALE'][0:d['LINE_SCALE'].rfind(' ')])
self.col_scale = float(d['SAMP_SCALE'][0:d['SAMP_SCALE'].rfind(' ')])
self.lat_scale = float(d['LAT_SCALE'][0:d['LAT_SCALE'].rfind(' ')])
self.lon_scale = float(d['LONG_SCALE'][0:d['LONG_SCALE'].rfind(' ')])
self.alt_scale = float(d['HEIGHT_SCALE'][0:d['HEIGHT_SCALE'].rfind(' ')])
self.row_num = list(map(float, d['LINE_NUM_COEFF'].split()))
self.row_den = list(map(float, d['LINE_DEN_COEFF'].split()))
self.col_num = list(map(float, d['SAMP_NUM_COEFF'].split()))
self.col_den = list(map(float, d['SAMP_DEN_COEFF'].split()))
if 'LON_NUM_COEFF' in d:
self.lon_num = list(map(float, d['LON_NUM_COEFF'].split()))
self.lon_den = list(map(float, d['LON_DEN_COEFF'].split()))
self.lat_num = list(map(float, d['LAT_NUM_COEFF'].split()))
self.lat_den = list(map(float, d['LAT_DEN_COEFF'].split()))
elif dict_format == "rpcm":
self.__dict__ = d
else:
raise ValueError(
"dict_format '{}' not supported. "
"Should be {{'geotiff','rpcm'}}".format(dict_format)
)
def projection(self, lon, lat, alt):
"""
Convert geographic coordinates of 3D points into image coordinates.
正投影：从地理坐标到图像坐标
Args:
lon (float or list): longitude(s) of the input 3D point(s)
lat (float or list): latitude(s) of the input 3D point(s)
alt (float or list): altitude(s) of the input 3D point(s)
Returns:
float or list: horizontal image coordinate(s) (column index, ie x)
float or list: vertical image coordinate(s) (row index, ie y)
"""
nlon = (np.asarray(lon) - self.lon_offset) / self.lon_scale
nlat = (np.asarray(lat) - self.lat_offset) / self.lat_scale
nalt = (np.asarray(alt) - self.alt_offset) / self.alt_scale
col = apply_rfm(self.col_num, self.col_den, nlat, nlon, nalt)
row = apply_rfm(self.row_num, self.row_den, nlat, nlon, nalt)
col = col * self.col_scale + self.col_offset
row = row * self.row_scale + self.row_offset
return col, row
def localization(self, col, row, alt, return_normalized=False):
"""
Convert image coordinates plus altitude into geographic coordinates.
反投影：从图像坐标到地理坐标
Args:
col (float or list): x image coordinate(s) of the input point(s)
row (float or list): y image coordinate(s) of the input point(s)
alt (float or list): altitude(s) of the input point(s)
Returns:
float or list: longitude(s)
float or list: latitude(s)
"""
ncol = (np.asarray(col) - self.col_offset) / self.col_scale
nrow = (np.asarray(row) - self.row_offset) / self.row_scale
nalt = (np.asarray(alt) - self.alt_offset) / self.alt_scale
if not hasattr(self, 'lat_num'):
lon, lat = self.localization_iterative(ncol, nrow, nalt)
else:
lon = apply_rfm(self.lon_num, self.lon_den, nrow, ncol, nalt)
lat = apply_rfm(self.lat_num, self.lat_den, nrow, ncol, nalt)
if not return_normalized:
lon = lon * self.lon_scale + self.lon_offset
lat = lat * self.lat_scale + self.lat_offset
return lon, lat
def localization_iterative(self, col, row, alt):
"""
Iterative estimation of the localization function (image to ground),
for a list of image points expressed in image coordinates.
逆投影时的迭代函数
Args:
col, row: normalized image coordinates (between -1 and 1)
alt: normalized altitude (between -1 and 1) of the corresponding 3D
point
Returns:
lon, lat: normalized longitude and latitude
Raises:
MaxLocalizationIterationsError: if the while loop exceeds the max
number of iterations, which is set to 100.
"""
# target point: Xf (f for final)
Xf = np.vstack([col, row]).T
# use 3 corners of the lon, lat domain and project them into the image
# to get the first estimation of (lon, lat)
# EPS is 2 for the first iteration, then 0.1.
lon = -col ** 0 # vector of ones
lat = -col ** 0
EPS = 2
x0 = apply_rfm(self.col_num, self.col_den, lat, lon, alt)
y0 = apply_rfm(self.row_num, self.row_den, lat, lon, alt)
x1 = apply_rfm(self.col_num, self.col_den, lat, lon + EPS, alt)
y1 = apply_rfm(self.row_num, self.row_den, lat, lon + EPS, alt)
x2 = apply_rfm(self.col_num, self.col_den, lat + EPS, lon, alt)
y2 = apply_rfm(self.row_num, self.row_den, lat + EPS, lon, alt)
n = 0
while not np.all((x0 - col) ** 2 + (y0 - row) ** 2 < 1e-18):
if n > 100:
raise MaxLocalizationIterationsError("Max localization iterations (100) exceeded")
X0 = np.vstack([x0, y0]).T
X1 = np.vstack([x1, y1]).T
X2 = np.vstack([x2, y2]).T
e1 = X1 - X0
e2 = X2 - X0
u = Xf - X0
# project u on the base (e1, e2): u = a1*e1 + a2*e2
# the exact computation is given by:
# M = np.vstack((e1, e2)).T
# a = np.dot(np.linalg.inv(M), u)
# but I don't know how to vectorize this.
# Assuming that e1 and e2 are orthogonal, a1 is given by
# <u, e1> / <e1, e1>
num = np.sum(np.multiply(u, e1), axis=1)
den = np.sum(np.multiply(e1, e1), axis=1)
a1 = np.divide(num, den).squeeze()
num = np.sum(np.multiply(u, e2), axis=1)
den = np.sum(np.multiply(e2, e2), axis=1)
a2 = np.divide(num, den).squeeze()
# use the coefficients a1, a2 to compute an approximation of the
# point on the gound which in turn will give us the new X0
lon += a1 * EPS
lat += a2 * EPS
# update X0, X1 and X2
EPS = .1
x0 = apply_rfm(self.col_num, self.col_den, lat, lon, alt)
y0 = apply_rfm(self.row_num, self.row_den, lat, lon, alt)
x1 = apply_rfm(self.col_num, self.col_den, lat, lon + EPS, alt)
y1 = apply_rfm(self.row_num, self.row_den, lat, lon + EPS, alt)
x2 = apply_rfm(self.col_num, self.col_den, lat + EPS, lon, alt)
y2 = apply_rfm(self.row_num, self.row_den, lat + EPS, lon, alt)
n += 1
return lon, lat
#############################################
# 校正影像 输出影像转换表
#############################################
"""
[pool.putRequest(req) for req in requests]等同于
　　for req in requests:
　　　　 pool.putRequest(req)
"""
def rpc_row_col(params):
rpc_mdl,r_block,c_block,ati_block,i,block_shape=params
r_block ,c_block = rpc_mdl.localization(r_block.reshape(-1).astype(np.int32),c_block.reshape(-1).astype(np.int32) ,ati_block )
return [i,r_block ,c_block,block_shape]
def getRCImageRC(inputPath,out_rpc_rc_path,rpc_file_path):
rpc_tool = read_rpc_file(rpc_file_path)
#shutil.copy2(inputPath, out_rpc_rc_path)
# 创建重采样输出文件（设置投影及六参数）
input_rpc_sar=gdal.Open(inputPath)
driver = gdal.GetDriverByName('GTiff')
output = driver.Create(out_rpc_rc_path, input_rpc_sar.RasterXSize,input_rpc_sar.RasterYSize, 2,gdal.GDT_Float32)
output.SetGeoTransform(list(input_rpc_sar.GetGeoTransform()))
output.SetProjection(input_rpc_sar.GetProjection())
# 参数说明 输入数据集、输出文件、输入投影、参考投影、重采样方法(最邻近内插\双线性内插\三次卷积等)、回调函数
rpc_rc_img=output
# 计算行列号
rc_height=rpc_rc_img.RasterYSize
rc_width=rpc_rc_img.RasterXSize
with ThreadPoolExecutor(max_workers=8) as t:
plist=[]
for i in range(0,rc_height,100):
c_block=np.ones((100,rc_width)).astype(np.float32)*(np.array(list(range(rc_width))).reshape(1,rc_width))
r_block=np.ones((100,rc_width)).astype(np.float32)*(np.array(list(range(100))).reshape(100,1))
r_block=r_block+i
if not rc_height-i>100:
num=rc_height-i
r_block=r_block[:num,:].astype(np.float32)
c_block=c_block[:num,:].astype(np.float32)
block_shape=r_block.shape
plist.append(t.submit(rpc_row_col,[rpc_tool,r_block.reshape(-1).astype(np.int32),c_block.reshape(-1).astype(np.int32) ,c_block.reshape(-1)*0+0,i,block_shape]))
for future in as_completed(plist):
data = future.result()
i,r_block ,c_block,block_shape=data
rpc_rc_img.GetRasterBand(1).WriteArray(r_block.reshape(block_shape).astype(np.float32),0,i)
rpc_rc_img.GetRasterBand(2).WriteArray(c_block.reshape(block_shape).astype(np.float32),0,i)
del rpc_rc_img,output,input_rpc_sar