# -*- coding: utf-8 -*-
"""
Created on Mon Sep 14 14:10:47 2020
去地平
@author: yang
"""
import matplotlib.pyplot as plt #nice 2D plots
import numpy as np #matrix/array calculations etc
from matplotlib import cm #colormap shortcut
from phase_unwrap_1D import phase_unwrap_1D
####################
### data import ####
####################
def unwrap_and_residues(phase_sea,c,fc,H,Fr,daitah,theta_B,R0,B,Naz,Offset):
## 参数计算
Nrg=phase_sea.shape[0]
lamda = c/fc
#########################
### unwrap one line ####
#########################
# selection of one line
rw, cl = phase_sea.shape
cl_vec = np.array(range(0,cl,1))
rw_selected =np.int(rw/2)# selected line number (row)
# unwrapping one line
phase_1D = phase_sea[rw_selected,]
uw_phase_1D_dict = phase_unwrap_1D(0, phase_1D)
# get data from dictionary
uw_phase_1D = uw_phase_1D_dict['uw_phase']
uw_cycles_1D = uw_phase_1D_dict['uw_cycles']
# show wrapped and unwrapped phases
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
cax1 = ax.plot(cl_vec, phase_1D, 'k.', label='initial phase')
cax2 = ax.plot(cl_vec, uw_phase_1D, 'r.', label='unwrapped phase')
leg = plt.legend(fancybox=True, loc='best')#设置图例参数fancybox: 是否将图例框的边角设为圆形，loc:图例位置
ax.axhline(y=-np.pi, xmin=0, xmax=cl-1, color='b')#绘制平行于x轴的水平参考线
ax.axhline(y=+np.pi, xmin=0, xmax=cl-1, color='b')
ax.set_title('Unwrapping phase of one line')
plt.axis('tight')
plt.show()
# show_cycle plot
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
cax = ax.plot(cl_vec, uw_cycles_1D)
ax.set_title('Total number of cycles used for unwrapping')
plt.axis('tight')
plt.show()
rw_all, cl_all = np.shape(phase_sea);
###########
##二维相位解缠
###########
# phase unwrapping line wise
# 1. do phase unwrapping for first column to get initial values for each line
phase_1D = phase_sea[:,0]
uw_phase_1D_dict = phase_unwrap_1D(0, phase_1D)
final_phase_1st_cl = uw_phase_1D_dict['uw_phase']
cycles_1st_cl = uw_phase_1D_dict['uw_cycles']
# 2. apply PhU for each line in interferogram, using initial values for each line
final_phase_2D_rw = np.zeros( (rw_all, cl_all) )
for k in range(0,rw_all):
phase_1D = phase_sea[k,:]
tmp_dict = phase_unwrap_1D( cycles_1st_cl[k], phase_1D )
final_phase_2D_rw[k,:] = tmp_dict['uw_phase']
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
cax = ax.imshow(final_phase_2D_rw, interpolation='nearest', cmap=cm.jet)
ax.set_title('Unwrapped phases (row wise)')
plt.axis('tight')
plt.show()
phase_1D = phase_sea[0,:]
uw_phase_1D_dict = phase_unwrap_1D(0, phase_1D)
final_phase_1st_rw = uw_phase_1D_dict['uw_phase']
cycles_1st_rw = uw_phase_1D_dict['uw_cycles']
# 2. apply PhU for each line in interferogram, using initial values for each line
final_phase_2D_cl = np.zeros( (rw_all, cl_all) )
for k in range(0,cl_all):
phase_1D = phase_sea[:,k]
tmp_dict = phase_unwrap_1D( cycles_1st_rw[k], phase_1D )
final_phase_2D_cl[:,k] = tmp_dict['uw_phase']
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
cax = ax.imshow(final_phase_2D_cl, interpolation='nearest', cmap=cm.jet)
ax.set_title('Unwrapped phases (column wise)')
plt.axis('tight')
plt.show()
#生成相对海平面高程
# tr_RCMC = np.zeros((1,Naz), dtype=np.float )
# for i in range(Naz):
# tr_RCMC[0,0]=2*R0/c-Naz/2/Fr
# tr_RCMC[0,i]=tr_RCMC[0,i-1]+Naz/2/Fr/Naz*2
# R0_RCMC = (c/2)*tr_RCMC #随距离线变化的R0，记为R0_RCMC，用来计算RCM
R0_RCMC = R0 + np.arange(Naz) / Fr *c/2
R0_RCMC_mtx = np.ones((Nrg,1))*R0_RCMC #形成斜距矩阵
#裁剪同一区域
if Offset>=0:
r_range_A = R0_RCMC_mtx[:,Offset:]
else :
r_range_A = R0_RCMC_mtx[:,:Offset] # 截取所计算区域
cos_theta=np.zeros((r_range_A.shape[0],r_range_A.shape[1]), dtype=np.float )
theta_range_A=np.zeros((r_range_A.shape[0],r_range_A.shape[1]), dtype=np.float )
cos_theta= (H-daitah)/r_range_A # 对应于每一个天线A的最近斜距时，波束下视角的余弦
theta_range_A=np.arccos(cos_theta) #计算出波束下视角（弧度）；
#计算相对海平面高度
h=np.zeros((phase_sea.shape[0],phase_sea.shape[1]), dtype=np.float )
ad=-1*lamda*r_range_A/2/np.pi/B
print(theta_range_A.shape[1])
print(final_phase_2D_cl.shape[1])
h=final_phase_2D_cl*ad*np.sin(theta_range_A)/np.cos(theta_range_A-theta_B)
h1=h[20:-20,20:-40]
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
cax = ax.imshow(h1, aspect='equal',interpolation='nearest', cmap='ocean')
cbar = fig.colorbar(cax)
ax.set_title('Height map [m]')
plt.axis('tight')
plt.show()
return h
评论5