416 lines
14 KiB
Python
416 lines
14 KiB
Python
# 一米正射辅助数据处理类
|
||
import time
|
||
import math
|
||
import numpy as np
|
||
from osgeo import gdal
|
||
from xml.etree.ElementTree import ElementTree
|
||
from scipy.optimize import leastsq
|
||
|
||
|
||
class OrthoAuxData:
|
||
def __init__(self):
|
||
pass
|
||
|
||
@staticmethod
|
||
def time_stamp(tm):
|
||
list = tm.split(':')
|
||
sec = math.ceil(float(list[2]))
|
||
tm1 = list[0] + ':' + list[1] + ':' + str(sec)
|
||
tmArr = time.strptime(tm1, "%Y-%m-%d %H:%M:%S")
|
||
# tmArr = time.strptime(tm1, "%Y-%m-%d %H:%M:%S.%f")
|
||
ts = float(time.mktime(tmArr)) # 转换为时间戳
|
||
return ts
|
||
|
||
@staticmethod
|
||
def read_meta(meta_file_path):
|
||
tree = ElementTree()
|
||
tree.parse(meta_file_path)
|
||
root = tree.getroot()
|
||
T = []
|
||
Xs = []
|
||
Ys = []
|
||
Zs = []
|
||
Vsx = []
|
||
Vsy = []
|
||
Vsz = []
|
||
GPS_data = root.find('GPS')
|
||
for child in GPS_data:
|
||
Xs.append(float(child.find('xPosition').text))
|
||
Ys.append(float(child.find('yPosition').text))
|
||
Zs.append(float(child.find('zPosition').text))
|
||
Vsx.append(float(child.find('xVelocity').text))
|
||
Vsy.append(float(child.find('yVelocity').text))
|
||
Vsz.append(float(child.find('zVelocity').text))
|
||
tm = child.find('TimeStamp').text
|
||
ts = OrthoAuxData.time_stamp(tm)
|
||
T.append(ts)
|
||
meta_data = [Xs, Ys, Zs, Vsx, Vsy, Vsz]
|
||
return T, meta_data
|
||
|
||
@staticmethod
|
||
def read_control_points(meta_file_path):
|
||
tree = ElementTree()
|
||
tree.parse(meta_file_path)
|
||
root = tree.getroot()
|
||
imageinfo = root.find('imageinfo')
|
||
center = imageinfo.find('center')
|
||
corner = imageinfo.find('corner')
|
||
ctrl_pts = [[] for i in range(2)]
|
||
ctrl_pts[0].append(float(center.find('longitude').text))
|
||
ctrl_pts[1].append(float(center.find('latitude').text))
|
||
|
||
for child in corner:
|
||
ctrl_pts[0].append(float(child.find('longitude').text))
|
||
ctrl_pts[1].append(float(child.find('latitude').text))
|
||
return ctrl_pts
|
||
|
||
@staticmethod
|
||
def read_dem(dem_resampled_path, flag=1):
|
||
in_ds = gdal.Open(dem_resampled_path)
|
||
gt = list(in_ds.GetGeoTransform())
|
||
bands_num = in_ds.RasterCount
|
||
x_size = in_ds.RasterXSize
|
||
y_size = in_ds.RasterYSize
|
||
pstn_arr = np.zeros([y_size, x_size, 3], dtype=np.float)
|
||
for i in range(1, bands_num + 1):
|
||
data = in_ds.GetRasterBand(i).ReadAsArray(0, 0, x_size, y_size)
|
||
for y in range(y_size):
|
||
for x in range(x_size):
|
||
longitude = gt[0] + x * gt[1]
|
||
latitude = gt[3] + y * gt[5]
|
||
altitude = data[y, x]
|
||
if flag == 1:
|
||
pstn = OrthoAuxData.LLA2XYZ(longitude, latitude, altitude)
|
||
else:
|
||
pstn = [longitude, latitude, altitude]
|
||
pstn_arr[y, x, 0] = pstn[0]
|
||
pstn_arr[y, x, 1] = pstn[1]
|
||
pstn_arr[y, x, 2] = pstn[2]
|
||
|
||
del in_ds, data
|
||
return pstn_arr
|
||
|
||
@staticmethod
|
||
def read_demM(dem_resampled_path, part_cnt, r_cnt, c_cnt, flag=1):
|
||
in_ds = gdal.Open(dem_resampled_path)
|
||
gt = list(in_ds.GetGeoTransform())
|
||
bands_num = in_ds.RasterCount
|
||
x_size = in_ds.RasterXSize // part_cnt
|
||
y_size = in_ds.RasterYSize // part_cnt
|
||
x = [[i] * y_size for i in range(x_size)]
|
||
y = [[i] * x_size for i in range(y_size)]
|
||
x = np.array(x)
|
||
x = x.T
|
||
y = np.array(y)
|
||
x_off = c_cnt * x_size
|
||
y_off = r_cnt * y_size
|
||
gt[0] = gt[0] + c_cnt * x_size * gt[1]
|
||
gt[3] = gt[3] + r_cnt * y_size * gt[5]
|
||
for i in range(1, bands_num + 1):
|
||
data = in_ds.GetRasterBand(i).ReadAsArray(x_off, y_off, x_size, y_size)
|
||
altitude = data / 255 * 1024
|
||
longitude = gt[0] + x * gt[1]
|
||
latitude = gt[3] + y * gt[5]
|
||
if flag == 1:
|
||
pstn = OrthoAuxData.LLA2XYZM(longitude, latitude, altitude)
|
||
else:
|
||
pstn = [longitude, latitude, altitude]
|
||
|
||
del in_ds, data
|
||
return pstn
|
||
|
||
@staticmethod
|
||
def read_dem_row(dem_resampled_path, p, flag=1):
|
||
in_ds = gdal.Open(dem_resampled_path)
|
||
gt = list(in_ds.GetGeoTransform())
|
||
bands_num = in_ds.RasterCount
|
||
x_size = in_ds.RasterXSize
|
||
y_size = in_ds.RasterYSize
|
||
|
||
x = [[i] for i in range(x_size)]
|
||
x = np.array(x)
|
||
x = x.T
|
||
y = np.ones((1, x_size)) * p
|
||
x_off = 0
|
||
y_off = p
|
||
for i in range(1, bands_num + 1):
|
||
data = in_ds.GetRasterBand(i).ReadAsArray(x_off, y_off, x_size, 1)
|
||
altitude = data
|
||
longitude = gt[0] + x * gt[1]
|
||
latitude = gt[3] + y * gt[5]
|
||
if flag == 1:
|
||
pstn = OrthoAuxData.LLA2XYZM(longitude, latitude, altitude)
|
||
else:
|
||
pstn = [longitude, latitude, altitude]
|
||
|
||
del in_ds, data
|
||
return pstn
|
||
|
||
@staticmethod
|
||
def orbit_fitting(time_array, meta_data):
|
||
# 最小二乘法求解轨道参数
|
||
T0 = (time_array[0] + time_array[len(time_array)-1]) / 2
|
||
t = []
|
||
for i in range(len(time_array)):
|
||
t.append(time_array[i]-T0)
|
||
|
||
def func(p, x):
|
||
w3, w2, w1, w0 = p
|
||
return w3*x**3 + w2*x**2 + w1*x + w0
|
||
|
||
def error(p, x, y):
|
||
return func(p, x) - y
|
||
|
||
orbital_paras = []
|
||
for j in range(len(meta_data)):
|
||
p0 = [1, 2, 3, 4]
|
||
x = np.array(t)
|
||
y = np.array(meta_data[j])
|
||
Para = leastsq(error, p0, args=(x, y))
|
||
orbital_paras.append(Para[0])
|
||
print(Para[0], Para[1])
|
||
|
||
return orbital_paras, T0
|
||
|
||
@staticmethod
|
||
def get_PRF(meta_file_path):
|
||
tree = ElementTree()
|
||
tree.parse(meta_file_path)
|
||
root = tree.getroot()
|
||
sensor = root.find('sensor')
|
||
waveParams = sensor.find('waveParams')
|
||
PRF = float(waveParams.find('wave').find('prf').text)
|
||
return PRF
|
||
|
||
@staticmethod
|
||
def get_delta_R(meta_file_path):
|
||
tree = ElementTree()
|
||
tree.parse(meta_file_path)
|
||
root = tree.getroot()
|
||
sensor = root.find('sensor')
|
||
pulseWidth = float(sensor.find('waveParams').find('wave').find('pulseWidth').text)
|
||
bandWidth = float(sensor.find('waveParams').find('wave').find('bandWidth').text)
|
||
c = 299792458
|
||
delta_R = c / (1000000 * 2 * bandWidth)
|
||
return delta_R
|
||
|
||
@staticmethod
|
||
def get_doppler_rate_coef(meta_file_path):
|
||
tree = ElementTree()
|
||
tree.parse(meta_file_path)
|
||
root = tree.getroot()
|
||
processinfo = root.find('processinfo')
|
||
doppler = processinfo.find('DopplerRateValuesCoefficients')
|
||
t0 = float(processinfo.find('DopplerParametersReferenceTime').text)
|
||
r0 = float(doppler.find('r0').text)
|
||
r1 = float(doppler.find('r1').text)
|
||
r2 = float(doppler.find('r2').text)
|
||
r3 = float(doppler.find('r3').text)
|
||
r4 = float(doppler.find('r4').text)
|
||
|
||
return t0, np.array([r0, r1, r2, r3, r4]).reshape(5, 1)
|
||
|
||
@staticmethod
|
||
def get_doppler_center_coef(meta_file_path):
|
||
tree = ElementTree()
|
||
tree.parse(meta_file_path)
|
||
root = tree.getroot()
|
||
processinfo = root.find('processinfo')
|
||
doppler = processinfo.find('DopplerCentroidCoefficients')
|
||
b0 = float(doppler.find('d0').text)
|
||
b1 = float(doppler.find('d1').text)
|
||
b2 = float(doppler.find('d2').text)
|
||
return b0, b1, b2
|
||
|
||
@staticmethod
|
||
def get_lamda(meta_file_path):
|
||
tree = ElementTree()
|
||
tree.parse(meta_file_path)
|
||
root = tree.getroot()
|
||
sensor = root.find('sensor')
|
||
λ = float(sensor.find('lamda').text)
|
||
return λ
|
||
|
||
@staticmethod
|
||
def get_t0(meta_file_path):
|
||
tree = ElementTree()
|
||
tree.parse(meta_file_path)
|
||
root = tree.getroot()
|
||
imageinfo = root.find('imageinfo')
|
||
tm = imageinfo.find('imagingTime').find('start').text
|
||
t0 = OrthoAuxData.time_stamp(tm)
|
||
return t0
|
||
|
||
@staticmethod
|
||
def get_start_and_end_time(meta_file_path):
|
||
tree = ElementTree()
|
||
tree.parse(meta_file_path)
|
||
root = tree.getroot()
|
||
imageinfo = root.find('imageinfo')
|
||
tm0 = imageinfo.find('imagingTime').find('start').text
|
||
tm1 = imageinfo.find('imagingTime').find('end').text
|
||
starttime = OrthoAuxData.time_stamp(tm0)
|
||
endtime = OrthoAuxData.time_stamp(tm1)
|
||
return starttime, endtime
|
||
|
||
@staticmethod
|
||
def get_width_and_height(meta_file_path):
|
||
tree = ElementTree()
|
||
tree.parse(meta_file_path)
|
||
root = tree.getroot()
|
||
imageinfo = root.find('imageinfo')
|
||
width = int(imageinfo.find('width').text)
|
||
height = int(imageinfo.find('height').text)
|
||
return width, height
|
||
|
||
@staticmethod
|
||
def get_R0(meta_file_path):
|
||
tree = ElementTree()
|
||
tree.parse(meta_file_path)
|
||
root = tree.getroot()
|
||
imageinfo = root.find('imageinfo')
|
||
R0 = float(imageinfo.find('nearRange').text)
|
||
return R0
|
||
|
||
@staticmethod
|
||
def get_h():
|
||
h = 6.6
|
||
return h
|
||
|
||
@staticmethod
|
||
def LLA2XYZ(longitude, latitude, altitude):
|
||
'''
|
||
WGS-84坐标系下：经纬度坐标转空间直角坐标
|
||
'''
|
||
# 经纬度余弦值
|
||
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 #地球赤道平均半径
|
||
f = 1.0/298.257224 #WGS84椭球扁率
|
||
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)
|
||
return [X, Y, Z]
|
||
|
||
@staticmethod
|
||
def LLA2XYZM(longitude, latitude, altitude):
|
||
# 经纬度余弦值
|
||
cosLat = np.cos(latitude * math.pi / 180).reshape(-1,1)
|
||
sinLat = np.sin(latitude * math.pi / 180).reshape(-1,1)
|
||
cosLon = np.cos(longitude * math.pi / 180).reshape(-1,1)
|
||
sinLon = np.sin(longitude * math.pi / 180).reshape(-1,1)
|
||
# WGS84坐标系参数
|
||
rad = 6378137.0 #地球赤道平均半径
|
||
f = 1.0/298.257224 #WGS84椭球扁率
|
||
C = 1.0/(np.sqrt(cosLat*cosLat + (1-f)*(1-f)*sinLat*sinLat)).reshape(-1,1)
|
||
S = (1-f)*(1-f)*C
|
||
h = altitude.reshape(-1,1)
|
||
# 计算XYZ坐标
|
||
X = (rad * C + h) * cosLat * cosLon
|
||
Y = (rad * C + h) * cosLat * sinLon
|
||
Z = (rad * S + h) * sinLat
|
||
return [X, Y, Z]
|
||
|
||
@staticmethod
|
||
def XYZ2LLA(X, Y, Z):
|
||
''' 大地坐标系转经纬度
|
||
适用于WGS84坐标系
|
||
args:
|
||
x,y,z
|
||
return:
|
||
lat,long,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])
|
||
return [np.degrees(longitude), np.degrees(latitude), altitude]
|
||
|
||
@staticmethod
|
||
def XYZ2LLAM(X, Y, Z):
|
||
''' 大地坐标系转经纬度
|
||
适用于WGS84坐标系
|
||
args:
|
||
x,y,z
|
||
return:
|
||
lat,long,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])
|
||
return [np.degrees(longitude), np.degrees(latitude), altitude]
|
||
|
||
@staticmethod
|
||
def world2Pixel(geoMatrix, x, y):
|
||
"""
|
||
使用GDAL库的geomatrix对象((gdal.GetGeoTransform()))计算地理坐标的像素位置
|
||
"""
|
||
ulx = geoMatrix[0]
|
||
uly = geoMatrix[3]
|
||
xDist = geoMatrix[1]
|
||
yDist = geoMatrix[5]
|
||
rtnX = geoMatrix[2]
|
||
rtnY = geoMatrix[4]
|
||
pixel = int((x - ulx) / xDist)
|
||
line = int((uly - y) / abs(yDist))
|
||
return pixel, line
|
||
|
||
@staticmethod
|
||
def sar_intensity_synthesis(in_sar_tif, out_sar_tif):
|
||
# 获取SLC格式SAR影像的相关信息
|
||
in_ds = gdal.Open(in_sar_tif)
|
||
bands_num = in_ds.RasterCount
|
||
rows = in_ds.RasterYSize
|
||
columns = in_ds.RasterXSize
|
||
proj = in_ds.GetProjection()
|
||
geotrans = in_ds.GetGeoTransform()
|
||
|
||
# 创建输出的SAR强度图
|
||
gtiff_driver = gdal.GetDriverByName('GTiff')
|
||
out_ds = gtiff_driver.Create(out_sar_tif, columns, rows, 1)
|
||
out_ds.SetProjection(proj)
|
||
out_ds.SetGeoTransform(geotrans)
|
||
|
||
# 输出SAR强度图
|
||
in_data1 = in_ds.GetRasterBand(1).ReadAsArray(0, 0, columns, rows)
|
||
in_data1 = in_data1/10
|
||
in_data1 = np.power(10, in_data1)
|
||
in_data2 = in_ds.GetRasterBand(2).ReadAsArray(0, 0, columns, rows)
|
||
in_data2 = in_data2 / 10
|
||
in_data2 = np.power(10, in_data2)
|
||
out_data = np.sqrt(in_data1**2 + in_data2**2)
|
||
out_ds.GetRasterBand(1).WriteArray(out_data)
|
||
|
||
del in_ds, out_ds
|