# -*- coding: UTF-8 -*- """ @Project ：microproduct @File ：OneOrthoAlg.py @Function ：正射校正算法 @Author ：KHZ @Contact： @Date ：2021/8/14 @Version ：1.0.0 """ # from re import I, T, X, match # from types import DynamicClassAttribute # from numpy.core.einsumfunc import _parse_possible_contraction # from numpy.core.fromnumeric import shape # from numpy.core.shape_base import block # from numpy.lib import row_stack # from numpy.lib.function_base import delete import psutil import os import gc # import sys # import scipy # from scipy.sparse import data # from scipy.sparse.construct import random import xmltodict import numpy as np import json import datetime, time import yaml import math import time #import cv2 #import cv2 as cv from osgeo import gdal, gdalconst # from yaml.events import AliasEvent, NodeEvent from OrthoAuxData import OrthoAuxData # from OrthoImage import ImageHandler from tool.algorithm.image.ImageHandle import ImageHandler from tool.algorithm.algtools.MetaDataHandler import MetaDataHandler # from logHandler import LogHandler from osgeo import osr from scipy import interpolate import scipy.io as scipyio import copy import scipy.sparse as ss import shutil import pandas import uuid from concurrent.futures._base import as_completed, wait from concurrent.futures.thread import ThreadPoolExecutor from multiprocessing import Pool 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) # 计算后向散射系数 #对数形式 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, 0) tif_array = np.power(10.0, coef_arr / 10.0) # dB --> 线性值 后向散射系数 tif_array[np.isnan(tif_array)] = 0 tif_array[np.isinf(tif_array)] = 0 tif_array[where_9999_0] = 0 tif_array[where_9999_1] = 0 ImageHandler.write_img(out_sar_tif, proj, geotrans, tif_array, 0) return True @staticmethod def lin_to_db(lin_path, db_path): proj, geotrans, in_data = ImageHandler.read_img(lin_path) db_arr = 10 * np.log10(in_data) # db_arr[np.isnan(db_arr)] = -9999 # db_arr[np.isinf(db_arr)] = -9999 ImageHandler.write_img(db_path, proj, geotrans, db_arr, -9999) def FindInfomationFromJson(HeaderFile_dom_json, node_path_list): """ 在Json文件中，按照指定路径解析出制定节点 """ result_node = HeaderFile_dom_json for nodename in node_path_list: result_node = result_node[nodename] return result_node def GetVectorNorm(Vecter): """ 得到向量的模 """ Vecter = Vecter.reshape(-1,1) Vecter_Norm_pow = np.matmul(Vecter.T,Vecter) return np.sqrt(Vecter_Norm_pow) def XYZOuterM2(A, B): """ 外积（叉乘）,日后版本换成可以任意维度的外积运算方程 args: A:nx3 B:nx3 """ cnt = A.shape[0] C = np.zeros((cnt, 3)) C[:, 0] = A[:, 1] * B[:, 2] - A[:, 2] * B[:, 1] C[:, 1] = A[:, 2] * B[:, 0] - A[:, 0] * B[:, 2] C[:, 2] = A[:, 0] * B[:, 1] - A[:, 1] * B[:, 0] return C def LinearSampling_raster(source_raster,tar_raster,source_start_row,source_start_col,tar_start_row,tar_start_col): ''' 双线性重采样 agrs: source_raster:原影像 tar_raster:目标影像 source_start_row:原影像的起始行号 source_start_col:原影像的起始列号 tar_start_row:目标影像的起始行号 tar_start_col:目标影像的起始列号 return: tar_raster:目标影像 ''' # 根据原影像的起始行列号与目标影像行列号，计算目标影像行列号的校正值 d_row=tar_start_row-source_start_row # d_col=tar_start_col-source_start_col # source_height=source_raster.shape[0] source_width=source_raster.shape[1] tar_height=tar_raster.shape[0] tar_width=tar_raster.shape[1] source_row=(np.ones((source_width,source_height))*range(source_height)).transpose(1,0) source_col=np.ones((source_height,source_width))*range(source_width) tar_row=(np.ones((tar_width,tar_height))*range(tar_height)).transpose(1,0) tar_col=np.ones((tar_height,tar_width))*range(tar_width) tar_row=tar_row+d_row # 坐标系变换 tar_col=tar_col+d_col # 坐标系变换 # 确定行列号 last_source_row=np.floor(tar_row) last_source_col=np.floor(tar_col) next_source_row=np.ceil(tar_row) next_source_col=np.ceil(tar_col) # 双线性重采样模型示意图， # y1 r1 y2 # # r # # y3 r2 y4 # 计算重采样 r1=source_raster[last_source_row,last_source_col]*(tar_col-last_source_col)+source_raster[last_source_row,next_source_col]*(next_source_col-tar_col) r2=source_raster[next_source_row,last_source_col]*(tar_col-last_source_col)+source_raster[next_source_row,next_source_col]*(next_source_col-tar_col) tar_raster=r1*(tar_row-last_source_row)+r2*(next_source_row-tar_row) tar_raster=tar_raster.reshape(tar_height,tar_width) # 目标影像 return tar_raster.copy() '''-----双线性重采样-----''' """ 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 # / 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,c_block,0) print("\t{}\t".format(i)) return [i,r_block ,c_block,block_shape] def get_RPC_lon_lat(in_rpc_tiff,out_rpc_rc_path): exe_path=r".\baseTool\x64\Release\SIMOrthoProgram.exe" exe_cmd=r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} ".format(exe_path,8,in_rpc_tiff,out_rpc_rc_path) print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") def getRCImageRC2(inputPath,out_rpc_rc_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 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 rc_height-i>100: rpc_rc_img.GetRasterBand(1).WriteArray(r_block.astype(np.float32),0,i) rpc_rc_img.GetRasterBand(2).WriteArray(c_block.astype(np.float32),0,i) else: num=rc_height-i rpc_rc_img.GetRasterBand(1).WriteArray(r_block[:num,:].astype(np.float32),0,i) rpc_rc_img.GetRasterBand(2).WriteArray(c_block[:num,:].astype(np.float32),0,i) del rpc_rc_img,output,input_rpc_sar 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 Pool(8) as t: plist=[] for i in range(0,rc_height,1000): c_block=np.ones((1000,rc_width)).astype(np.float32)*(np.array(list(range(rc_width))).reshape(1,rc_width)) r_block=np.ones((1000,rc_width)).astype(np.float32)*(np.array(list(range(1000))).reshape(1000,1)) r_block=r_block+i if not rc_height-i>1000: 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.apply_async(rpc_row_col,args=([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],))) t.close() t.join() for future in plist: data = future.get() 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 """ RPC 结束 """ class SatelliteOrbit(object): def __init__(self) -> None: super().__init__() self.starttime = 1262275200.0 self.modelName="" def get_starttime(self): ''' 返回卫星轨道时间起算点 ''' return self.starttime def ReconstructionSatelliteOrbit(self, GPSPoints_list): ''' 重建卫星轨道，使用多项式拟合法 args: GPSPoints_list:GPS 卫星轨道点 return： SatelliteOrbitModel 卫星轨道模型 ''' self.SatelliteOrbitModel = None def SatelliteSpaceState(self, time_float): ''' 根据时间戳，返回对应时间的卫星的轨迹状态 args: time_float:时间戳 return： State_list:[time,Xp,Yp,Zp,Vx,Vy,Vz] ''' return None class SatelliteOrbitFitPoly(SatelliteOrbit): ''' 继承于SatelliteOribit类，为拟合多项式实现方法 ''' def __init__(self) -> None: super().__init__() self.modelName="多项式" self.polynum=4 def ReconstructionSatelliteOrbit(self, GPSPoints_list, starttime): if len(GPSPoints_list)==2: self.polynum=1 self.starttime = starttime record_count = len(GPSPoints_list) time_arr = np.zeros((record_count, 1), dtype=np.float64) # 使用np.float64只是为了精度高些；如果32位也能满足需求，请用32位 state_arr = np.zeros((record_count, 6), dtype=np.float64) A_arr = np.zeros((self.polynum+1, 6), dtype=np.float64) # 四次项 X=np.ones((record_count,self.polynum+1),dtype=np.float64) # 记录时间坐标 # 将点记录转换为自变量矩阵、因变量矩阵 for i in range(record_count): GPSPoint = GPSPoints_list[i] time_ = GPSPoint[0] - self.starttime # 为了保证精度，对时间进行缩放 X[i,:]=np.array([1,time_]) state_arr[i, :] = np.array(GPSPoint[1:],dtype=np.float64).reshape(1,6) # 空间坐标 self.model_f=[] for i in range(6): Y = state_arr[:, i].reshape(-1,1) A_arr[:,i]=np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),Y)[:,0] self.A_arr=copy.deepcopy(A_arr.copy()) return self.A_arr elif len(GPSPoints_list) > 6: self.polynum=4 # 多项式的节点数，理论上是超过5个可以起算，这里为了精度选择10个点起算。 # 多项式 XA=Y ==> A=(X'X)^X'Y，其中 A 为待求系数，X为变量，Y为因变量 # 这里使用三次项多项式，共有6组参数。 # 声明自变量，因变量，系数矩阵 self.starttime = starttime record_count = len(GPSPoints_list) time_arr = np.zeros((record_count, 1), dtype=np.float64) # 使用np.float64只是为了精度高些；如果32位也能满足需求，请用32位 state_arr = np.zeros((record_count, 6), dtype=np.float64) A_arr = np.zeros((self.polynum+1, 6), dtype=np.float64) # 四次项 X=np.ones((record_count,self.polynum+1),dtype=np.float64) # 记录时间坐标 # 将点记录转换为自变量矩阵、因变量矩阵 for i in range(record_count): GPSPoint = GPSPoints_list[i] time_ = GPSPoint[0] - self.starttime # 为了保证精度，对时间进行缩放 X[i,:]=np.array([1,time_,time_**2,time_**3,time_**4]) state_arr[i, :] = np.array(GPSPoint[1:],dtype=np.float64).reshape(1,6) # 空间坐标 self.model_f=[] for i in range(6): Y = state_arr[:, i].reshape(-1,1) A_arr[:,i]=np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),Y)[:,0] self.A_arr=copy.deepcopy(A_arr.copy()) ''' 测试误差 from matplotlib import pyplot label_list=['x','y','z','vx','vy','vz'] color_list=['r','g','b','gold','gray','pink'] pyplot.figure() for i in range(6): Y = state_arr[:, i] Y_predict=self.model_f[i](X) pyplot.subplot(int("23{}".format(i+1))) d=Y-Y_predict pyplot.plot(X,d,label=label_list[i],color=color_list[i]) pyplot.title("max:{}".format(np.max(d))) #self.model_f.append(interpolate.interp1d(X,Y,kind='cubic',fill_value='extrapolate')) pyplot.legend() pyplot.show() ''' return self.A_arr else: self.A_arr = None return None def SatelliteSpaceState(self, time_float): ''' 逐像素求解 根据时间戳，返回对应时间的卫星的轨迹状态，会自动计算与起算时间之差 args: time_float:时间戳 return： State_list:[time,Xp,Yp,Zp,Vx,Vy,Vz] ''' if self.model_f is None: return None result_arr=np.zeros((1,7)) time_float = time_float - self.starttime result_arr[0,0]=time_float #time_arr[0, 4] = time_arr[0, 3] * time_float ** 4 time_float=np.array([1,time_float,time_float**2,time_float**3,time_float**4]).reshape(1,5) result_arr=np.matmul(time_float,self.A_arr) return [time_float,result_arr] def getSatelliteSpaceState(self, time_array): ''' 矩阵求解 根据时间戳矩阵，返回对应时刻的卫星空间状态（位置，速度），且会自动计算与起算时间之差 args: time_array:nparray nx1 时间戳 return: SatellitSpaceStateArray:nparray nx6 状态信息 ''' if self.model_f is None: return None # 返回None,表示没有结果 if self.polynum==4: n=time_array.shape[0] result_arr_=np.zeros((n,6),dtype=np.float64) time_float = time_array - self.starttime time_float=time_float.reshape(-1) # nx1 time_arr=np.ones((time_float.shape[0],5)) # nx5 time_arr[:,1]=time_float time_arr[:,2]=time_float**2 time_arr[:,3]=time_float**3 time_arr[:,4]=time_float**4 result_arr_=np.matmul(time_arr,self.A_arr) # nx5 5x6 #time_arr[0, 4] = time_arr[0, 3] * time_float ** 4 #result_arr=result_arr_ return result_arr_ # 位置矩阵 else: n=time_array.shape[0] result_arr_=np.zeros((n,6),dtype=np.float64) time_float = time_array - self.starttime time_float=time_float.reshape(-1) # nx1 time_arr=np.ones((time_float.shape[0],self.polynum+1)) # nx5 time_arr[:,1]=time_float result_arr_=np.matmul(time_arr,self.A_arr) # nx5 5x6 #time_arr[0, 4] = time_arr[0, 3] * time_float ** 4 #result_arr=result_arr_ return result_arr_ # 位置矩阵 def ReconstructionSatelliteOrbit(GPSPoints_list, starttime): ''' 构建卫星轨道 args: GPSPoints_list:卫星轨道点 starttime:起算时间 ''' if len(GPSPoints_list) > 10: SatelliteOrbitModel = SatelliteOrbitFitPoly() if SatelliteOrbitModel.ReconstructionSatelliteOrbit(GPSPoints_list, starttime=starttime) is None: return None return SatelliteOrbitModel elif len(GPSPoints_list)==2: SatelliteOrbitModel = SatelliteOrbitFitPoly() if SatelliteOrbitModel.ReconstructionSatelliteOrbit(GPSPoints_list, starttime=starttime) is None: return None return SatelliteOrbitModel class DEMProcess(object): def __init__(self): pass @staticmethod def get_extent(fn): ''' 原文链接：https://blog.csdn.net/XBR_2014/article/details/85255412 ''' ds = gdal.Open(fn) rows = ds.RasterYSize cols = ds.RasterXSize # 获取图像角点坐标 gt = ds.GetGeoTransform() minx = gt[0] maxy = gt[3] maxx = gt[0] + gt[1] * rows miny = gt[3] + gt[5] * cols return (minx, maxy, maxx, miny) @staticmethod def img_mosaic(in_files, out_dem_path): # 通过两两比较大小,将最终符合条件的四个角点坐标保存 # 即为拼接图像的四个角点坐标 minX, maxY, maxX, minY = DEMProcess.get_extent(in_files[0]) for fn in in_files[1:]: minx, maxy, maxx, miny = DEMProcess.get_extent(fn) minX = min(minX, minx) maxY = max(maxY, maxy) maxX = max(maxX, maxx) minY = min(minY, miny) # 获取输出图像的行列数 in_ds = gdal.Open(in_files[0]) bands_num = in_ds.RasterCount gt = in_ds.GetGeoTransform() rows = int((maxX - minX) / abs(gt[5])) cols = int((maxY - minY) / gt[1]) # 判断栅格数据的数据类型 datatype = gdal.GDT_UInt16 # 创建输出图像 driver = gdal.GetDriverByName('GTiff') out_dem = os.path.join(out_dem_path, 'mosaic0.tif') out_ds = driver.Create(out_dem, cols, rows, bands_num, datatype) out_ds.SetProjection(in_ds.GetProjection()) gt = list(in_ds.GetGeoTransform()) gt[0], gt[3] = minX, maxY out_ds.SetGeoTransform(gt) for fn in in_files: in_ds = gdal.Open(fn) x_size = in_ds.RasterXSize y_size = in_ds.RasterYSize trans = gdal.Transformer(in_ds, out_ds, []) success, xyz = trans.TransformPoint(False, 0, 0) x, y, z = map(int, xyz) for i in range(1, bands_num + 1): data = in_ds.GetRasterBand(i).ReadAsArray() out_band = out_ds.GetRasterBand(i) out_data = out_band.ReadAsArray(x, y, x_size, y_size) data = np.maximum(data, out_data) out_band.WriteArray(data, x, y) del in_ds, out_band, out_ds @staticmethod def dem_clip(OutFilePath, DEMFilePath, SelectArea): ''' 根据选择范围裁剪DEM,并输出 agrs: outFilePath:裁剪DEM输出地址 DEMFilePath:被裁减DEM地址 SelectArea:list [(xmin,ymax),(xmax,ymin)] 框选范围 左上角，右下角 ''' DEM_ptr = gdal.Open(DEMFilePath) DEM_GeoTransform = DEM_ptr.GetGeoTransform() # 读取影像的投影变换 DEM_InvGeoTransform = gdal.InvGeoTransform(DEM_GeoTransform) SelectAreaArrayPoints = [gdal.ApplyGeoTransform(DEM_InvGeoTransform, p[0], p[1]) for p in SelectArea] SelectAreaArrayPoints = list(map(lambda p: (int(p[0]), int(p[1])), SelectAreaArrayPoints)) # 确定坐标 [(ulx, uly), (brx, bry)] = SelectAreaArrayPoints rowCount, colCount = bry - uly, brx - ulx # 输出DEM的桌面坐标转换 Out_Transfrom = list(DEM_GeoTransform) Out_Transfrom[0] = SelectArea[0][0] Out_Transfrom[3] = SelectArea[0][1] # 构建输出DEM Bands_num = DEM_ptr.RasterCount gtiff_driver = gdal.GetDriverByName('GTiff') datatype = gdal.GDT_UInt16 out_dem = gtiff_driver.Create(OutFilePath, colCount, rowCount, Bands_num, datatype) out_dem.SetProjection(DEM_ptr.GetProjection()) out_dem.SetGeoTransform(Out_Transfrom) for i in range(1, Bands_num + 1): data_band = DEM_ptr.GetRasterBand(i) out_band = out_dem.GetRasterBand(i) data = data_band.ReadAsArray(ulx, uly, colCount, rowCount) out_band.WriteArray(data) del out_dem @staticmethod def dem_merged(in_dem_path, meta_file_path, out_dem_path): ''' DEM重采样函数，默认坐标系为WGS84 agrs: in_dem_path: 输入的DEM文件夹路径 meta_file_path: 输入的xml元文件路径 out_dem_path: 输出的DEM文件夹路径 ''' # 读取文件夹中所有的DEM dem_file_paths=[os.path.join(in_dem_path,dem_name) for dem_name in os.listdir(in_dem_path) if dem_name.find(".tif")>=0 and dem_name.find(".tif.")==-1] spatialreference=osr.SpatialReference() spatialreference.SetWellKnownGeogCS("WGS84") # 设置地理坐标，单位为度 degree # 设置投影坐标，单位为度 degree spatialproj=spatialreference.ExportToWkt() # 导出投影结果 # 将DEM拼接成一张大图 mergeFile =gdal.BuildVRT(os.path.join(out_dem_path,"mergedDEM_VRT.tif"),dem_file_paths) out_DEM=os.path.join(out_dem_path,"mergedDEM.tif") gdal.Warp(out_DEM, mergeFile, format="GTiff", dstSRS=spatialproj, dstNodata=-9999, outputType=gdal.GDT_Float32) time.sleep(3) #gdal.CloseDir(out_DEM) return out_DEM @staticmethod def dem_resampled(in_dem_path,out_dem_path,samling_f): in_dem=gdal.Open(in_dem_path,gdalconst.GA_ReadOnly) width=in_dem.RasterXSize height=in_dem.RasterYSize gt=list(in_dem.GetGeoTransform()) width=width*samling_f height=height*samling_f gt=[gt[0],gt[1]/samling_f,gt[2]/samling_f,gt[3],gt[4]/samling_f,gt[5]/samling_f] driver = gdal.GetDriverByName('GTiff') output = driver.Create(out_dem_path, width,height, 1,in_dem.GetRasterBand(1).DataType) output.SetGeoTransform(gt) output.SetProjection(in_dem.GetProjection()) # 参数说明 输入数据集、输出文件、输入投影、参考投影、重采样方法(最邻近内插\双线性内插\三次卷积等)、回调函数 gdal.ReprojectImage(in_dem, output, in_dem.GetProjection(), output.GetProjection(), gdalconst.GRA_CubicSpline,0.0,0.0,) return out_dem_path @staticmethod def dem_resampled_RPC(in_dem_path,refrence_img_path,out_dem_path): in_dem=gdal.Open(in_dem_path,gdalconst.GA_ReadOnly) refrence_img=gdal.Open(refrence_img_path,gdalconst.GA_ReadOnly) width=refrence_img.RasterXSize height=refrence_img.RasterYSize gt=list(refrence_img.GetGeoTransform()) driver = gdal.GetDriverByName('GTiff') output = driver.Create(out_dem_path, width,height, 1,in_dem.GetRasterBand(1).DataType) output.SetGeoTransform(gt) output.SetProjection(refrence_img.GetProjection()) # 参数说明 输入数据集、输出文件、输入投影、参考投影、重采样方法(最邻近内插\双线性内插\三次卷积等)、回调函数 gdal.ReprojectImage(in_dem, output, in_dem.GetProjection(), output.GetProjection(), gdalconst.GRA_CubicSpline,0.0,0.0,) return out_dem_path # referencefile = gdal.Open(referencefilefilePath, gdal.GA_ReadOnly) # referencefileProj = referencefile.GetProjection() # referencefileTrans = referencefile.GetGeoTransform() # bandreferencefile = referencefile.GetRasterBand(1) # Width= referencefile.RasterXSize # Height = referencefile.RasterYSize # nbands = referencefile.RasterCount # # 创建重采样输出文件（设置投影及六参数） # driver = gdal.GetDriverByName('GTiff') # output = driver.Create(outputfilePath, Width,Height, nbands, bandreferencefile.DataType) # output.SetGeoTransform(referencefileTrans) # output.SetProjection(referencefileProj) # # 参数说明 输入数据集、输出文件、输入投影、参考投影、重采样方法(最邻近内插\双线性内插\三次卷积等)、回调函数 # gdal.ReprojectImage(inputrasfile, output, inputProj, referencefileProj, gdalconst.GRA_Bilinear,0.0,0.0,) class Orthorectification(object): """ 正射校正模块类 """ def __init__(self, configPath="./config.ymal") -> None: super().__init__() # 常量声明 # 加载配置文件 d=os.listdir(".") with open(configPath, 'r', encoding='utf-8') as f: const = f.read() self.config = yaml.load(const, Loader=yaml.FullLoader) self.config['SatelliteOribit']['StartTime'] = datetime.datetime.strptime( self.config['SatelliteOribit']['StartTime']['Value'], self.config['SatelliteOribit']['StartTime']['format']).timestamp() # 转化成UTC时间格式 self.R_e = self.config['SatelliteOribit']['ReferenceSpheroid']['Re'] # m self.R_p = self.config['SatelliteOribit']['ReferenceSpheroid']['Rp'] # m self.w_e = self.config['SatelliteOribit']['ReferenceSpheroid']['we'] # rad/s def ParseHearderFile(self, HeaderFilePath_str): """ 解析头文件，返回计算所需要的必要信息。 args: HeaderFilePath_str: 头文件地址 return: HeaderFileInfomation_json: 头文件信息包 """ with open(HeaderFilePath_str, 'r', encoding='utf-8') as fp: HeaderFile_dom_str = fp.read() HeaderFile_dom = xmltodict.parse(HeaderFile_dom_str) # 将XML转成json文本 HeaderFile_dom_json = json.loads(json.dumps(HeaderFile_dom)) # 将json字符串，转成json对象（对应于python中的dict） # 获取头文件 HeaderInformation_json = {} # 解析轨道节点，假定是GPS节点 HeaderInformation_json['GPS'] = [] # GPS 解析信息 GPSNode_Path = self.config['GPS']['GPSNode_Path'] GPSNodeNames_list = self.config['GPS']['NodeInfomation_Name'] GPSTime_Format = self.config['GPS']['Time_format'] GPSPoints = FindInfomationFromJson(HeaderFile_dom_json, GPSNode_Path) for GPSPoint in GPSPoints: GPSPoint = [ datetime.datetime.strptime(GPSPoint[GPSNodeNames_list[0]], GPSTime_Format).timestamp(), # TimeStamp float(GPSPoint[GPSNodeNames_list[1]]), # Xp float(GPSPoint[GPSNodeNames_list[2]]), # Yp float(GPSPoint[GPSNodeNames_list[3]]), # Zp float(GPSPoint[GPSNodeNames_list[4]]), # Vx float(GPSPoint[GPSNodeNames_list[5]]), # Vy float(GPSPoint[GPSNodeNames_list[6]])] # VZ HeaderInformation_json['GPS'].append(GPSPoint) # 提取成像相关信息 HeaderInformation_json['ImageInformation'] = {} # 1、开始成像时间 HeaderInformation_json['ImageInformation']['StartTime'] = datetime.datetime.strptime( FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['StartImageTime']['NodePath']), self.config['imageinfo']['StartImageTime']['Format'] ).timestamp() # 2、结束成像时间 HeaderInformation_json['ImageInformation']['EndTime'] = datetime.datetime.strptime( FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['EndImageTime']['NodePath']), self.config['imageinfo']['StartImageTime']['Format'] ).timestamp() # 3、影像的宽高 HeaderInformation_json['ImageInformation']['height'] = int( FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['ImageHeight']['NodePath'])) HeaderInformation_json['ImageInformation']['width'] = int( FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['ImageWidth']['NodePath'])) # 4、影像的近斜距 HeaderInformation_json['ImageInformation']['NearRange'] = float( FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['NearRange']['NodePath'])) # 5、入射角 HeaderInformation_json['ImageInformation']['incidenceAngle'] = { 'NearRange': float(FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['incidenceAngle']['NearRange'][ 'NodePath'])), 'FarRange': float(FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['incidenceAngle']['FarRange'][ 'NodePath'])) } # 6、成像时刻影像带宽 HeaderInformation_json['ImageInformation']['bandWidth'] = float( FindInfomationFromJson(HeaderFile_dom_json, self.config['sensor']['bandWidth']['NodePath'])) # 7、多普勒质心常数 DopplerCentroidCoefficients_list = FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo'][ 'DopplerCentroidCoefficients']['NodePath']) HeaderInformation_json['ImageInformation']['DopplerCentroidCoefficients'] = [ float(DopplerCentroidCoefficients_list[ self.config['imageinfo']['DopplerCentroidCoefficients']['DopplerCentroidCoefficients_Name'][0]]), float(DopplerCentroidCoefficients_list[ self.config['imageinfo']['DopplerCentroidCoefficients']['DopplerCentroidCoefficients_Name'][1]]), float(DopplerCentroidCoefficients_list[ self.config['imageinfo']['DopplerCentroidCoefficients']['DopplerCentroidCoefficients_Name'][2]]), float(DopplerCentroidCoefficients_list[ self.config['imageinfo']['DopplerCentroidCoefficients']['DopplerCentroidCoefficients_Name'][3]]), float(DopplerCentroidCoefficients_list[ self.config['imageinfo']['DopplerCentroidCoefficients']['DopplerCentroidCoefficients_Name'][4]]), ] # 8、波长 HeaderInformation_json['ImageInformation']['lambda'] = float( FindInfomationFromJson(HeaderFile_dom_json, self.config['sensor']['lambda']['NodePath'])) # 9、中心像元对应的中心经纬度 CenterPoint = FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['CenterImagePositon']['NodePath']) HeaderInformation_json['ImageInformation']['centerPosition'] = [ float(CenterPoint[self.config['imageinfo']['CenterImagePositon']['Value'][0]]), # 纬度 float(CenterPoint[self.config['imageinfo']['CenterImagePositon']['Value'][1]]), # 经度 ] # 10、多普勒参数参考时间 HeaderInformation_json['ImageInformation']['DopplerParametersReferenceTime'] = float( FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['DopplerParametersReferenceTime']['NodePath'])) # 11、参考斜距 HeaderInformation_json['ImageInformation']['refRange'] = float( FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['ReferenceRange']['NodePath'])) # 12、PRF HeaderInformation_json['PRF'] = float( FindInfomationFromJson(HeaderFile_dom_json, self.config['sensor']['PRF']['NodePath'])) HeaderInformation_json['Fs'] = float( FindInfomationFromJson(HeaderFile_dom_json, self.config['sensor']['Fs']['NodePath'])) # 13、中心时间 HeaderInformation_json['ImageInformation']['CenterTime'] = datetime.datetime.strptime( FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['CenterImageTime']['NodePath']), self.config['imageinfo']['CenterImageTime']['Format'] ).timestamp() #14. 影像的空间间隔分辨率 HeaderInformation_json['ImageInformation']['ImageWidthSpace']=float( FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['ImageWidthSpace']['NodePath'])) HeaderInformation_json['ImageInformation']['ImageHeightSpace']=float( FindInfomationFromJson(HeaderFile_dom_json, self.config['imageinfo']['ImageHeightSpace']['NodePath'])) # 部分需要解析 self.lamda=HeaderInformation_json['ImageInformation']['lambda'] self.PRF=HeaderInformation_json['PRF'] self.imgwidth=HeaderInformation_json['ImageInformation']['width'] self.imgheight=HeaderInformation_json['ImageInformation']['height'] self.imgstarttime=HeaderInformation_json['ImageInformation']['StartTime'] self.widthspace=HeaderInformation_json['ImageInformation']['ImageWidthSpace'] self.heightspace=HeaderInformation_json['ImageInformation']['ImageHeightSpace'] self.refrange=HeaderInformation_json['ImageInformation']['refRange'] self.nearrange=HeaderInformation_json['ImageInformation']['NearRange'] self.Fs=HeaderInformation_json['Fs']*1e6 # Mhz return HeaderInformation_json pass def LLA_to_XYZ(self, latitude, longitude, altitude): """ 经纬度转大地坐标系 args: latitude:纬度 longitude:经纬 altitude:海拔高程 refrence: https://blog.csdn.net/dulingwen/article/details/96868530 """ # 经纬度的余弦值 cosLat = math.cos(latitude * math.pi / 180) sinLat = math.sin(latitude * math.pi / 180) cosLon = math.cos(longitude * math.pi / 180) sinLon = math.sin(longitude * math.pi / 180) # WGS84坐标系的参数 rad = 6378137.0 # 地球赤道平均半径（椭球长半轴：a） f = 1.0 / 298.257224 # WGS84椭球扁率 :f = (a-b)/a C = 1.0 / math.sqrt(cosLat * cosLat + (1 - f) * (1 - f) * sinLat * sinLat) S = (1 - f) * (1 - f) * C h = altitude # 计算XYZ坐标 X = (rad * C + h) * cosLat * cosLon Y = (rad * C + h) * cosLat * sinLon Z = (rad * S + h) * sinLat return np.array([X, Y, Z]).reshape(1, 3) def XYZ_to_LLA(self, X, Y, Z): """ 大地坐标系转经纬度 适用于WGS84坐标系 args: x,y,z return: lat,lon,altitude """ # WGS84坐标系的参数 a = 6378137.0 # 椭球长半轴 b = 6356752.314245 # 椭球短半轴 ea = np.sqrt((a ** 2 - b ** 2) / a ** 2) eb = np.sqrt((a ** 2 - b ** 2) / b ** 2) p = np.sqrt(X ** 2 + Y ** 2) theta = np.arctan2(Z * a, p * b) # 计算经纬度及海拔 longitude = np.arctan2(Y, X) latitude = np.arctan2(Z + eb ** 2 * b * np.sin(theta) ** 3, p - ea ** 2 * a * np.cos(theta) ** 3) N = a / np.sqrt(1 - ea ** 2 * np.sin(latitude) ** 2) altitude = p / np.cos(latitude) - N return np.array([np.degrees(latitude), np.degrees(longitude), altitude]) def getRByColumnCode(self, c): """ 根据列号计算斜距 args: c: 列号 """ return self.R0 + c * self.delta_R pass def getTimeByLineCode(self, r): """ 根据行号计算成像时间 args: r: 行号 """ return np.matmul( np.concatenate([np.array([r ** i]).reshape(1, 1) for i in range(self.r2t_A_arr.shape[0])], axis=1), self.r2t_A_arr) def ReconstuctionTimesOflyDirectionPositionModel(self, timePoints_list): """ 构建飞行向坐标与时间的转换模型,注意这里没有调整时间的起算点。 args: timePoints_list:时间点坐标 [r,t] """ # 根据点数确定模型，最高为3次项模型 Y = np.zeros((len(timePoints_list), 1)) X = np.zeros((len(timePoints_list), 1)) for i in range(len(timePoints_list)): Y[i, 0] = timePoints_list[i][1] X[i, 0] = timePoints_list[i][0] if len(timePoints_list) == 2: # 一次项 X = np.concatenate([np.ones((len(timePoints_list), 1)), X], axis=1) A = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), Y) elif len(timePoints_list) >= 3: # 二次项 X = np.concatenate([np.ones(len(timePoints_list), 1), X, np.power(X, 2)], axis=1) A = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), Y) elif len(timePoints_list) >= 8: # 三次项 X = np.concatenate([np.ones(len(timePoints_list), 1), X, np.power(X, 2), np.power(X, 3)], axis=1) A = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), Y) self.r2t_A_arr = A def PrepareConvertSystem(self): """ 计算常量 在数据计算之前，提前处理部分常量信息 构建r,c到坐标斜距的计算公式 """ self.R0 = self.header_info['ImageInformation']['NearRange'] # 起始斜距/近斜距 self.ReconstuctionTimesOflyDirectionPositionModel([ [0, self.header_info['ImageInformation']['StartTime']], [self.header_info['ImageInformation']['height'] - 1, self.header_info['ImageInformation']['EndTime']] ]) # 构建坐标到时间的变换，这里仅用两个数据点构建 self.delta_R = self.widthspace def ConvertCoordinary(self): """ 批量求解点坐标 """ pass def Orthorectification(self, FilePath_str): """ 正射校正组件 正射校正整体可以分为两个步骤： Step 1：计算出栅格坐标对应的真实的大地坐标系坐标 Step 2：根据大地坐标系将影像重采样到大地坐标系空间中 args: FilePath_str:影像所在文件夹 """ # 分离需要校正对象 file_name_list = os.listdir(FilePath_str) header_name = [file_name for file_name in file_name_list if file_name.rfind(".meta.xml") == len(file_name) - len('.meta.xml')][0] # xml头文件，存储在list中 tiff_name = [file_name for file_name in file_name_list if file_name.rfind(".tiff") == len(file_name) - len('.tiff')] # tif影像文件，存储在list中 # 解析头文件 self.header_info = self.ParseHearderFile(os.path.join(FilePath_str, header_name)) # 构建轨道模型 self.SatelliteOrbitModel = ReconstructionSatelliteOrbit(self.header_info['GPS'], starttime=self.config['SatelliteOribit']['StartTime']) # 求解模型前计算常数变换模型 self.PrepareConvertSystem() # 批量计算空间坐标 self.ConvertCoordinary() ''' 卫星轨道定位法，修改为静态类 间接定法法，模拟SAR方法的静态方法，与类独立，方便使用多进程 ''' def getSatelliteSpaceState(time_array,SatelliteOrbitModelstarttime,SatelliteOrbitModelParameters=None,SatelliteOrbitModelName="多项式"): if SatelliteOrbitModelName=="多项式": if SatelliteOrbitModelParameters is None: return None # 返回None,表示没有结果 n=time_array.shape[0] result_arr_=np.zeros((n,6),dtype=np.float64) time_float = time_array - SatelliteOrbitModelstarttime time_float=time_float.reshape(-1) # nx1 time_arr=np.ones((time_float.shape[0],5)) # nx5 time_arr[:,1]=time_float time_arr[:,2]=time_float**2 time_arr[:,3]=time_float**3 time_arr[:,4]=time_float**4 result_arr_=np.matmul(time_arr,SatelliteOrbitModelParameters) # nx5 5x6 #time_arr[0, 4] = time_arr[0, 3] * time_float ** 4 #result_arr=result_arr_ return result_arr_ # 位置矩阵 ''' 双线性重采样 ''' def LinearSampling(ori_dem,sampling_f,offset_row=1,offset_col=1): ''' 简化双线性重采样模型 根据重采样率对原始对象进行重采样。 双线性重采样，这里去栅格所涉及的四个角点作为起算点。 模型示意图： y1 r1 y2 r y3 r2 y4 args: ori_dem:原始DEM sampling_f:采样率 int return: new_dem:采样结果 ''' # 原始dem的行列号 ori_height=ori_dem.shape[0]-offset_row ori_width=ori_dem.shape[1]-offset_col # 采样之后的影像的范围大小 new_height=int(ori_height*sampling_f) new_width=int(ori_width*sampling_f) # 获取采样之后的栅格单元对应的原DEM的行列号 row_index_list=list(range(new_height)) col_index_list=list(range(new_width)) new_dem=np.ones((new_height,new_width)) row_index_arr=(new_dem.transpose(1,0)*row_index_list).transpose(1,0) col_index_arr=new_dem*col_index_list row_index_arr=row_index_arr*1.0/sampling_f col_index_arr=col_index_arr*1.0/sampling_f # 初始化 new_dem=0 # 计算权重 row_index_arr=row_index_arr.reshape(-1) col_index_arr=col_index_arr.reshape(-1) last_row_index_arr=np.floor(row_index_arr).astype(np.int32) next_row_index_arr=np.ceil(row_index_arr).astype(np.int32) last_col_index_arr=np.floor(col_index_arr).astype(np.int32) next_col_index_arr=np.ceil(col_index_arr).astype(np.int32) # 双线性重采样模型示意图， # y1 r1 y2 # # r # # y3 r2 y4 # 计算重采样 R1=ori_dem[last_row_index_arr,last_col_index_arr]*(col_index_arr-last_col_index_arr)+ori_dem[last_row_index_arr,next_col_index_arr]*(next_col_index_arr-col_index_arr) R2=ori_dem[next_row_index_arr,last_col_index_arr]*(col_index_arr-last_col_index_arr)+ori_dem[next_row_index_arr,next_col_index_arr]*(next_col_index_arr-col_index_arr) new_dem=R1*(row_index_arr-last_row_index_arr)+R2*(next_row_index_arr-row_index_arr) # 双线性重采样 new_dem=new_dem.reshape(new_height,new_width) del R1,R2,next_row_index_arr,next_col_index_arr,last_row_index_arr,last_col_index_arr,row_index_arr,col_index_arr result_dem=np.ones((new_height,new_width,3)) result_dem[:,:,2]=new_dem result_dem[:,:,0]=(np.ones((new_height,new_width)).transpose(1,0)*list(range(new_height))).transpose(1,0) result_dem[:,:,1]=np.ones((new_height,new_width))*col_index_list return result_dem pass class IndirectOrthorectification(Orthorectification): """ 间接定位法 """ def outParameterText(self,outparameter_path): ''' 生成配置文件 ''' with open(outparameter_path,"w",encoding='utf-8') as fp: # 输出文件 fp.write("{}\n".format(self.header_info['ImageInformation']['height'])) fp.write("{}\n".format(self.header_info['ImageInformation']['width'])) fp.write("{}\n".format(self.header_info['ImageInformation']['NearRange'])) # 近斜距 fp.write("{}\n".format(self.header_info['ImageInformation']['incidenceAngle']['NearRange'])) # 近斜距入射角 fp.write("{}\n".format(self.header_info['ImageInformation']['incidenceAngle']['FarRange'])) # 远距入射角 fp.write("{}\n".format(self.config['LightSpeed'])) # 光速 fp.write("{}\n".format(self.header_info['ImageInformation']['lambda'])) fp.write("{}\n".format(self.header_info['ImageInformation']['StartTime'])) fp.write("{}\n".format(self.header_info['PRF'])) fp.write("{}\n".format(self.refrange)) fp.write("{}\n".format(self.Fs)) fp.write("{}\n".format(self.header_info['ImageInformation']['DopplerParametersReferenceTime'])) # fp.write("{}\n".format(self.widthspace)) # 多普勒系数 fp.write("{}\n".format(len(self.header_info['ImageInformation']['DopplerCentroidCoefficients']))) for i in range(len(self.header_info['ImageInformation']['DopplerCentroidCoefficients'])): fp.write("{}\n".format(self.header_info['ImageInformation']['DopplerCentroidCoefficients'][i])) fp.write("{}\n".format(1)) fp.write("{}\n".format(self.SatelliteOrbitModel.polynum)) fp.write("{}\n".format(self.SatelliteOrbitModel.get_starttime())) # X for i in range(self.SatelliteOrbitModel.polynum+1): fp.write("{}\n".format(self.SatelliteOrbitModel.A_arr[i,0])) # Y for i in range(self.SatelliteOrbitModel.polynum+1): fp.write("{}\n".format(self.SatelliteOrbitModel.A_arr[i,1])) # Z for i in range(self.SatelliteOrbitModel.polynum+1): fp.write("{}\n".format(self.SatelliteOrbitModel.A_arr[i,2])) # Vx for i in range(self.SatelliteOrbitModel.polynum+1): fp.write("{}\n".format(self.SatelliteOrbitModel.A_arr[i,3])) # Vy for i in range(self.SatelliteOrbitModel.polynum+1): fp.write("{}\n".format(self.SatelliteOrbitModel.A_arr[i,4])) # vz for i in range(self.SatelliteOrbitModel.polynum+1): fp.write("{}\n".format(self.SatelliteOrbitModel.A_arr[i,5])) # UTC参数 startTime=time.localtime(self.header_info['ImageInformation']["StartTime"]) fp.write("{}\n".format(startTime.tm_year)) fp.write("{}\n".format(startTime.tm_mon)) fp.write("{}".format(startTime.tm_mday)) self.paramterFile_path=outparameter_path def IndirectOrthorectification(self, FilePath_str,workspace_dir): """ 正射校正组件 正射校正整体可以分为两个步骤： Step 1：计算出栅格坐标对应的真实的大地坐标系坐标 Step 2：根据大地坐标系将影像重采样到大地坐标系空间中 args: FilePath_str：影像所在文件夹 dem_resampled_path：重采样后的DEM路径 lut_tif_path：输出的r,c,ti数据影像 work_temp_path: 输出的临时文件地址路径，方便存放计算临时文件 """ # 分离需要校正对象 file_name_list = os.listdir(FilePath_str) # header_name = [file_name for file_name in file_name_list if # file_name.rfind(".meta.xml") == len(file_name) - len('.meta.xml')][0] # 头文件 header_name = [file_name for file_name in file_name_list if file_name.rfind(".xml") == len(file_name) - len('.xml')][0] # 头文件 tiff_name = [file_name for file_name in file_name_list if file_name.rfind(".tiff") == len(file_name) - len('.tiff')] # 影像文件 tiff_name = [file_name for file_name in file_name_list if file_name.rfind(".tiff") == len(file_name) - len('.tiff')] # 影像文件 # 解析头文件 self.header_info = self.ParseHearderFile(os.path.join(FilePath_str, header_name)) # 构建轨道模型 self.SatelliteOrbitModel = ReconstructionSatelliteOrbit(self.header_info['GPS'], starttime=self.header_info['ImageInformation']['CenterTime']) # 构建卫星轨道模型,取第0个节点的时间 # 获取所有轨道节点时间 gpslist=np.array(self.header_info['GPS']).reshape(-1,7) # 验证轨道模型 statepoints=self.SatelliteOrbitModel.getSatelliteSpaceState(gpslist[:,0].reshape(-1)) # 计算轨道差距 statepoints_dist=statepoints-gpslist[:,1:] statepoints_dist_line=statepoints_dist[:,:3] statepoints_dist_line=np.sum(statepoints_dist_line**2,axis=1)**0.5 statepoints_dist_ver=np.sum(statepoints_dist[:,3:]**2,axis=1)**0.5 print("sata Point distance:\t{}\t -\t{}\t ".format(np.min(statepoints_dist_line),np.max(statepoints_dist_line))) # 求解模型前计算常数变换模型 self.PrepareConvertSystem() self.outParameterText(os.path.join(workspace_dir,"orth_para.txt")) def preCaldem_sar_rc(self,dem_path,ori_sar_path,work_path,taget_path): #.\SIMOrthoProgram\x64\Release\SIMOrthoProgram.exe 1 parameter_path dem_path ori_sar_path work_path taget_path # SIMOrthoProgram.exe 1 in_parameter_path in_dem_path in_ori_sar_path in_work_path in_taget_path exe_path=r".\baseTool\x64\Release\SIMOrthoProgram.exe" exe_cmd=r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5} {6} {7}".format(exe_path,1,self.paramterFile_path,dem_path,ori_sar_path,work_path ,taget_path,"\\") print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") def getRPC_incidenceAngle(self,in_dem_path,in_rpc_rc_path,out_rpc_dem_path,out_incident_angle_path,out_local_incident_angle_path): ''' std::cout << "mode 4: get RPC incident and local incident angle sar model:"; std::cout << "SIMOrthoProgram.exe 4 in_parameter_path in_dem_path in_rpc_rc_path out_rpc_dem_path out_incident_angle_path out_local_incident_angle_path"; ''' exe_path=r".\baseTool\x64\Release\SIMOrthoProgram.exe" exe_cmd=r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5} {6} {7}".format(exe_path,4,self.paramterFile_path,in_dem_path,in_rpc_rc_path,out_rpc_dem_path,out_incident_angle_path,out_local_incident_angle_path) print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") def getRPC_incidenceAngle_lon_lat(self,in_dem_path,in_gec_lon_lat_path,work_path,taget_path,out_incident_angle_path,out_local_incident_angle_path,out_incangle_geo_path,out_localincangle_geo_path): ''' std::cout << "mode 7: get RPC incident and local incident angle sar model:"; std::cout << "SIMOrthoProgram.exe 7 in_parameter_path in_dem_path in_gec_lon_lat_path work_path taget_path out_incident_angle_path out_local_incident_angle_path"; ''' cwd_path = os.getcwd() print("cwd_path:" + cwd_path) exe_path=r"{}\baseTool\x64\Release\SIMOrthoProgram.exe".format(cwd_path) exe_cmd=r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10}".format(exe_path,7,self.paramterFile_path,in_dem_path,in_gec_lon_lat_path,work_path,taget_path,out_incident_angle_path,out_local_incident_angle_path,out_incangle_geo_path,out_localincangle_geo_path) print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") def get_incidenceAngle(self,in_dem_path,in_rc_wgs84_path,out_incident_angle_path,out_local_incident_angle_path): ''' std::cout << "mode 2: get incident angle and local incident angle by rc_wgs84 and dem and statellite model:\n"; std::cout << "SIMOrthoProgram.exe 2 in_parameter_path in_dem_path in_rc_wgs84_path out_incident_angle_path out_local_incident_angle_path"; ''' exe_path=r".\baseTool\x64\Release\SIMOrthoProgram.exe" exe_cmd=r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5} {6}".format(exe_path,2,self.paramterFile_path,in_dem_path,in_rc_wgs84_path,out_incident_angle_path,out_local_incident_angle_path) print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") def get_GEC_incidenceAngle(self,in_dem_path,in_gec_lon_lat_path,out_incident_angle_path,out_local_incident_angle_path): ''' std::cout << "mode 6: get gec incident and local incident angle sar model:"; std::cout << "SIMOrthoProgram.exe 6 in_parameter_path in_dem_path in_gec_lon_lat_path out_incident_angle_path out_local_incident_angle_path"; ''' exe_path=r".\baseTool\x64\Release\SIMOrthoProgram.exe" exe_cmd=r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5} {6}".format(exe_path,6,self.paramterFile_path,in_dem_path,in_gec_lon_lat_path,out_incident_angle_path,out_local_incident_angle_path) print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") def get_GTC_SAR(self,in_rc_wgs84_path,in_ori_sar_path,out_orth_sar_path,modecode=3): ''' std::cout << "mode 3: interpolation(cubic convolution) orth sar value by rc_wgs84 and ori_sar image and model:\n "; std::cout << "SIMOrthoProgram.exe 3 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path"; ''' exe_path=r".\baseTool\x64\Release\SIMOrthoProgram.exe" exe_cmd=r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5}".format(exe_path,modecode,self.paramterFile_path,in_rc_wgs84_path,in_ori_sar_path,out_orth_sar_path) print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") def inter_Range2Geo(self,lon_lat_path , data_tiff , grid_path , space): ''' # std::cout << "mode 10"; # std::cout << "SIMOrthoProgram.exe 10 lon_lat_path data_tiff grid_path space"; ''' exe_path=r".\baseTool\x64\Release\SIMOrthoProgram.exe" exe_cmd=r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5}".format(exe_path,10,lon_lat_path , data_tiff , grid_path , space) print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") def calInterpolation_cubic_Wgs84_rc_sar_sigma(self, parameter_path, dem_rc, in_sar, out_sar): ''' # std::cout << "mode 9"; # std::cout << "SIMOrthoProgram.exe 9 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path"; ''' exe_path = r".\baseTool\x64\Release\SIMOrthoProgram.exe" exe_cmd = r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5}".format(exe_path, 9, parameter_path, dem_rc, in_sar, out_sar) print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") def calInterpolation_bil_Wgs84_rc_sar_sigma(self, parameter_path, dem_rc, in_sar, out_sar): ''' # std::cout << "mode 11"; # std::cout << "SIMOrthoProgram.exe 11 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path"; ''' exe_path = r".\baseTool\x64\Release\SIMOrthoProgram.exe" exe_cmd = r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5}".format(exe_path, 11, parameter_path, dem_rc, in_sar, out_sar) print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") def getPowerTif(self,in_ori_path,out_power_path): ''' std::cout << "mode 5: convert ori tiff to power tiff:"; std::cout << "SIMOrthoProgram.exe 5 in_ori_path out_power_path"; ''' exe_path=r".\baseTool\x64\Release\SIMOrthoProgram.exe" exe_cmd=r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} ".format(exe_path,5,in_ori_path,out_power_path) print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") def test_PSTN(self): landpoint=np.array([-1945072.5,5344083.0,2878316.0]).reshape(1,3) # 间接定位法所需的参数 Doppler_d = np.array(self.header_info['ImageInformation']['DopplerCentroidCoefficients']).reshape(-1, 1) LightSpeed = self.config['LightSpeed'] T0 = self.header_info['ImageInformation']['refRange'] * 2 / LightSpeed lamda = self.header_info['ImageInformation']['lambda'] StartTime = self.header_info['ImageInformation']['StartTime'] PRF = self.header_info['PRF'] R0 = self.header_info['ImageInformation']['NearRange'] r,c,ti=self.PSTN_block(landpoint, Doppler_d, StartTime, lamda, T0, LightSpeed, PRF, R0) pass def findInitPoint(self,sim_mask,dem_gt,sim_r_tif,sim_c_tif,sampling_f): # 根据采样点，大致计算出影像的仿射矩阵，并以此进行外扩 sampling_f=sampling_f//10 [r_ids,c_ids]=np.where(sim_mask==1) sim_r_ids=sim_r_tif[r_ids,c_ids] sim_c_ids=sim_c_tif[r_ids,c_ids] X=dem_gt[0]+dem_gt[1]*c_ids+dem_gt[2]*r_ids Y=dem_gt[3]+dem_gt[4]*c_ids+dem_gt[5]*r_ids # 首先计算X 变换式 sim_r_ids=sim_r_ids.reshape(-1,1) sim_c_ids=sim_c_ids.reshape(-1,1) X=X.reshape(-1,1) Y=Y.reshape(-1,1) point_count=sim_c_ids.shape[0] X_input=np.ones((point_count,3)) X_input[:,1]=sim_c_ids[:,0] X_input[:,2]=sim_r_ids[:,0] A_X=np.matmul(np.matmul(np.linalg.inv(np.matmul(X_input.T,X_input)), X_input.T), X).reshape(-1,1) A_Y=np.matmul(np.matmul(np.linalg.inv(np.matmul(X_input.T,X_input)), X_input.T), Y).reshape(-1,1) sar_gt=[A_X[0,0],A_X[1,0],A_X[2,0],A_Y[0,0],A_Y[1,0],A_Y[2,0]] width=self.header_info['ImageInformation']['width'] height=self.header_info['ImageInformation']['height'] width_f=self.header_info['ImageInformation']['width']/width height_f=self.header_info['ImageInformation']['height']/height sar_gt[1]=sar_gt[1]*width_f sar_gt[2]=sar_gt[2]*height_f sar_gt[4]=sar_gt[4]*width_f sar_gt[5]=sar_gt[5]*height_f # 外扩 sar_gt[0]=sar_gt[0]+sar_gt[1]*-60+sar_gt[2]*-60 sar_gt[3]=sar_gt[3]+sar_gt[4]*-60+sar_gt[5]*-60 inv_sar_gt=gdal.InvGeoTransform(sar_gt) for dem_height in [0,sim_mask.shape[0]-1]: for dem_width in [0,sim_mask.shape[1]-1]: dem_edge_x=dem_gt[0]+dem_gt[1]*dem_width+dem_gt[2]*dem_height dem_edge_y=dem_gt[3]+dem_gt[4]*dem_width+dem_gt[5]*dem_height c_width=inv_sar_gt[0]+inv_sar_gt[1]*dem_edge_x+inv_sar_gt[2]*dem_edge_y r_height=inv_sar_gt[3]+inv_sar_gt[4]*dem_edge_x+inv_sar_gt[5]*dem_edge_y if c_width<=0 or r_height<=0: continue width=width+120 if width+120=-3000)&(sim_r_tif=0)&(sim_c_tif=0 else 0 c_min=c_min-len_col if c_min-len_col>=0 else 0 r_max=r_max+len_row if r_max+len_row=0)&(r_int=0)&(c_int0: # 可以进行行前推 start_row=1 # 因为外推了一行 else: # 不可以进行行前推 start_row=0 # 此时没有外推一行 if i+block_height_width>=row_count: # 不存在行后推 if i>0: # 存在行前推， end_row=start_row+block_row_count-1 else: # 不存在行前推 end_row=start_row+block_row_count else: # 存在行后推 if i>0: # 存在行前推 end_row=start_row+block_row_count-2 else: # 不存在行前推 end_row=start_row+block_row_count-1 pass # 确定列的范围 if j>0: # 存在列前推 start_col=1 else: # 不存在列前推 start_col=0 if j+block_height_width>=col_count: # 不存在列后推 if j>0: # 存在列前推 end_col=start_col+block_col_count-1 else: # 不存在列前推 end_col=start_col+block_col_count else: # 存在列后推 if j>0: # 存在列前推 end_col=start_col+block_col_count-2 else: # 不存在列前推 end_col=start_col+block_col_count-1 pass return [start_row,start_col,end_row,end_col] def CalIntensityResampled(self,TargetPosition,ti,min_theta=0,max_theta=1): '''计算入射角 agrs: TargetPosition: 地面坐标 nxmx3 ti:对应的卫星空间位置 nxm return: Intensity:入射角阵，-9999 不参与运算 ''' # 1 # 2 0 4 # 3 # n=(n12+n23+n34+n41)/4 # n=(n24xn31)/4 # 计算 n24=n4-n2 n24=TargetPosition[:,2:,:]-TargetPosition[:,:-2,:] # 4-2 2: 5070x5068 n24=n24[1:-1,:,:] # 计算 n31=n1-n3 n31=TargetPosition[:-2,:,:]-TargetPosition[2:,:,:] # 1-3 5068x5070 n31=n31[:,1:-1,:] # 计算坡度平面向量 N=np.zeros((n31.shape[0],n31.shape[1],n31.shape[2])) # n24=[a,b,c] # n31=[d,f,g] N[:,:,0]=n24[:,:,1]*n31[:,:,2]-n24[:,:,2]*n31[:,:,1] # bg-cf N[:,:,1]=n24[:,:,2]*n31[:,:,0]-n24[:,:,0]*n31[:,:,2] # cd-ag N[:,:,2]=n24[:,:,0]*n31[:,:,1]-n24[:,:,1]*n31[:,:,0] # af-bd N=N*0.25 # 平面坡度向量 N=N.reshape(-1,3) del n24,n31 n=ti.shape[0] m=ti.shape[1] # 分块计算 theta=np.ones((n-2,m-2)).astype(bool) ti=ti[1:-1,1:-1] TargetPosition=TargetPosition[1:-1,1:-1] TargetPosition=TargetPosition.reshape(-1,3) ti=ti.reshape(-1) theta=theta.reshape(-1) # 分块计算 le_iter=4000 for i in range(0,theta.shape[0],le_iter): len_ti=le_iter if len_ti+i>theta.shape[0]: len_ti=theta.shape[0]-i end_i=i+len_ti # 计算卫星空间坐标 R_sp = self.SatelliteOrbitModel.getSatelliteSpaceState(ti[i:end_i])[:, :3]-TargetPosition[i:end_i] # 夹角=arccos((R_sp*N)/(|R_sp|*|N|)) theta_=np.sum(R_sp*N[i:end_i],axis=1)/(np.sum(R_sp**2,axis=1)*np.sum(N[i:end_i],axis=1)) theta[i:end_i]=(theta_>min_theta)&(theta_ self.header_info['ImageInformation']['height']+sampling_f*2: return False,ti[-1,0] if np.max(c)<-1*sampling_f*2 or np.min(c) > self.header_info['ImageInformation']['width']+sampling_f*2: return False,ti[-1,0] return True def ConvertCoordinary_test(self, dem_merged_path, sim_out_path,work_temp_path, flag=1): print("当前python程序，内存占用 {}G".format(psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024)) """ 批量求解点坐标的性能改进： 主要针对数据频繁开关文件读写 flag:是否想要经纬度转大地坐标系 """ start = time.time() # 间接定位法所需的参数 Doppler_d = np.array(self.header_info['ImageInformation']['DopplerCentroidCoefficients']).reshape(-1, 1) LightSpeed = self.config['LightSpeed'] T0 = self.header_info['ImageInformation']['refRange'] * 2 / LightSpeed lamda = self.header_info['ImageInformation']['lambda'] StartTime = self.header_info['ImageInformation']['StartTime'] PRF = self.header_info['PRF'] R0 = self.header_info['ImageInformation']['NearRange'] # 计算采样率 # 读取影像 dem_resampled_tiff = gdal.Open(dem_merged_path) row_count = dem_resampled_tiff.RasterYSize # 读取图像的行数 col_count = dem_resampled_tiff.RasterXSize # 读取图像的列数 im_proj = dem_resampled_tiff.GetProjection() im_geotrans = list(dem_resampled_tiff.GetGeoTransform()) gt = im_geotrans # 坐标变换参数 # 计算采样率 sampling_f=self.CalSamplingRateOfDEM(dem_resampled_tiff) # 创建目标虚拟影像 SAR_width=self.header_info['ImageInformation']['width'] SAR_height=self.header_info['ImageInformation']['height'] SAR_row_col_array=np.zeros((row_count,col_count,2),np.float32) for ii in range(0,row_count,3000): if ii==row_count: continue altii=dem_resampled_tiff.GetRasterBand(1).ReadAsArray(0, ii, col_count, 3000) if altii is None: # 长度过长，设置行数 altii=dem_resampled_tiff.GetRasterBand(1).ReadAsArray(0, ii, col_count, int(row_count-ii)) pass h=int(altii.shape[0]) w=int(altii.shape[1]) y_pixel=(np.ones((w,h))*range(h)+ii).transpose(1,0) x_pixel=np.ones((h,w))*range(w) lat=gt[0]+gt[1]*x_pixel+gt[2]*y_pixel # 经度 lon=gt[3]+gt[4]*x_pixel+gt[5]*y_pixel # 纬度 #del x_pixel,y_pixel #gc.collect() [X, Y, Z] = OrthoAuxData.LLA2XYZM(lat.reshape(-1,1), lon.reshape(-1,1), altii.reshape(-1,1)) # 计算错误 TargetPosition = np.ones((int(h*w), 3), dtype=np.float32) TargetPosition[:, 0] = X.reshape(-1) TargetPosition[:, 1] = Y.reshape(-1) TargetPosition[:, 2] = Z.reshape(-1) r, c, ti = self.PSTN_M(TargetPosition, Doppler_d, StartTime, lamda, T0, LightSpeed, PRF, R0) # 间接定位法，核心求解代码 SAR_row_col_array[ii:ii+h,:,0]=r.reshape(-1,col_count) SAR_row_col_array[ii:ii+h,:,1]=c.reshape(-1,col_count) gtiff_driver = gdal.GetDriverByName('GTiff') SAR_row_col=gtiff_driver.Create(sim_out_path.replace(".tif","rc.tif"),col_count,row_count,2,gdal.GDT_Float32) SAR_row_col.SetGeoTransform(im_geotrans) # 写入仿射变换参数 SAR_row_col.SetProjection(im_proj) # 写入投影 # 写入tiff SAR_row_col.GetRasterBand(1).WriteArray(SAR_row_col_array[:,:,0]) SAR_row_col.GetRasterBand(2).WriteArray(SAR_row_col_array[:,:,1]) print(sim_out_path.replace(".tif","rc.tif")) return sim_out_path.replace(".tif","rc.tif") pass def CalSamplingRateOfDEM(self,DEM_tiff): '''计算DEM的重采样率f，必须在解析头文件之后，参考陈尔学 博士学位论文 《 星载合成孔径雷达影像正射校正方法研究 》 args: DEM_tiff:tif DEM的影像 return: f:采样率，向上取整 ''' # 获取DEM_tiff的高宽、投影、仿射矩阵 Trans_array=list(DEM_tiff.GetGeoTransform()) # 仿射变换矩阵 width=DEM_tiff.RasterXSize height=DEM_tiff.RasterYSize # x=trans[0]+trans[1]*c+trans[2]*r # y=trans[3]+trans[4]*c+trans[5]*r # 获得其高宽分辨率 delta_x=Trans_array[1]+Trans_array[2] # x轴上的分辨率 经度 delta_y=Trans_array[4]+Trans_array[5] # y轴上的分辨率 纬度 # 因为单位为度，所以需要转换为 米 y_resultion=abs(111194.926644558737*delta_y) x_resultion=abs(111194.926644*math.cos(Trans_array[3])*delta_x) resultion=x_resultion if x_resultion>y_resultion else y_resultion delta_R=self.delta_R # 斜距分辨率 incidenceAngle=self.header_info['ImageInformation']['incidenceAngle']["NearRange"] # 入射角 f=math.sin(incidenceAngle*math.pi/180)*1.4142135623730951*resultion/delta_R f=math.ceil(f) return f def PSTN_M2(self, TargetPostion, Doppler_d, StartTime, lamda, T0, LightSpeed, PRF, R0): """ 间接求解方法，使用矩阵方法,分段计算，方便加速 args: TargetPosition：nparray,nx3 地面坐标 Doppler_d:多普勒系数 -1，1 lamda:波长 T0:多普勒参考时间 LightSpeed:光速 return: rc:nparray shape nx2 行列号 """ # n = int(TargetPostion.shape[0]) # 坐标点数据 nx3 r=np.zeros((n,1)) c=np.zeros((n,1)) ti=np.zeros((n,1),dtype=np.float64) len_iter=2000 pool= ThreadPoolExecutor() pool_ls=[] SatelliteOrbitModelstarttime=self.SatelliteOrbitModel.get_starttime() SatelliteOrbitModelParameters=self.SatelliteOrbitModel.A_arr SatelliteOrbitModelName=self.SatelliteOrbitModel.modelName delta_R=self.delta_R for i in range(0, n, len_iter): start_i=i len_block=len_iter if start_i+len_block>n: len_block=n-start_i end_i=start_i+len_block temp_TargetPosition=TargetPostion[start_i:end_i,:] pool_ls.append( pool.submit(PSTN_block_MuilThread,copy.deepcopy(start_i),copy.deepcopy(end_i), temp_TargetPosition.copy(), copy.deepcopy(Doppler_d), copy.deepcopy(StartTime), copy.deepcopy(lamda), copy.deepcopy(T0), copy.deepcopy(LightSpeed), copy.deepcopy(PRF),copy.deepcopy(R0), copy.deepcopy(delta_R), copy.deepcopy(SatelliteOrbitModelstarttime), SatelliteOrbitModelParameters.copy(), copy.deepcopy(SatelliteOrbitModelName),)) pool.close() pool.join() for p in pool_ls: r_temp, c_temp, ti_temp,start_i,end_i=p.get() r[start_i:end_i,:]=r_temp c[start_i:end_i,:]=c_temp ti[start_i:end_i,:]=ti_temp return r,c,ti def PSTN_M(self, TargetPostion, Doppler_d, StartTime, lamda, T0, LightSpeed, PRF, R0): """ 间接求解方法，使用矩阵方法,分段计算，方便加速 args: TargetPosition：nparray,nx3 地面坐标 Doppler_d:多普勒系数 -1，1 lamda:波长 T0:多普勒参考时间 LightSpeed:光速 return: rc:nparray shape nx2 行列号 """ # n = int(TargetPostion.shape[0]) # 坐标点数据 nx3 r=np.zeros((n,1)) c=np.zeros((n,1)) ti=np.zeros((n,1),dtype=np.float64) init_size=3000 for i in range(0,n,init_size): start_i=i len_block=init_size if start_i+len_block>n: len_block=n-start_i end_i=start_i+len_block temp_TargetPosition=TargetPostion[start_i:end_i,:] r_temp, c_temp, ti_temp=self.PSTN_block(temp_TargetPosition, Doppler_d, StartTime, lamda, T0, LightSpeed, PRF, R0) r[start_i:end_i,:]=r_temp.copy() c[start_i:end_i,:]=c_temp.copy() ti[start_i:end_i,:]=ti_temp.copy() return r,c,ti def PSTN_block(self, TargetPostion, Doppler_d, StartTime, lamda, T0, LightSpeed, PRF, R0): ''' 间接求解方法，使用矩阵方法 args: TargetPosition：nparray,nx3 地面坐标 Doppler_d:多普勒系数 -1，1 lamda:波长 T0:多普勒参考时间 LightSpeed:光速 return: rc:nparray shape nx2 行列号 ''' n = TargetPostion.shape[0] # 坐标点数据 nx3 ti = np.ones((n, 1),dtype=np.float64) * StartTime # 初始值 delta_R = self.delta_R delta_t = 1 / PRF R=np.zeros((n,1)) Rs = np.zeros((n, 3)) R_sc = np.zeros((n, 3)) # 地面-卫星空间位置矢量 V_sc = np.zeros((n, 3)) # 地面-卫星空间速率矢量 FdNumerical = np.ones((n, 1)) # 多普勒公式的数值解法 FdTheory = np.ones((n, 1)) # 多普勒公式的理论解法 FdTheory_grad = np.ones((n, 1)) # 多普勒公式的导数 # 多普勒参数系数 Doppler_d = Doppler_d.reshape(-1, 1) # 5x1 TSR=np.ones((n, 1)) # nx1 TSRx = np.ones((n, Doppler_d.shape[0])) # nx5 inc_t = np.ones((n, 1)) # 计算增加的导数 dt = 0.0001 # 数值法求解多普勒距离公式的导数 dt_rec=1/dt # 默认最高迭代100次，如果还未求解出结果，就将所有的未收敛的时间值记录为0 ，注意时间不能为0 for i in range(100): tempT = ti + dt # 计算增加时间 # ------- 求解导数 ---------------- #del Rs,R_sc,V_sc,TSR,TSRx,satelliteSpaceState,R satelliteSpaceState = self.SatelliteOrbitModel.getSatelliteSpaceState(tempT) # 卫星坐标与速度 nx6 Rs = satelliteSpaceState[:, :3] R_sc = satelliteSpaceState[:, :3] - TargetPostion # 地面-卫星空间位置矢量 V_sc = satelliteSpaceState[:, 3:] # 地面-卫星空间位置矢量 R =np.sqrt(np.sum(R_sc ** 2, axis=1)).reshape(-1, 1) # nx1 #FdTheory_grad = (-2 / lamda) * (1 / R) * np.sum(R_sc * V_sc, axis=1).reshape(-1, 1) # nx1 FdTheory_grad = -2*np.reciprocal(R*lamda) * np.sum(R_sc * V_sc, axis=1).reshape(-1, 1) # nx1 # 获取卫星轨道状态信息 satelliteSpaceState = self.SatelliteOrbitModel.getSatelliteSpaceState(ti) # 卫星坐标与速度 nx6 # 卫星轨道坐标 Rs =satelliteSpaceState[:, :3] # 卫星相对地面的位置向量、卫星速度向量 R_sc=Rs - TargetPostion # 地面-卫星空间位置矢量 V_sc =satelliteSpaceState[:, 3:] # 地面-卫星空间位置矢量 # 斜距 R = np.sqrt(np.sum(R_sc ** 2, axis=1)).reshape(-1, 1) # nx1 # 根据多普勒数值求解公式求解多普勒频移值 TSR= R * (2 / LightSpeed) - T0 # nx1 TSRx[:,0]=TSR[:,0]**0 # nx5 TSRx[:,1]=TSR[:,0]**1 TSRx[:,2]=TSR[:,0]**2 TSRx[:,3]=TSR[:,0]**3 TSRx[:,4]=TSR[:,0]**4 # 根据多普勒物理公式求解求解多普勒频移值 FdTheory=-2*np.reciprocal(lamda*R)*np.sum(R_sc * V_sc, axis=1).reshape(-1, 1) # nx1 FdNumerical= np.matmul(TSRx, Doppler_d) # nx1 FdTheory_grad = np.reciprocal((FdTheory_grad - FdTheory) * dt_rec) # nx1 inc_t = (FdTheory - FdNumerical) * FdTheory_grad #(FdTheory - FdNumerical) * (1 / FdTheory_grad) # nx1 # inc_t = inc_t - inc_t * ti_mask # 允许继续迭代 nx1 if np.max(np.abs(inc_t)) < delta_t * 0.001: break ti = ti - inc_t # 更新坐标 #del Rs,R_sc,V_sc,R,satelliteSpaceState,FdTheory_grad,FdNumerical,inc_t #print("当前python程序，{} step， inc_t {} 求解内存占用 {}G".format(i,np.max(inc_t),psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024)) # 计算行号 ti_=copy.deepcopy(ti.copy()) r_ = (ti_ - StartTime)*PRF delta_t_light=delta_t*LightSpeed r=copy.deepcopy(r_.copy()) del r_ c_ = (R - R0) * (1 / delta_R) c=copy.deepcopy(c_.copy()) del c_ ti = copy.deepcopy(ti_.copy()) del ti_ return r, c, ti pass def sim_intensity(self, in_tif_path, dem_resampled_path, lut_tif_path, incidence_tif_path, sim_out_path): """ 局部入射角计算 原方法存在的问题：多次文件读写 """ dem_resampled_tiff = gdal.Open(dem_resampled_path) # 获取重采样的DEM 数据，为了减少文件读写操作，这里做了持久化 im_proj = dem_resampled_tiff.GetProjection() # 构建投影 im_geotrans = list(dem_resampled_tiff.GetGeoTransform()) # 构建投影变换对象 [a1, b1, c1, a2, b2, c2] = im_geotrans # 坐标变换参数：c1=0,b2=0 start = time.time() # 度量算法 dem_resampled = dem_resampled_tiff.GetRasterBand(1).ReadAsArray(0, 0, dem_resampled_tiff.RasterXSize, dem_resampled_tiff.RasterYSize) del dem_resampled_tiff # 释放内存 incidence = dem_resampled * 0.0 + 1.8 # 读取时间 in_ds = gdal.Open(lut_tif_path) # 读取间接法 row_cnt = in_ds.RasterYSize col_cnt = in_ds.RasterXSize t_arr = in_ds.GetRasterBand(3).ReadAsArray(0, 0, col_cnt, row_cnt).astype(np.float64) # 获取时间阵列 t_arr = t_arr[1:-1, 1:-1] + self.header_info['ImageInformation']['StartTime'] # 去除最外围一圈时间数据 row_data = in_ds.GetRasterBand(1).ReadAsArray(0, 0, col_cnt, row_cnt) # 读取行号 col_data = in_ds.GetRasterBand(2).ReadAsArray(0, 0, col_cnt, row_cnt) # 读取列号 row_data = np.round(row_data) col_data = np.round(col_data) row_data = row_data.astype(np.int16) col_data = col_data.astype(np.int16) del in_ds # 释放内存 c = range(0, col_cnt) c = np.array(c).reshape(-1, 1) for r0 in range(1, row_cnt - 1): T0 = np.zeros((col_cnt, 3)) # colcount-2 ，3 T0[:, 0] = a1 + b1 * c[:, 0] + c1 * r0 T0[:, 1] = a2 + b2 * c[:, 0] + c2 * r0 T0[:, 2] = dem_resampled[r0, :] T0 = np.concatenate(OrthoAuxData.LLA2XYZM(T0[:, 0].reshape(-1, 1), T0[:, 1].reshape(-1, 1), T0[:, 2].reshape(-1, 1)), axis=1) n1 = np.concatenate(OrthoAuxData.LLA2XYZM(a1 + b1 * c + c1 * (r0 - 1), a2 + b2 * c + c2 * (r0 - 1), dem_resampled[r0 - 1, :]), axis=1)[1:-1, :] - T0[1:-1, :] # colcount-2 ，3 n3 = np.concatenate(OrthoAuxData.LLA2XYZM(a1 + b1 * c + c1 * (r0 + 1), a2 + b2 * c + c2 * (r0 + 1), dem_resampled[r0 + 1, :]), axis=1)[1:-1, :] - T0[1:-1, :] # colcount-2 ，3 n2 = T0[:-2, :] - T0[1:-1, :] n4 = T0[2:, :] - T0[1:-1, :] T0 = T0[1:-1, :] n31 = n1 - n3 n42 = n2 - n4 N = XYZOuterM2(n31, n42) / 4 ti = t_arr[r0 - 1, :].reshape(-1, 1) # rowcount-2,colcount-2 satelliteSpaceState = self.SatelliteOrbitModel.getSatelliteSpaceState(ti) Rps = satelliteSpaceState[:, :3] - T0 # nx3 # 计算向量模之积 N_Norm = np.sqrt(np.sum(N * N, axis=1)).reshape(-1, 1) Rps_Norm = np.sqrt(np.sum(Rps * Rps, axis=1)).reshape(-1, 1) N_Norm = np.concatenate((N_Norm, N_Norm, N_Norm), axis=1) Rps_Norm = np.concatenate((Rps_Norm, Rps_Norm, Rps_Norm), axis=1) N = N / N_Norm Rps = Rps / Rps_Norm # 计算向量相乘 N_Rps = np.sum(N * Rps, axis=1) # col_count-2 theta = np.arccos(N_Rps) # 入射角 incidence[r0, 1:-1] = theta incidence = np.degrees(incidence) print('整张图入射角计算完成，耗时{}秒'.format(time.time() - start)) ImageHandler.write_img(incidence_tif_path, im_proj, im_geotrans, incidence) # 读取原始影像相关信息 in_tif = gdal.Open(in_tif_path) col_num = in_tif.RasterXSize row_num = in_tif.RasterYSize im_proj1 = in_tif.GetProjection() # 构建投影 im_geotrans1 = in_tif.GetGeoTransform() out_arr = np.zeros((row_num, col_num), dtype=np.float32) for h in range(row_cnt): r_min = min(row_data[h]) r_max = max(row_data[h]) c_min = min(col_data[h]) c_max = max(col_data[h]) if 0 <= r_min < row_num or 0 <= r_max < row_num or (r_min < 0 and r_max >= row_num): if 0 <= c_min < col_num or 0 <= c_max < col_num or (c_min < 0 and c_max >= col_num): for w in range(col_cnt): r = row_data[h, w] c = col_data[h, w] if 0 <= r < row_num and 0 <= c < col_num: if 0 <= incidence[h, w] < 90: out_arr[r, c] += 1 K = 10 out_arr = K * out_arr max_delta = np.max(out_arr) min_delta = np.min(out_arr) out_arr = (out_arr - min_delta) / (max_delta - min_delta) ImageHandler.write_img(sim_out_path, im_proj1, im_geotrans1, out_arr) @staticmethod def sar_intensity_synthesis(sar_in_tif, sar_out_tif): # 获取SLC格式SAR影像的相关信息 in_ds = gdal.Open(sar_in_tif) bands_num = in_ds.RasterCount rows = in_ds.RasterYSize columns = in_ds.RasterXSize proj = in_ds.GetProjection() geotrans = in_ds.GetGeoTransform() datatype = gdal.GDT_Float32 # 创建输出的SAR强度图 gtiff_driver = gdal.GetDriverByName('GTiff') out_ds = gtiff_driver.Create(sar_out_tif, columns, rows, 1, datatype) # 输出SAR强度图 out_data=np.zeros((rows,columns)) for i in range(0,rows,10000): for j in range(0,columns,10000): len_row=10000 len_col=10000 if i+len_row>rows: len_row=rows-i if j+len_col>columns: len_col=columns-j end_i=i+len_row end_j=j+len_col in_data1 = in_ds.GetRasterBand(1).ReadAsArray(j, i, len_col, len_row).astype(np.float32) in_data2 = in_ds.GetRasterBand(2).ReadAsArray(j, i, len_col, len_row).astype(np.float32) out_data[i:end_i,j:end_j] = np.log10((in_data1**2 + in_data2**2)+1)*10 out_ds.GetRasterBand(1).WriteArray(out_data) del in_ds, out_ds @staticmethod def img_match(sim_sar_path, ori_sar_path,work_sapce_path): ''' 影像匹配并根据配准点 args: sim_sar_path:模拟SAR强度图 orig_sar_path: 原始SAR强度图 out_path:配准结果输出点 work_space_path: 工作目录 return: point_correct: 点集映射结果 raise： None: 出现错误 ''' match_points_path=ImageMatchClass.ImageMatch(ori_sar_path,sim_sar_path,work_sapce_path) match_points_path=os.path.join(work_sapce_path,"GridMatch.npy") match_result_path=os.path.join(work_sapce_path,"match_result.npy") match_result_image_path=os.path.join(work_sapce_path,"match_result_show.tif") ImageMatchClass.DeNoiseMatch(match_points_path,match_result_path) ImageMatchClass.DrawMatchResult(ori_sar_path,sim_sar_path,match_result_path,match_result_image_path) return match_result_path def IndireOrtho_smi(self,sim_SAR_rc_path,ori_SAR_path,tar_sar_path): ''' 测试间接定位法计算的影像坐标是否有问题。 方法核心，插值出间接定位法坐标处的 实部和虚部，不考虑其他区域 使用插值方法：等权反距离权重法 考虑实际的数据量较大，为节省内存，使用分块处理方法。 块大小为1000。 ''' sim_rc_tif=gdal.Open(sim_SAR_rc_path) sim_width=sim_rc_tif.RasterXSize sim_height=sim_rc_tif.RasterYSize sim_r_tif=sim_rc_tif.GetRasterBand(1).ReadAsArray(0,0,sim_width,sim_height) sim_c_tif=sim_rc_tif.GetRasterBand(2).ReadAsArray(0,0,sim_width,sim_height) # 构建掩膜 sim_mask=(sim_r_tif>0)&(sim_r_tif0)&(sim_c_tif模拟 args: point_arr:[模拟row,模拟col，原始row,原始col] ''' # 由模拟->原始 point_cnt = point_arr.shape[0] # rs_arr = point_arr[:,2].reshape(point_cnt, 1) # 原始 行 cs_arr = point_arr[:,3].reshape(point_cnt, 1) rm_arr = point_arr[:,0].reshape(point_cnt, 1)-3000 # 模拟 行 还原偏移 cm_arr = point_arr[:,1].reshape(point_cnt, 1) X = np.concatenate( [np.ones((point_cnt, 1), dtype=np.float64), rm_arr, cm_arr, rm_arr ** 2, cm_arr ** 2, rm_arr * cm_arr], axis=1) sim2orirc_arr=np.ones((2,6)) sim2orirc_arr[0,:] = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), rs_arr).reshape(1, 6) sim2orirc_arr[1,:] = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), cs_arr).reshape(1, 6) # 翻转 sim->ori X = np.concatenate( [np.ones((point_cnt, 1), dtype=np.float64), rs_arr, cs_arr, rs_arr ** 2, cs_arr ** 2, rs_arr * cs_arr], axis=1) ori2simrc_arr=np.ones((2,6)) ori2simrc_arr[0,:] = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), rm_arr).reshape(1, 6) ori2simrc_arr[1,:] = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), cm_arr).reshape(1, 6) return ori2simrc_arr.copy(),sim2orirc_arr.copy() @staticmethod def ReSamplingSAR(ori_sar_path,out_sar_path,lon_lat_point_match_path,gt): ''' 根据影像点匹配点，插值得到正射图像 args: ori_sar_path: 输入sar地址 out_sar_path: 正射结果地址 lon_lat_point_match_path: CoordinateCorrection计算结果的路径 gt: 仿射矩阵 return: out_sar_path ''' ori_sar=gdal.Open(ori_sar_path) ori_height=ori_sar.RasterYSize # 获得高度 ori_width=ori_sar.RasterXSize # 获得宽度 # 开始插值 spatialreference=osr.SpatialReference() spatialreference.SetWellKnownGeogCS("WGS84") # 设置地理坐标，单位为度 degree # 设置投影坐标，单位为度 degree spatialproj=spatialreference.ExportToWkt() # 导出投影结果 gtiff_driver = gdal.GetDriverByName('GTiff') out_sar=gtiff_driver.Create(out_sar_path, ori_width, ori_height, 2, gdal.GDT_Float32) out_sar.SetGeoTransform(gt) # 写入仿射变换参数 out_sar.SetProjection(spatialproj) # 写入投影 ab=np.array([gt[1],gt[2],gt[4],gt[5]]).reshape(2,2) ab=np.np.matmul(np.linalg.inv(np.np.matmul(ab.T,ab)),ab.T) # 解算仿射矩阵，为后面的插值做准备 lon_lat_file_list=os.listdir(lon_lat_point_match_path) block_info_json=[] for lon_lat_file_name in lon_lat_file_list: lon_lat=np.load(os.path.join(lon_lat_point_match_path,lon_lat_file_name)) ori_r=lon_lat[:,0] # 原始影像的行坐标 ori_c=lon_lat[:,1] # 原始影像的列坐标 lon=lon_lat[:,2].copy() lat=lon_lat[:,3].copy() # 根据仿射矩阵阶段投影变换结果 lon=lon-gt[0] lat=lat-gt[3] lon_lat[:,2]=ab[0,0]*lon+ab[0,1]*lat # 行 lon_lat[:,3]=ab[1,0]*lon+ab[1,1]*lat # 列 block_info_json.append({ 'ori':[np.min(lon_lat[:,0]),np.max(lon_lat[:,0]),np.min(lon_lat[:,1]),np.max(lon_lat[:,1])],# r0 r1 c0 c1 'ortho':[np.min(lon_lat[:,2]),np.max(lon_lat[:,2]),np.min(lon_lat[:,3]),np.max(lon_lat[:,3])], 'path':os.path.join(lon_lat_point_match_path,lon_lat_file_name) }) for i in range(0,ori_height,6000): for j in range(0,ori_width,6000): # 判断需要检索的行列 len_row=6000 len_col=6000 if i+len_row>ori_height: len_row=ori_height-i if j+len_col>ori_width: len_col=ori_width-j end_i=i+len_row end_j=j+len_col select_region=[i-1000,end_i+1000,j-1000,end_i+1000] # 根据求交原则，解出需要可以参与插值计算的 方块 inter_block=[] for block in block_info_json: block_ortho =block['ortho'] if block_ortho[0]>select_region[0] and block_ortho[1]select_region[2] and block_ortho[3]模拟SAR # X=a0+a1*r+a2*c+a5*r*c+a3*r^2+a4*c^2 # Y=b0+b1*r+b2*c+b5*r*c+b3*r^2+b4*r^2 point_arr=np.load(point_match_path) rcori2simrc_arr,sim2orirc_arr=IndirectOrthorectification.CoordinateTransformation(point_arr) # 根据行列号裁剪DEM sim_rc_tif=gdal.Open(sim_out_path) sim_width=sim_rc_tif.RasterXSize sim_height=sim_rc_tif.RasterYSize sim_r_tif=sim_rc_tif.GetRasterBand(1).ReadAsArray(0,0,sim_width,sim_height) sim_c_tif=sim_rc_tif.GetRasterBand(2).ReadAsArray(0,0,sim_width,sim_height) # 计算实际行列号 r=sim2orirc_arr[0,0]+sim2orirc_arr[0,1]*sim_r_tif+sim2orirc_arr[0,2]*sim_c_tif+sim2orirc_arr[0,3]*sim_r_tif*sim_r_tif+sim2orirc_arr[0,4]*sim_c_tif*sim_c_tif+sim2orirc_arr[0,5]*sim_r_tif*sim_c_tif c=sim2orirc_arr[1,0]+sim2orirc_arr[1,1]*sim_r_tif+sim2orirc_arr[1,2]*sim_c_tif+sim2orirc_arr[1,3]*sim_r_tif*sim_r_tif+sim2orirc_arr[1,4]*sim_c_tif*sim_c_tif+sim2orirc_arr[1,5]*sim_r_tif*sim_c_tif sim_r_tif=r sim_c_tif=c del r,c # sim_mask=(sim_r_tif>=-1)&(sim_r_tif=-1)&(sim_c_tif原始 ori_rc=np.ones((r_ids.size,2)) ori_rc[:,0]=r_arr_trans[0,0]+r_ids*r_arr_trans[0,1]+c_ids*r_arr_trans[0,2]+r_arr_trans[0,3]*r_ids ** 2+r_arr_trans[0,4]*c_ids ** 2+r_arr_trans[0,5]*r_ids * c_ids ori_rc[:,1]=c_arr_trans[0,0]+r_ids*c_arr_trans[0,1]+c_ids*c_arr_trans[0,2]+c_arr_trans[0,3]*r_ids ** 2+c_arr_trans[0,4]*c_ids ** 2+c_arr_trans[0,5]*r_ids * c_ids # 插值x坐标 sim_x_data=gdal.Open(sim_sar_x_path) sim_x_arr=sim_x_data.GetRasterBand(1).ReadAsArray(0,0 , sim_x_data.RasterXSize, sim_x_data.RasterYSize) # 块 高程文件 sim_x=sim_x_arr[r_ids,c_ids] del sim_x_arr,sim_x_data tree=KDTree(ori_rc) # 原始sar点对应的坐标 gtiff_driver = gdal.GetDriverByName('GTiff') SAR_row_col=gtiff_driver.Create(ori_xyz_path,width,height,3,gdal.GDT_Float32) # 准备插值--分块插值 blocksize=5000 SAR_ori=np.zeros((height,width)) for i in range(0,height,blocksize): for j in range(0,width,blocksize): start_row=i start_col=j len_row=blocksize len_col=blocksize if(start_row+len_row>height): len_row=height-start_row if(start_col+len_col>width): len_col=width-start_col end_row=start_row+len_row end_col=start_col+len_col grid_xy= np.mgrid[start_row:end_row:1, start_col:end_col:1] grid_xy=grid_xy.reshape(2,-1).transpose(1,0) distance_weight,ori_select=tree.query(grid_xy,k=10) distance_weight=np.reciprocal(distance_weight) distance_weight=np.multiply(distance_weight.transpose(1,0),np.reciprocal(np.sum(distance_weight,axis=1))) distance_weight=distance_weight.transpose(1,0) SAR_ori[start_row:end_row,start_col:end_col]=np.sum( np.multiply(distance_weight,sim_x[ori_select]),axis=1).reshape(len_row,len_col) SAR_row_col.GetRasterBand(1).WriteArray(SAR_ori) # 插值y坐标 sim_y_data=gdal.Open(sim_sar_y_path) sim_y_arr=sim_y_data.GetRasterBand(1).ReadAsArray(0,0 , sim_y_data.RasterXSize, sim_y_data.RasterYSize) # 块 高程文件 sim_y=sim_y_arr[r_ids,c_ids] del sim_y_arr,sim_y_data,sim_x # 准备插值--分块插值 blocksize=5000 SAR_ori=np.zeros((height,width)) for i in range(0,height,blocksize): for j in range(0,width,blocksize): start_row=i start_col=j len_row=blocksize len_col=blocksize if(start_row+len_row>height): len_row=height-start_row if(start_col+len_col>width): len_col=width-start_col end_row=start_row+len_row end_col=start_col+len_col grid_xy= np.mgrid[start_row:end_row:1, start_col:end_col:1] grid_xy=grid_xy.reshape(2,-1).transpose(1,0) distance_weight,ori_select=tree.query(grid_xy,k=10) distance_weight=np.reciprocal(distance_weight) distance_weight=np.multiply(distance_weight.transpose(1,0),np.reciprocal(np.sum(distance_weight,axis=1))) distance_weight=distance_weight.transpose(1,0) SAR_ori[start_row:end_row,start_col:end_col]=np.sum(np.multiply(distance_weight,sim_y[ori_select]),axis=1).reshape(len_row,len_col) SAR_row_col.GetRasterBand(2).WriteArray(SAR_ori) # 插值z坐标 sim_z_data=gdal.Open(sim_sar_z_path) sim_z_arr=sim_z_data.GetRasterBand(1).ReadAsArray(0,0 , sim_z_data.RasterXSize, sim_z_data.RasterYSize) # 块 高程文件 sim_z=sim_z_arr[r_ids,c_ids] del sim_z_arr,sim_z_data,sim_y # 准备插值--分块插值 blocksize=5000 SAR_ori=np.zeros((height,width)) for i in range(0,height,blocksize): for j in range(0,width,blocksize): start_row=i start_col=j len_row=blocksize len_col=blocksize if(start_row+len_row>height): len_row=height-start_row if(start_col+len_col>width): len_col=width-start_col end_row=start_row+len_row end_col=start_col+len_col grid_xy= np.mgrid[start_row:end_row:1, start_col:end_col:1] grid_xy=grid_xy.reshape(2,-1).transpose(1,0) distance_weight,ori_select=tree.query(grid_xy,k=10) distance_weight=np.reciprocal(distance_weight) distance_weight=np.multiply(distance_weight.transpose(1,0),np.reciprocal(np.sum(distance_weight,axis=1))) distance_weight=distance_weight.transpose(1,0) SAR_ori[start_row:end_row,start_col:end_col]=np.sum(np.multiply(distance_weight,sim_z[ori_select]),axis=1).reshape(len_row,len_col) SAR_row_col.GetRasterBand(3).WriteArray(SAR_ori) return ori_xyz_path def outResamplingParas(self,work_space_file,ortho_sar_xyz_pow_path,sar_xyz_path,sar_incidence_path,ortho_sar_xyz_path,ortho_sar_incidence_path,out_sar_file_list,dem_clip_path,ori2sim,sim2ori): ''' sar_xyz_path:xyz插值结果 sar_incidence_path: 入射角 ortho_sar_xyz_path： 重采样结果文件夹 out_sar_file_list：需要重采样的文件 rc_trans_arr：变换结果 ''' with open(os.path.join(work_space_file,"resample_para.txt"),'w',encoding='utf-8') as fp: fp.write("{}\n".format(dem_clip_path)) fp.write("{}\n".format(os.path.join(work_space_file,"ori_sim.tif"))) fp.write("{}\n".format(sar_xyz_path)) fp.write("{}\n".format(sar_incidence_path)) fp.write("{}\n".format(ortho_sar_incidence_path)) fp.write("{}\n".format(ortho_sar_incidence_path.replace("orth_sar_incidence.tif","orth_sar_local_incidence.tif"))) fp.write("{}\n".format(ortho_sar_xyz_path)) fp.write("{}\n".format(len(out_sar_file_list))) for i in range(len(out_sar_file_list)): file_name=os.path.split(out_sar_file_list[i])[1] fp.write("{}\n".format(out_sar_file_list[i])) fp.write("{}\n".format(os.path.join(ortho_sar_xyz_path,file_name))) fp.write("{}\n".format(os.path.join(ortho_sar_xyz_pow_path,file_name))) #os.path.split(in_tif_path)[1] fp.write("{}\n".format(12)) for i in range(6): fp.write("{}\n".format(ori2sim[0,i])) for i in range(6): fp.write("{}\n".format(ori2sim[1,i])) fp.write("{}\n".format(12)) for i in range(6): fp.write("{}\n".format(sim2ori[0,i])) for i in range(6): fp.write("{}\n".format(sim2ori[1,i])) return os.path.join(work_space_file,"resample_para.txt") def Resampling(self,out_sar_xyz_path,inv_gt,gt,ori_resampling_file_path): out_xyz=gdal.Open(out_sar_xyz_path) ori_height=out_xyz.RasterYSize # 获得高度 ori_width=out_xyz.RasterXSize # 获得宽度 out_xyz_x=out_xyz.GetRasterBand(1).ReadAsArray(0,0 , ori_width, ori_height) out_xyz_y=out_xyz.GetRasterBand(2).ReadAsArray(0,0 , ori_width, ori_height) # 计算对应的行列号 out_xyz_rc=np.ones((ori_height,ori_width,2)) out_xyz_rc[:,:,1]=inv_gt[0]+inv_gt[1]*out_xyz_x+inv_gt[2]*out_xyz_y out_xyz_rc[:,:,0]=inv_gt[3]+inv_gt[4]*out_xyz_x+inv_gt[5]*out_xyz_y del out_xyz_x,out_xyz_y,out_xyz out_xyz_rc=out_xyz_rc.reshape(-1,2) blocksize=2000 ori_data=gdal.Open(ori_resampling_file_path[0][0][0]) width=ori_data.RasterXSize height=ori_data.RasterYSize SAR_ori=np.ones((width,height)) sim_x=ori_data.GetRasterBand(1).ReadAsArray(0,0 , width, height).reshape(-1) tree=KDTree(out_xyz_rc) for i in range(0,height,blocksize): for j in range(0,width,blocksize): start_row=i start_col=j len_row=blocksize len_col=blocksize if(start_row+len_row>height): len_row=height-start_row if(start_col+len_col>width): len_col=width-start_col end_row=start_row+len_row end_col=start_col+len_col grid_xy= np.mgrid[start_row:end_row:1, start_col:end_col:1] grid_xy=grid_xy.reshape(2,-1).transpose(1,0) distance_weight,ori_select=tree.query(grid_xy,k=4) distance_weight=np.reciprocal(distance_weight) distance_weight=np.multiply(distance_weight.transpose(1,0),np.reciprocal(np.sum(distance_weight,axis=1))) distance_weight=distance_weight.transpose(1,0) SAR_ori[start_row:end_row,start_col:end_col]=np.sum( np.multiply(distance_weight,sim_x[ori_select]),axis=1).reshape(len_row,len_col) pass ''' # 获得所有点集 del point_arr # 构建查询表 ori_sar=gdal.Open(ori_sar_path) ori_height=ori_sar.RasterYSize # 获得高度 ori_width=ori_sar.RasterXSize # 获得宽度 point_arr_path=[os.path.join(work_space_file,"point_refrence.npy")] for point_match_arr_filename in os.listdir(os.path.join(work_space_file,'temp_point_refrence')): point_arr_path.append(os.path.join(work_space_file,'temp_point_refrence',point_match_arr_filename)) point_arr=np.ones((1,5)) for point_path in point_arr_path: point_temp_arr=np.load(point_path) point_arr=np.row_stack((point_arr,point_temp_arr)) point_arr=point_arr[1:,:] ori_points=np.zeros((point_arr.shape[0],2)) # 求解对应的行列号 for i in range(r_arr_trans.shape[1]): ori_points[:,0]=ori_points[:,0]+r_arr_trans[0,i]*point_arr[:,1]**i ori_points[:,1]=ori_points[:,1]+c_arr_trans[0,i]*point_arr[:,1]**i point_arr[:,0]=ori_points[:,0].copy() point_arr[:,1]=ori_points[:,1].copy() del ori_points gc.collect() # 根据此进行散点插值，确定 lon_lat_point_match_path=os.path.join(work_space_file,"lon_lat_point_match") if os.path.exists(lon_lat_point_match_path): shutil.rmtree(lon_lat_point_match_path) os.makedirs(lon_lat_point_match_path) lon_lat=np.zeros((4,2)) block_info_json={"ori_rc":[]} # block_size=5000 for i in range(0,ori_height,block_size): for j in range(0,ori_width,block_size): len_row=block_size len_col=block_size if i+len_row>ori_height: len_row=ori_height-i if j+len_col>ori_width: len_col=ori_width-j end_i=i+len_row end_j=j+len_col r=(np.ones((len_col,len_row))*range(len_row)).transpose(1,0)+i c=np.ones((len_row,len_col))*range(len_col)+j # 格网插值 # scipy 规则格网插值 lon_lat_r_c=np.ones((len_col*len_row,5)) lon_lat_r_c[:,0]=r.reshape(-1) lon_lat_r_c[:,1]=c.reshape(-1) SAR_ori_X=interpolate.griddata(point_arr[:,:2], point_arr[:,2], (lon_lat_r_c[:,0], lon_lat_r_c[:,1]), method='linear').reshape(-1) SAR_ori_Y=interpolate.griddata(point_arr[:,:2], point_arr[:,3], (lon_lat_r_c[:,0], lon_lat_r_c[:,1]), method='linear').reshape(-1) SAR_ori_Z=interpolate.griddata(point_arr[:,:2], point_arr[:,4], (lon_lat_r_c[:,0], lon_lat_r_c[:,1]), method='linear').reshape(-1) [lon_lat_r_c[:,2],lon_lat_r_c[:,3],lon_lat_r_c[:,4]]=OrthoAuxData.XYZ2LLAM(SAR_ori_X,SAR_ori_Y,SAR_ori_Z) temp_path=os.path.join(lon_lat_point_match_path,"r_c_{}_{}.npy".format(i,j)) np.save(temp_path,lon_lat_r_c) block_info_json['ori_rc'].append({'ori':[i,j,end_i,end_j],"path":temp_path}) if i==0 and j==0: lon_lat[0,0]=lon_lat_r_c[:,2].reshape(len_row,len_col)[0,0] # 经度 lon_lat[0,1]=lon_lat_r_c[:,3].reshape(len_row,len_col)[0,0] # 纬度 if i+len_row==ori_height and j==0: lon_lat[1,0]=lon_lat_r_c[:,2].reshape(len_row,len_col)[-1,0] lon_lat[1,1]=lon_lat_r_c[:,3].reshape(len_row,len_col)[-1,0] if j+len_col==ori_width and i==0: lon_lat[2,0]=lon_lat_r_c[:,2].reshape(len_row,len_col)[0,-1] lon_lat[2,1]=lon_lat_r_c[:,3].reshape(len_row,len_col)[0,-1] if j+len_col==ori_width and i+len_row==ori_height: lon_lat[3,0]=lon_lat_r_c[:,2].reshape(len_row,len_col)[-1,-1] lon_lat[3,1]=lon_lat_r_c[:,3].reshape(len_row,len_col)[-1,-1] # 1 3 # 2 4 gt=[0,1,2,3,4,5] # 计算 经度 gt[0]= lon_lat[0,0] gt[1]=(lon_lat[2,0]-gt[0])/(ori_width-1) gt[2]=(lon_lat[3,0]-gt[0]-gt[1]*(ori_width-1))/(ori_height-1) # 计算 纬度 gt[3]= lon_lat[0,1] gt[5]=(lon_lat[1,1]-gt[3])/(ori_height-1) gt[4]=(lon_lat[3,1]-gt[3]-gt[5]*(ori_height-1))/(ori_width-1) return lon_lat_point_match_path,gt #return lon_lat_point_match_path return None pass '''