#!/usr/bin/env python3 import os,sys,struct,math import numpy as np import argparse from pyint import _utils as ut import subprocess #if len(sys.argv) != 5: # print("Usage:") # print("./applygacos.py inpfilename ztd1filename ztd2filename elevfilename") # print(" inpfilenmae : input interferogram, a inpfilename.rsc is needed!") # print(" ztd1filename: input GACOS ztd day1 filename, a ztd1filename.rsc is needed!") # print(" ztd2filename: input GACOS ztd day2 filename, a ztd2filename.rsc is needed!") # print(" elevfilename: elevation angle file, must be the same size with interferogram") # print(" elevation=90-incidence ") # print(" can be generated by look_vector in gamma") # print("outputfile will be saved as inpfilename.gacos") # print("!!Note, python3 is needed!!") # exit() def cmdLineParse(): parser = argparse.ArgumentParser(description='Coregister SM mode SLC to a reference SLC image using GAMMA.',\ formatter_class=argparse.RawTextHelpFormatter,\ epilog=INTRODUCTION+'\n'+EXAMPLE) parser.add_argument('projectName', help='Name of project.') parser.add_argument('mdate', help='date of the slave SLC image') parser.add_argument('sdate', help='date of the slave SLC image') parser.add_argument('ztd1', help='input GACOS ztd day1 filename, a ztd2filename.rsc is needed!') parser.add_argument('ztd2', help='input GACOS ztd day2 filename, a ztd2filename.rsc is needed!' ) inps = parser.parse_args() return inps INTRODUCTION = ''' ------------------------------------------------------------------- create_files for psokinv software running file by generated by Gamma ''' EXAMPLE = """Usage: coreg_gamma.py projectName coreg_gamma.py PacayaT163TsxHhA ------------------------------------------------------------------- """ #filename=sys.argv[1] #ztd1filename=sys.argv[2] #ztd2filename=sys.argv[3] #elevfilename=sys.argv[4] #print("start processing "+ filename) class HEADER: width = 0 length = 0 xfirst = 0.0 yfirst = 0.0 xstep = 0.0 ystep = 0.0 def read_header(filename): if not os.path.isfile(filename): print(filename+" file not exit") header = HEADER() with open(filename) as f: for line in f: data=line.split() if data[0] == "WIDTH": header.width=int(data[1]) if data[0] == "FILE_LENGTH": header.length=int(data[1]) if data[0] == "X_FIRST": header.xfirst=float(data[1]) if data[0] == "Y_FIRST": header.yfirst=float(data[1]) if data[0] == "X_STEP": header.xstep=float(data[1]) if data[0] == "Y_STEP": header.ystep=float(data[1]) return header def cut_image2(filename,headername,yfirst_new,length_new,xfirst_new,width_new): header=read_header(headername) #print(yfirst_new, header.yfirst) #print(length_new, header.length) #print(xfirst_new, header.xfirst) #print(width_new, header.width) with open(filename, 'rb') as f: data0 = np.fromfile(f, dtype=np.float32) data = np.reshape(data0, (header.length, header.width)) out=np.zeros((length_new,width_new),dtype=np.float32) #out=np.zeros(width_new*length_new,dtype=np.float32) for i in range(header.length): lat=header.yfirst+header.ystep*i row=int(round((lat-yfirst_new)/header.ystep)) if(row<0 or row>=length_new): continue for j in range(header.width): lon=header.xfirst+header.xstep*j col=int(round((lon-xfirst_new)/header.xstep)) if(col<0 or col>=width_new): continue #out[width_new*row+col]=data[header.width*i+j] out[row,col]=data[i,j] out=np.where(out==0,np.nan,out) #print(np.nanstd(out)) out.tofile(filename+".cut") f = open(filename+".cut.rsc", 'w') f.write("WIDTH "+str(width_new)+"\n") f.write("FILE_LENGTH "+str(length_new)+"\n") f.write("X_FIRST "+ str(xfirst_new)+"\n") f.write("Y_FIRST "+ str(yfirst_new)+"\n") f.write("X_STEP "+ str(header.xstep)+"\n") f.write("Y_STEP "+ str(header.ystep)+"\n") f.write("X_UNIT degree"+"\n") f.write("Y_UNIT degree"+"\n") f.write("Z_OFFSET 0"+"\n") f.write("Z_SCALE 1"+"\n") f.write("PROJECTION LATLON"+"\n") f.write("DATUM WGS84"+"\n") f.close() def make_correction(phsfilename,ztd1filename,ztd2filename,elevfilename): header=read_header(phsfilename+".rsc") #if not os.path.isfile(ztd1filename+".cut"): cut_image2(ztd1filename,ztd1filename+".rsc",header.yfirst,header.length,header.xfirst,header.width) #if not os.path.isfile(ztd2filename+".cut"): cut_image2(ztd2filename,ztd2filename+".rsc",header.yfirst,header.length,header.xfirst,header.width) with open(phsfilename, 'rb') as f: data = np.fromfile(f, dtype=np.float32) phase = np.reshape(data, [header.length, header.width]) with open(ztd1filename+".cut", 'rb') as f: data = np.fromfile(f, dtype=np.float32) ztd1 = np.reshape(data, [header.length, header.width]) with open(ztd2filename+".cut", 'rb') as f: data = np.fromfile(f, dtype=np.float32) ztd2 = np.reshape(data, [header.length, header.width]) with open(elevfilename, 'rb') as f: data = np.fromfile(f, dtype=np.float32) elev = np.reshape(data, [header.length, header.width]) #os.remove(ztd1filename+".cut") #os.remove(ztd2filename+".cut") dztd=ztd2-ztd1 #dztd=dztd/0.0044138251819503 dztd=dztd/np.sin(elev) index=np.where(phase==0) phase[index]=np.nan phasemean=np.nanmean(phase) print("before "+str(np.nanstd(phase))) phase=phase-phasemean #phase.tofile(phsfilename+".raw") phase=phase-dztd phase[index]=np.nan phasemean=np.nanmean(phase) print("after "+str(np.nanstd(phase))) phase=phase-phasemean phase[index]=0 phase.tofile(phsfilename+".gacos") def main(argv): inps = cmdLineParse() projectName = inps.projectName Sdate = inps.sdate Mdate = inps.mdate scratchDir = os.getenv('SCRATCHDIR') templateDir = os.getenv('TEMPLATEDIR') templateFile = templateDir + "/" + projectName + ".template" templateDict=ut.update_template(templateFile) rlks = templateDict['range_looks'] azlks = templateDict['azimuth_looks'] Smdate = templateDict['masterDate'] demDir = scratchDir + '/' + projectName + '/DEM' slcDir = scratchDir + '/' + projectName + "/SLC" rslcDir = scratchDir + '/' + projectName + "/RSLC" ifgramDir = scratchDir + '/' + projectName + "/ifgrams" workDir = ifgramDir + '/' + Mdate + '-' + Sdate MslcDir = rslcDir + '/' + Mdate SslcDir = rslcDir + '/' + Sdate offpar = SslcDir + '/' + Smdate + '_' + Sdate + '.off' slcpar = MslcDir + '/' + Smdate + '.rslc.par' dempar = demDir + '/' + Smdate + '_' + rlks + 'rlks.utm.dem.par' dem = demDir + '/' + Smdate + '_' + rlks + 'rlks.utm.dem' unwpar = demDir + '/' + Mdate + '_' + rlks + 'rlks.utm.dem.par' unw = workDir + '/geo_' + Mdate + '-' + Sdate + "_" + rlks+"rlks.diff_filt.unw" os.chdir(workDir) call_str = "look_vector " + slcpar+ " " + offpar + " " + dempar + " " + dem + ' lv_theta lv_phi' os.system(call_str) call_str = "grep 'corner_lat:' " + dempar + " | awk '{print $2}' " North=subprocess.getstatusoutput(call_str)[1] call_str = "grep 'corner_lon:' " + dempar + " | awk '{print $2}' " West=subprocess.getstatusoutput(call_str)[1] call_str = "grep 'post_lat:' " + dempar + " | awk '{print $2}' " posty=subprocess.getstatusoutput(call_str)[1] call_str = "grep 'post_lon:' " + dempar + " | awk '{print $2}' " postx=subprocess.getstatusoutput(call_str)[1] call_str = "grep 'width:' " + dempar + " | awk '{print $2}' " width=subprocess.getstatusoutput(call_str)[1] call_str = "grep 'nlines:' " + dempar + " | awk '{print $2}' " length=subprocess.getstatusoutput(call_str)[1] South=str(round(float(North)+(float(length)-1)*float(posty),7)) East=str(round(float(West)+(float(width)-1)*float(postx),7)) call_str = "swap_bytes lv_theta lv_theta.phase_swap 4 >dinsar.log " os.system(call_str) call_str = "xyz2grd lv_theta.phase_swap -Glv_theta.grd -Ddegree/degree/cm/1/0/=/= -R" + West + "/" + East + "/" +South + "/" +North + " -I" + postx + " -ZTLf -N0" os.system(call_str) call_str = "grdmath 90 lv_theta.grd 3.1415926 DIV 180 MUL SUB = lv_theta_final.grd" os.system(call_str) call_str = "grdmath 90 lv_theta_final.grd SUB = lv_elev.grd" os.system(call_str) call_str = "grdsample lv_elev.grd -Glv_elev_3c.grd -I3c" # call_str = "grdsample lv_theta_final.grd -Glv_elev_3c.grd -I3c" os.system(call_str) call_str = "grd2xyz lv_elev_3c.grd -ZTLf -N0 > " + Mdate + '-' + Sdate + '.gacos.elev' os.system(call_str) call_str = "swap_bytes " + unw+ " unw.phase_swap 4 >dinsar.log " os.system(call_str) call_str = "xyz2grd unw.phase_swap -Gunw_f.grd -Ddegree/degree/cm/1/0/=/= -R" + West + "/" + East + "/" +South + "/" +North + " -I" + postx +" -ZTLf -N0" os.system(call_str) call_str = "grdsample unw_f.grd -Gunw.grd -I3c" os.system(call_str) call_str = "grd2xyz unw.grd -ZTLf -N0 > " + Mdate + '-' + Sdate + '.gacos.unw' os.system(call_str) call_str = "grdinfo unw.grd -C | awk '{print $10}' " Width=subprocess.getstatusoutput(call_str)[1] call_str = "grdinfo unw.grd -C | awk '{print $11}' " line=subprocess.getstatusoutput(call_str)[1] call_str = "grdinfo unw.grd -C | awk '{print $2}' " West=subprocess.getstatusoutput(call_str)[1] call_str = "grdinfo unw.grd -C | awk '{print $5}' " North=subprocess.getstatusoutput(call_str)[1] ymax=str(int(round(float(line) ,7))) xmax=str(int(round(float(Width) ,7))) output=workDir + '/' + Mdate + '-' + Sdate + '.gacos.unw.rsc' if os.path.exists(output) is True: os.remove(output) fopen=open(output,'a+') fopen.write('WIDTH ' + Width + '\n') fopen.write('FILE_LENGTH ' + line + '\n') fopen.write('XMIN 1' + '\n') fopen.write('XMAX ' + xmax + '\n') fopen.write('YMIN 1' + '\n') fopen.write('YMAX ' + ymax + '\n') fopen.write('X_FIRST ' + West + '\n') fopen.write('Y_FIRST ' + North + '\n') fopen.write('X_STEP 8.33333333E-04' + '\n') fopen.write('Y_STEP -8.33333333E-04' + '\n') fopen.write('X_UNIT degrees' + '\n') fopen.write('Y_UNIT degrees' + '\n') fopen.write('Z_OFFSET 0' + '\n') fopen.write('Z_SCALE 1' + '\n') fopen.write('PROJECTION LATLON' + '\n') fopen.write('DATUM WGS84' + '\n') fopen.close() print("convert gamma products to GACOS is done!") print("start make gacos correction!") filename=workDir + '/' + Mdate + '-' + Sdate + '.gacos.unw' ztd1=inps.ztd1 ztd2=inps.ztd2 ztd1filename=workDir + '/' + ztd1 ztd2filename=workDir + '/' + ztd2 elevfilename=workDir + '/' + Mdate + '-' + Sdate + '.gacos.elev' make_correction(filename,ztd1filename,ztd2filename,elevfilename) sys.exit(1) if __name__ == '__main__': main(sys.argv[:]) #make_correction(filename,ztd1filename,ztd2filename,elevfilename) #exit()