diff --git a/README.md b/README.md index 412a0ac..a8787ba 100644 --- a/README.md +++ b/README.md @@ -152,6 +152,7 @@ py36-scipy +gcc7 +openblas py36-matplotlib +cairo +tkinter py36-matplotlib-basemap py36-h5py +py36-gdal ``` ### Python3 Convention diff --git a/SConstruct b/SConstruct index 973f7f2..a979b3c 100644 --- a/SConstruct +++ b/SConstruct @@ -216,10 +216,11 @@ else: ### End of GPU branch-specific modifications -os.makedirs(inst, exist_ok=True) env.Install(inst, '__init__.py') env.Install(inst, 'release_history.py') +if not os.path.exists(inst): + os.makedirs(inst) v = 0 if isrerun == 'no': diff --git a/applications/dem.py b/applications/dem.py index a4f5fc6..9cb307f 100755 --- a/applications/dem.py +++ b/applications/dem.py @@ -82,7 +82,7 @@ def main(): parser.add_argument('-n', '--uname', type = str, dest = 'uname', default = None, help = 'User name if using a server that requires authentication') parser.add_argument('-w', '--password', type = str, dest = 'password', default = None, help = 'Password if using a server that requires authentication') parser.add_argument('-t', '--type', type = str, dest = 'type', default = 'version3', help = \ - 'Use version 3 or version 2 SRTM') + 'Use version 3 or version 2 SRTM, or nasadem') parser.add_argument('-x', '--noextras', action = 'store_true', dest = 'noextras', help = 'Use this flag if the filenames do not have extra part') parser.add_argument('-u', '--url', type = str, dest = 'url', default = None, help = \ 'If --type=version2 then this is part of the url where the DEM files are located. The actual location must be' + \ @@ -91,10 +91,13 @@ def main(): args = parser.parse_args() #first get the url,uname and password since are needed in the constructor - ds = createDemStitcher(args.type) ds.configure() + #NASADEM only available in 1-arc sec resolution + if(args.type == 'nasadem'): + args.source == 1 + if(args.url): if(args.type == 'version3'): if(args.source == 1): diff --git a/applications/gdal2isce_xml.py b/applications/gdal2isce_xml.py index 9477c01..3eab342 100755 --- a/applications/gdal2isce_xml.py +++ b/applications/gdal2isce_xml.py @@ -3,12 +3,14 @@ # Author: Bekaert David # Year: 2017 +import os +import sys +import argparse +from osgeo import gdal + import isce import isceobj -import sys -from osgeo import gdal -import argparse -import os + # command line parsing of input file def cmdLineParse(): @@ -20,17 +22,14 @@ def cmdLineParse(): return parser.parse_args() -# main script -if __name__ == '__main__': - ''' - Main driver. - ''' - - # Parse command line - inps = cmdLineParse() - # check if the input file exist - if not os.path.isfile(inps.fname): - raise Exception('Input file is not found ....') +def gdal2isce_xml(fname): + """ + Generate ISCE xml file from gdal supported file + + Example: import isce + from applications.gdal2isce_xml import gdal2isce_xml + xml_file = gdal2isce_xml(fname+'.vrt') + """ # open the GDAL file and get typical data informationi GDAL2ISCE_DATATYPE = { @@ -57,47 +56,46 @@ if __name__ == '__main__': # } # check if the input file is a vrt - filename, file_extension = os.path.splitext(inps.fname) - print(file_extension) - if file_extension == ".vrt": - inps.outname = filename + fbase, fext = os.path.splitext(fname) + print(fext) + if fext == ".vrt": + outname = fbase else: - inps.outname = inps.fname - print(inps.outname) + outname = fname + print(outname) - # open the GDAL file and get typical data informationi - data = gdal.Open(inps.fname, gdal.GA_ReadOnly) - width = data.RasterXSize - length = data.RasterYSize - bands = data.RasterCount - # output to user - print("width: " + "\t" + str(width)) - print("length: " + "\t" + str(length)) - print("nof bands:" + "\t" + str(bands)) + # open the GDAL file and get typical ds information + ds = gdal.Open(fname, gdal.GA_ReadOnly) + width = ds.RasterXSize + length = ds.RasterYSize + bands = ds.RasterCount + print("width: " + "\t" + str(width)) + print("length: " + "\t" + str(length)) + print("num of bands:" + "\t" + str(bands)) # getting the datatype information - raster = data.GetRasterBand(1) + raster = ds.GetRasterBand(1) dataTypeGdal = raster.DataType + # user look-up dictionary from gdal to isce format dataType= GDAL2ISCE_DATATYPE[dataTypeGdal] - # output to user print("dataType: " + "\t" + str(dataType)) - # transformation contains gridcorners (lines/pixels or lonlat and the spacing 1/-1 or deltalon/deltalat) - transform = data.GetGeoTransform() + transform = ds.GetGeoTransform() # if a complex data type, then create complex image # if a real data type, then create a regular image img = isceobj.createImage() - img.setFilename(os.path.abspath(inps.outname)) + img.setFilename(os.path.abspath(outname)) img.setWidth(width) img.setLength(length) img.setAccessMode('READ') img.bands = bands img.dataType = dataType - - md = data.GetMetadata('IMAGE_STRUCTURE') + + # interleave + md = ds.GetMetadata('IMAGE_STRUCTURE') sch = md.get('INTERLEAVE', None) if sch == 'LINE': img.scheme = 'BIL' @@ -110,10 +108,29 @@ if __name__ == '__main__': print('Assuming default, BIP') img.scheme = 'BIP' - img.firstLongitude = transform[0] img.firstLatitude = transform[3] img.deltaLatitude = transform[5] - img.deltaLongitude = transform[1] - img.dump(inps.outname + ".xml") + img.deltaLongitude = transform[1] + + xml_file = outname + ".xml" + img.dump(xml_file) + + return xml_file + + +# main script +if __name__ == '__main__': + ''' + Main driver. + ''' + + # Parse command line + inps = cmdLineParse() + + # check if the input file exist + if not os.path.isfile(inps.fname): + raise Exception('Input file is not found ....') + + gdal2isce_xml(inps.fname) diff --git a/components/isceobj/Alos2Proc/runDownloadDem.py b/components/isceobj/Alos2Proc/runDownloadDem.py index 36b9ae6..40d2d4e 100644 --- a/components/isceobj/Alos2Proc/runDownloadDem.py +++ b/components/isceobj/Alos2Proc/runDownloadDem.py @@ -34,16 +34,24 @@ def runDownloadDem(self): os.makedirs(demDir, exist_ok=True) os.chdir(demDir) - downloadUrl = 'http://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11' - cmd = 'dem.py -a stitch -b {} -k -s 1 -c -f -u {}'.format( - bboxStr, - downloadUrl - ) - runCmd(cmd) - cmd = 'fixImageXml.py -i demLat_*_*_Lon_*_*.dem.wgs84 -f' - runCmd(cmd) - cmd = 'rm *.hgt* *.log demLat_*_*_Lon_*_*.dem demLat_*_*_Lon_*_*.dem.vrt demLat_*_*_Lon_*_*.dem.xml' - runCmd(cmd) + # downloadUrl = 'http://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11' + # cmd = 'dem.py -a stitch -b {} -k -s 1 -c -f -u {}'.format( + # bboxStr, + # downloadUrl + # ) + # runCmd(cmd) + # cmd = 'fixImageXml.py -i demLat_*_*_Lon_*_*.dem.wgs84 -f' + # runCmd(cmd) + # cmd = 'rm *.hgt* *.log demLat_*_*_Lon_*_*.dem demLat_*_*_Lon_*_*.dem.vrt demLat_*_*_Lon_*_*.dem.xml' + # runCmd(cmd) + + #replace the above system calls with function calls + downloadDem(list(bbox), demType='version3', resolution=1, fillingValue=-32768, outputFile=None, userName=None, passWord=None) + imagePathXml((glob.glob('demLat_*_*_Lon_*_*.dem.wgs84'))[0], fullPath=True) + filesRemoved = glob.glob('*.hgt*') + glob.glob('*.log') + glob.glob('demLat_*_*_Lon_*_*.dem') + glob.glob('demLat_*_*_Lon_*_*.dem.vrt') + glob.glob('demLat_*_*_Lon_*_*.dem.xml') + for filex in filesRemoved: + os.remove(filex) + os.chdir('../') self.dem = glob.glob(os.path.join(demDir, 'demLat_*_*_Lon_*_*.dem.wgs84'))[0] @@ -54,16 +62,24 @@ def runDownloadDem(self): os.makedirs(demGeoDir, exist_ok=True) os.chdir(demGeoDir) - downloadUrl = 'http://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL3.003/2000.02.11' - cmd = 'dem.py -a stitch -b {} -k -s 3 -c -f -u {}'.format( - bboxStr, - downloadUrl - ) - runCmd(cmd) - cmd = 'fixImageXml.py -i demLat_*_*_Lon_*_*.dem.wgs84 -f' - runCmd(cmd) - cmd = 'rm *.hgt* *.log demLat_*_*_Lon_*_*.dem demLat_*_*_Lon_*_*.dem.vrt demLat_*_*_Lon_*_*.dem.xml' - runCmd(cmd) + # downloadUrl = 'http://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL3.003/2000.02.11' + # cmd = 'dem.py -a stitch -b {} -k -s 3 -c -f -u {}'.format( + # bboxStr, + # downloadUrl + # ) + # runCmd(cmd) + # cmd = 'fixImageXml.py -i demLat_*_*_Lon_*_*.dem.wgs84 -f' + # runCmd(cmd) + # cmd = 'rm *.hgt* *.log demLat_*_*_Lon_*_*.dem demLat_*_*_Lon_*_*.dem.vrt demLat_*_*_Lon_*_*.dem.xml' + # runCmd(cmd) + + #replace the above system calls with function calls + downloadDem(list(bbox), demType='version3', resolution=3, fillingValue=-32768, outputFile=None, userName=None, passWord=None) + imagePathXml((glob.glob('demLat_*_*_Lon_*_*.dem.wgs84'))[0], fullPath=True) + filesRemoved = glob.glob('*.hgt*') + glob.glob('*.log') + glob.glob('demLat_*_*_Lon_*_*.dem') + glob.glob('demLat_*_*_Lon_*_*.dem.vrt') + glob.glob('demLat_*_*_Lon_*_*.dem.xml') + for filex in filesRemoved: + os.remove(filex) + os.chdir('../') self.demGeo = glob.glob(os.path.join(demGeoDir, 'demLat_*_*_Lon_*_*.dem.wgs84'))[0] @@ -77,10 +93,17 @@ def runDownloadDem(self): #cmd = 'wbd.py {}'.format(bboxStr) #runCmd(cmd) download_wbd(np.int(np.floor(bbox[0])), np.int(np.ceil(bbox[1])), np.int(np.floor(bbox[2])), np.int(np.ceil(bbox[3]))) - cmd = 'fixImageXml.py -i swbdLat_*_*_Lon_*_*.wbd -f' - runCmd(cmd) - cmd = 'rm *.log' - runCmd(cmd) + #cmd = 'fixImageXml.py -i swbdLat_*_*_Lon_*_*.wbd -f' + #runCmd(cmd) + #cmd = 'rm *.log' + #runCmd(cmd) + + #replace the above system calls with function calls + imagePathXml((glob.glob('swbdLat_*_*_Lon_*_*.wbd'))[0], fullPath=True) + filesRemoved = glob.glob('*.log') + for filex in filesRemoved: + os.remove(filex) + os.chdir('../') self.wbd = glob.glob(os.path.join(wbdDir, 'swbdLat_*_*_Lon_*_*.wbd'))[0] @@ -94,6 +117,54 @@ def runDownloadDem(self): self._insar.procDoc.addAllFromCatalog(catalog) +def downloadDem(bbox, demType='version3', resolution=1, fillingValue=-32768, outputFile=None, userName=None, passWord=None): + ''' + bbox: [s, n, w, e] + demType: can be 'version3' or 'nasadem'. nasadem is also tested. + resolution: 1 or 3, NASADEM only available in 1-arc sec resolution + ''' + import numpy as np + import isceobj + from contrib.demUtils import createDemStitcher + + ds = createDemStitcher(demType) + ds.configure() + + if demType == 'version3': + if resolution == 1: + ds._url1 = 'http://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11' + else: + ds._url3 = 'http://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL3.003/2000.02.11' + elif demType == 'nasadem': + resolution = 1 + #this url is included in the module + #ds._url1 = 'http://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_HGT.001/2000.02.11' + else: + raise Exception('unknown DEM type, currently supported DEM types: version3 and nasadem') + + ds.setUsername(userName) + ds.setPassword(passWord) + + ds._keepAfterFailed = True + ds.setCreateXmlMetadata(True) + ds.setUseLocalDirectory(False) + ds.setFillingValue(fillingValue) + ds.setFilling() + + bbox = [np.int(np.floor(bbox[0])), np.int(np.ceil(bbox[1])), np.int(np.floor(bbox[2])), np.int(np.ceil(bbox[3]))] + if outputFile==None: + outputFile = ds.defaultName(bbox) + + if not(ds.stitchDems(bbox[0:2],bbox[2:4],resolution,outputFile,'./',keep=True)): + print('Could not create a stitched DEM. Some tiles are missing') + else: + #Apply correction EGM96 -> WGS84 + demImg = ds.correct() + + #report downloads + for k,v in list(ds._downloadReport.items()): + print(k,'=',v) + def download_wbd(s, n, w, e): ''' @@ -199,3 +270,21 @@ def download_wbd(s, n, w, e): return outputFile + + +def imagePathXml(imageFile, fullPath=True): + import os + import isceobj + from isceobj.Util.ImageUtil import ImageLib as IML + + img = IML.loadImage(imageFile)[0] + + dirname = os.path.dirname(imageFile) + if fullPath: + fname = os.path.abspath( os.path.join(dirname, os.path.basename(imageFile))) + else: + fname = os.path.basename(imageFile) + + img.filename = fname + img.setAccessMode('READ') + img.renderHdr() diff --git a/components/isceobj/Alos2Proc/runIonFilt.py b/components/isceobj/Alos2Proc/runIonFilt.py index 330faf2..41c531d 100644 --- a/components/isceobj/Alos2Proc/runIonFilt.py +++ b/components/isceobj/Alos2Proc/runIonFilt.py @@ -44,7 +44,8 @@ def runIonFilt(self): ################################### #SET PARAMETERS HERE #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self) - corThresholdAdj = 0.85 + corThresholdAdj = 0.97 + corOrderAdj = 20 ################################### print('\ncomputing ionosphere') @@ -75,11 +76,22 @@ def runIonFilt(self): upperUnw[area[0]:area[1], area[2]:area[3]] = 0 cor[area[0]:area[1], area[2]:area[3]] = 0 + #remove possible wired values in coherence + cor[np.nonzero(cor<0)] = 0.0 + cor[np.nonzero(cor>1)] = 0.0 + + #remove water body + wbd = np.fromfile('wbd'+ml2+'.wbd', dtype=np.int8).reshape(length, width) + cor[np.nonzero(wbd==-1)] = 0.0 + + #remove small values + cor[np.nonzero(cor= size_max: + print('\n\nWARNING: minimum window size for filtering ionosphere phase {} >= maximum window size {}'.format(size_min, size_max)) + print(' resetting maximum window size to {}\n\n'.format(size_min+5)) + size_max = size_min + 5 + #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self) - corThresholdIon = 0.85 + #corThresholdFit = 0.85 + + #Now changed to use lower band coherence. crl, 23-apr-2020. + useDiffCoherence = False + if useDiffCoherence: + #parameters for using diff coherence + corfile = 'diff'+ml2+'.cor' + corThresholdFit = 0.95 + # 1 is not good for low coherence case, changed to 20 + #corOrderFit = 1 + corOrderFit = 20 + corOrderFilt = 14 + else: + #parameters for using lower/upper band coherence + corfile = subbandPrefix[0]+ml2+'.cor' + corThresholdFit = 0.4 + corOrderFit = 10 + corOrderFilt = 4 + ################################################# print('\nfiltering ionosphere') ionfile = 'ion'+ml2+'.ion' - corfile = 'diff'+ml2+'.cor' + #corfile = 'diff'+ml2+'.cor' ionfiltfile = 'filt_ion'+ml2+'.ion' img = isceobj.createImage() @@ -137,6 +172,10 @@ def runIonFilt(self): cor[np.nonzero(cor<0)] = 0.0 cor[np.nonzero(cor>1)] = 0.0 + #remove water body + wbd = np.fromfile('wbd'+ml2+'.wbd', dtype=np.int8).reshape(length, width) + cor[np.nonzero(wbd==-1)] = 0.0 + # #applying water body mask here # waterBodyFile = 'wbd'+ml2+'.wbd' # if os.path.isfile(waterBodyFile): @@ -145,14 +184,17 @@ def runIonFilt(self): # cor[np.nonzero(wbd!=0)] = 0.00001 if fit: - ion_fit = weight_fitting(ion, cor, width, length, 1, 1, 1, 1, 2, corThresholdIon) + import copy + wgt = copy.deepcopy(cor) + wgt[np.nonzero(wgt=corThresholdAdj) + flag = (lowerUnw!=0)*(wgt!=0) index = np.nonzero(flag!=0) mv = np.mean((lowerUnw - upperUnw)[index], dtype=np.float64) print('mean value of phase difference: {}'.format(mv)) diff = mv #adjust phase using a surface else: - diff = weight_fitting(lowerUnw - upperUnw, cor, width, length, 1, 1, 1, 1, 2, corThresholdAdj) + diff = weight_fitting(lowerUnw - upperUnw, wgt, width, length, 1, 1, 1, 1, 2) flag2 = (lowerUnw!=0) index2 = np.nonzero(flag2) @@ -435,10 +476,10 @@ def cal_surface(x, y, c, order): return z -def weight_fitting(ionos, cor, width, length, nrli, nali, nrlo, nalo, order, coth): +def weight_fitting(ionos, weight, width, length, nrli, nali, nrlo, nalo, order): ''' ionos: input ionospheric phase - cor: coherence of the interferogram + weight: weight width: file width length: file length nrli: number of range looks of the input interferograms @@ -446,7 +487,6 @@ def weight_fitting(ionos, cor, width, length, nrli, nali, nrlo, nalo, order, cot nrlo: number of range looks of the output ionosphere phase nalo: number of azimuth looks of the ioutput ionosphere phase order: the order of the polynomial for fitting ionosphere phase estimates - coth: coherence threshhold for ionosphere phase estimation ''' from isceobj.Alos2Proc.Alos2ProcPublic import create_multi_index2 @@ -460,11 +500,8 @@ def weight_fitting(ionos, cor, width, length, nrli, nali, nrlo, nalo, order, cot rgindex = create_multi_index2(widtho, nrli, nrlo) azindex = create_multi_index2(lengtho, nali, nalo) - #convert coherence to weight - cor = cor**2/(1.009-cor**2) - #look for data to use - flag = (cor>coth)*(ionos!=0) + flag = (weight!=0)*(ionos!=0) point_index = np.nonzero(flag) m = point_index[0].shape[0] @@ -475,7 +512,7 @@ def weight_fitting(ionos, cor, width, length, nrli, nali, nrlo, nalo, order, cot x = x0[point_index].reshape(m, 1) y = y0[point_index].reshape(m, 1) z = ionos[point_index].reshape(m, 1) - w = cor[point_index].reshape(m, 1) + w = weight[point_index].reshape(m, 1) #convert to higher precision type before use x=np.asfarray(x,np.float64) diff --git a/components/isceobj/Alos2Proc/runIonUwrap.py b/components/isceobj/Alos2Proc/runIonUwrap.py index f9a35af..b8cb059 100644 --- a/components/isceobj/Alos2Proc/runIonUwrap.py +++ b/components/isceobj/Alos2Proc/runIonUwrap.py @@ -25,6 +25,7 @@ def runIonUwrap(self): masterTrack = self._insar.loadTrack(master=True) slaveTrack = self._insar.loadTrack(master=False) + wbdFile = os.path.abspath(self._insar.wbd) from isceobj.Alos2Proc.runIonSubband import defineIonDir ionDir = defineIonDir() @@ -40,6 +41,7 @@ def runIonUwrap(self): ############################################################ from isceobj.Alos2Proc.Alos2ProcPublic import create_xml from contrib.alos2proc.alos2proc import look + from isceobj.Alos2Proc.Alos2ProcPublic import waterBodyRadar ml2 = '_{}rlks_{}alks'.format(self._insar.numberRangeLooks1*self._insar.numberRangeLooksIon, self._insar.numberAzimuthLooks1*self._insar.numberAzimuthLooksIon) @@ -69,6 +71,14 @@ def runIonUwrap(self): # look(wbdOutFile, 'wbd'+ml2+'.wbd', width, self._insar.numberRangeLooksIon, self._insar.numberAzimuthLooksIon, 0, 0, 1) # create_xml('wbd'+ml2+'.wbd', width2, length2, 'byte') + #water body + if k == 0: + look(os.path.join(fullbandDir, self._insar.latitude), 'lat'+ml2+'.lat', width, self._insar.numberRangeLooksIon, self._insar.numberAzimuthLooksIon, 3, 0, 1) + look(os.path.join(fullbandDir, self._insar.longitude), 'lon'+ml2+'.lon', width, self._insar.numberRangeLooksIon, self._insar.numberAzimuthLooksIon, 3, 0, 1) + create_xml('lat'+ml2+'.lat', width2, length2, 'double') + create_xml('lon'+ml2+'.lon', width2, length2, 'double') + waterBodyRadar('lat'+ml2+'.lat', 'lon'+ml2+'.lon', wbdFile, 'wbd'+ml2+'.wbd') + ############################################################ # STEP 2. compute coherence diff --git a/components/isceobj/Alos2Proc/runLook.py b/components/isceobj/Alos2Proc/runLook.py index 70c4aa0..f6ff049 100644 --- a/components/isceobj/Alos2Proc/runLook.py +++ b/components/isceobj/Alos2Proc/runLook.py @@ -50,8 +50,21 @@ def runLook(self): create_xml(self._insar.multilookLongitude, width2, length2, 'double') create_xml(self._insar.multilookHeight, width2, length2, 'double') #los has two bands, use look program in isce instead - cmd = "looks.py -i {} -o {} -r {} -a {}".format(self._insar.los, self._insar.multilookLos, self._insar.numberRangeLooks2, self._insar.numberAzimuthLooks2) - runCmd(cmd) + #cmd = "looks.py -i {} -o {} -r {} -a {}".format(self._insar.los, self._insar.multilookLos, self._insar.numberRangeLooks2, self._insar.numberAzimuthLooks2) + #runCmd(cmd) + + #replace the above system call with function call + from mroipac.looks.Looks import Looks + from isceobj.Image import createImage + inImage = createImage() + inImage.load(self._insar.los+'.xml') + + lkObj = Looks() + lkObj.setDownLooks(self._insar.numberAzimuthLooks2) + lkObj.setAcrossLooks(self._insar.numberRangeLooks2) + lkObj.setInputImage(inImage) + lkObj.setOutputFilename(self._insar.multilookLos) + lkObj.looks() #water body #this looking operation has no problems where there is only water and land, but there is also possible no-data area diff --git a/components/isceobj/Alos2Proc/runSwathMosaic.py b/components/isceobj/Alos2Proc/runSwathMosaic.py index 8dad041..9d6b175 100644 --- a/components/isceobj/Alos2Proc/runSwathMosaic.py +++ b/components/isceobj/Alos2Proc/runSwathMosaic.py @@ -378,33 +378,94 @@ def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num #get difference - corth = 0.87 - if filt: - corth = 0.90 - diffMean0 = 0.0 - for k in range(5): - dataDiff = data1 * np.conj(data2) - cor = cal_coherence_1(dataDiff, win=3) - index = np.nonzero(np.logical_and(cor>corth, dataDiff!=0)) + dataDiff = data1 * np.conj(data2) + cor = cal_coherence_1(dataDiff, win=5) + index = np.nonzero(np.logical_and(cor>0.85, dataDiff!=0)) + + DEBUG=False + if DEBUG: + from isceobj.Alos2Proc.Alos2ProcPublic import create_xml + (length7, width7)=dataDiff.shape + filename = 'diff_ori_s{}-s{}.int'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber) + dataDiff.astype(np.complex64).tofile(filename) + create_xml(filename, width7, length7, 'int') + filename = 'cor_ori_s{}-s{}.cor'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber) + cor.astype(np.float32).tofile(filename) + create_xml(filename, width7, length7, 'float') + + print('\ncompute phase difference between subswaths {} and {}'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber)) + print('number of pixels with coherence > 0.85: {}'.format(index[0].size)) + + #if already filtered the subswath overlap interferograms (MAI), do not filtered differential interferograms + if (filt == False) and (index[0].size < 4000): + #coherence too low, filter subswath overlap differential interferogram + diffMean0 = 0.0 + breakFlag = False + for (filterStrength, filterWinSize) in zip([3.0, 9.0], [64, 128]): + dataDiff = data1 * np.conj(data2) + dataDiff /= (np.absolute(dataDiff)+(dataDiff==0)) + dataDiff = filterInterferogram(dataDiff, filterStrength, filterWinSize, 1) + cor = cal_coherence_1(dataDiff, win=7) + + DEBUG=False + if DEBUG: + from isceobj.Alos2Proc.Alos2ProcPublic import create_xml + (length7, width7)=dataDiff.shape + filename = 'diff_filt_s{}-s{}_strength_{}_winsize_{}.int'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, filterStrength, filterWinSize) + dataDiff.astype(np.complex64).tofile(filename) + create_xml(filename, width7, length7, 'int') + filename = 'cor_filt_s{}-s{}_strength_{}_winsize_{}.cor'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, filterStrength, filterWinSize) + cor.astype(np.float32).tofile(filename) + create_xml(filename, width7, length7, 'float') + + for corth in [0.99999, 0.9999]: + index = np.nonzero(np.logical_and(cor>corth, dataDiff!=0)) + if index[0].size > 30000: + breakFlag = True + break + if breakFlag: + break + if index[0].size < 100: diffMean0 = 0.0 - print('\n\nWARNING: too few high coherence pixels for swath phase difference estimation between swath {} and {}'.format(i-1, i)) - print(' : first swath swath number: 0\n\n') - break - angle = np.mean(np.angle(dataDiff[index]), dtype=np.float64) - diffMean0 += angle - data2 *= np.exp(np.complex64(1j) * angle) - print('phase offset: %15.12f rad after loop: %3d'%(diffMean0, k)) + print('\n\nWARNING: too few high coherence pixels for swath phase difference estimation') + print(' number of high coherence pixels: {}\n\n'.format(index[0].size)) + else: + print('filtered coherence threshold used: {}, number of pixels used: {}'.format(corth, index[0].size)) + angle = np.mean(np.angle(dataDiff[index]), dtype=np.float64) + diffMean0 += angle + data2 *= np.exp(np.complex64(1j) * angle) + print('phase offset: %15.12f rad with filter strength: %f, window size: %3d'%(diffMean0, filterStrength, filterWinSize)) + else: + diffMean0 = 0.0 + for k in range(30): + dataDiff = data1 * np.conj(data2) + cor = cal_coherence_1(dataDiff, win=5) + if filt: + index = np.nonzero(np.logical_and(cor>0.95, dataDiff!=0)) + else: + index = np.nonzero(np.logical_and(cor>0.85, dataDiff!=0)) + if index[0].size < 100: + diffMean0 = 0.0 + print('\n\nWARNING: too few high coherence pixels for swath phase difference estimation') + print(' number of high coherence pixels: {}\n\n'.format(index[0].size)) + break + angle = np.mean(np.angle(dataDiff[index]), dtype=np.float64) + diffMean0 += angle + data2 *= np.exp(np.complex64(1j) * angle) + print('phase offset: %15.12f rad after loop: %3d'%(diffMean0, k)) + + DEBUG=False + if DEBUG and (k==0): + from isceobj.Alos2Proc.Alos2ProcPublic import create_xml + (length7, width7)=dataDiff.shape + filename = 'diff_ori_s{}-s{}_loop_{}.int'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, k) + dataDiff.astype(np.complex64).tofile(filename) + create_xml(filename, width7, length7, 'int') + filename = 'cor_ori_s{}-s{}_loop_{}.cor'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, k) + cor.astype(np.float32).tofile(filename) + create_xml(filename, width7, length7, 'float') - DEBUG=False - if DEBUG and (k==0): - from isceobj.Alos2Proc.Alos2ProcPublic import create_xml - (lengthxx, widthxx)=dataDiff.shape - filtnamePrefix = 'subswath{}_subswath{}_loop{}'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, k) - cor.astype(np.float32).tofile(filtnamePrefix+'.cor') - create_xml(filtnamePrefix+'.cor', widthxx, lengthxx, 'float') - dataDiff.astype(np.complex64).tofile(filtnamePrefix+'.int') - create_xml(filtnamePrefix+'.int', widthxx, lengthxx, 'int') diffMean.append(diffMean0) print('phase offset: subswath{} - subswath{}: {}'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, diffMean0)) diff --git a/configuration/buildHelper.py b/configuration/buildHelper.py index 70a502c..138bdb3 100755 --- a/configuration/buildHelper.py +++ b/configuration/buildHelper.py @@ -7,7 +7,11 @@ tmpdump = 'tmpdump.json' def createHelp(env,factoryFile,installDir): #jng: try to have scons handle all the creation but could not figure out how # so handled dir creation manually - os.makedirs(env['HELPER_BUILD_DIR'], exist_ok=True) + try: + os.makedirs(env['HELPER_BUILD_DIR']) + except: + # already exists + pass try: #one could probably also use __import__ but needs to make sure the #the cwd is prepended to the sys.path otherwise if factoryFile = __init__.py diff --git a/contrib/alos2filter/src/psfilt1.c b/contrib/alos2filter/src/psfilt1.c index 941e082..3051f1f 100644 --- a/contrib/alos2filter/src/psfilt1.c +++ b/contrib/alos2filter/src/psfilt1.c @@ -2,6 +2,7 @@ #include #include #include +#include #define PLUS 1 #define MINU 2 @@ -75,6 +76,11 @@ int psfilt1(char *inputfile, char *outputfile, int width, double alpha, int fftw FILE *int_file, *sm_file; + double psd,psd_sc; /* power spectrum, scale factor */ + int ii, jj; + fftwf_plan p_forward; + fftwf_plan p_backward; + int k; float sf; // scale factor for FFT, otherwise FFT will magnify the data by FFT length // usually the magnitude of the interferogram is very large, so the data are @@ -212,6 +218,9 @@ int psfilt1(char *inputfile, char *outputfile, int width, double alpha, int fftw } lc=0; + p_forward = fftwf_plan_dft_2d(fftw, fftw, (fftw_complex *)seg_fft[0], (fftw_complex *)seg_fft[0], FFTW_FORWARD, FFTW_MEASURE); + p_backward = fftwf_plan_dft_2d(fftw, fftw, (fftw_complex *)seg_fft[0], (fftw_complex *)seg_fft[0], FFTW_BACKWARD, FFTW_MEASURE); + for (i=0; i < yh; i += step){ for(i1=fftw - step; i1 < fftw; i1++){ fread((char *)cmp[i1], sizeof(fcomplex), width, int_file); @@ -229,8 +238,30 @@ int psfilt1(char *inputfile, char *outputfile, int width, double alpha, int fftw if(i%(2*step) == 0)fprintf(stderr,"\rprogress: %3d%%", (int)(i*100/yh + 0.5)); for (j=0; j < width; j += step){ - psd_wgt(cmp, seg_fft, alpha, j, i, fftw); - fourn((float *)seg_fft[0]-1,nfft,ndim,isign); /* 2D inverse FFT of region, get back filtered fringes */ + //psd_wgt(cmp, seg_fft, alpha, j, i, fftw); + //////////////////////////////////////////////////////////////////////////////////////// + //replace function psd_wgt with the following to call FFTW, crl, 23-APR-2020 + for (ii=0; ii < fftw; ii++){ /* load up data array */ + for (jj=j; jj < j+fftw; jj++){ + seg_fft[ii][jj-j].re = cmp[ii][jj].re; + seg_fft[ii][jj-j].im = cmp[ii][jj].im; + } + } + + //fourn((float *)seg_fft[0]-1, nfft, ndim, -1); /* 2D forward FFT of region */ + fftwf_execute(p_forward); + + for (ii=0; ii < fftw; ii++){ + for (jj=0; jj < fftw; jj++){ + psd = seg_fft[ii][jj].re * seg_fft[ii][jj].re + seg_fft[ii][jj].im * seg_fft[ii][jj].im; + psd_sc = pow(psd,alpha/2.); + seg_fft[ii][jj].re *= psd_sc; + seg_fft[ii][jj].im *= psd_sc; + } + } + ///////////////////////////////////////////////////////////////////////////////////////// + //fourn((float *)seg_fft[0]-1,nfft,ndim,isign); /* 2D inverse FFT of region, get back filtered fringes */ + fftwf_execute(p_backward); for (i1=0; i1 < fftw; i1++){ /* save filtered output values */ for (j1=0; j1 < fftw; j1++){ @@ -285,6 +316,9 @@ int psfilt1(char *inputfile, char *outputfile, int width, double alpha, int fftw fprintf(stdout,"\nnumber of lines written to file: %d\n",lc); stop_timing(); + fftwf_destroy_plan(p_forward); + fftwf_destroy_plan(p_backward); + free(bufcz); //free_cmatrix(cmp, 0, fftw-1, -fftw,width+fftw); //free_cmatrix(sm, 0,fftw-1,-fftw,width+fftw); diff --git a/contrib/demUtils/__init__.py b/contrib/demUtils/__init__.py index 999dd89..4ecb1df 100755 --- a/contrib/demUtils/__init__.py +++ b/contrib/demUtils/__init__.py @@ -32,6 +32,8 @@ def createDemStitcher(type='version3',name = ''): from contrib.demUtils.DemStitcherV3 import DemStitcher if(type == 'version2'): from contrib.demUtils.DemStitcher import DemStitcher + if(type == 'nasadem'): + from contrib.demUtils.DemStitcherND import DemStitcher return DemStitcher(name=name) def createSWBDStitcher(name = ''): from contrib.demUtils.SWBDStitcher import SWBDStitcher @@ -56,4 +58,4 @@ def getFactoriesInfo(): { 'factory':'createCorrect_geoid_i2_srtm' } - } \ No newline at end of file + } diff --git a/contrib/demUtils/demstitcher/DemStitcher.py b/contrib/demUtils/demstitcher/DemStitcher.py index ff36cce..0cbc6f1 100755 --- a/contrib/demUtils/demstitcher/DemStitcher.py +++ b/contrib/demUtils/demstitcher/DemStitcher.py @@ -763,6 +763,7 @@ class DemStitcher(Component): def getUnzippedName(self,name,source = None): return name.replace(self._zip,'') def stitchDems(self,lat,lon,source, outname, downloadDir = None,region = None, keep = None, swap = None): + import glob if downloadDir is None: downloadDir = self._downloadDir else: @@ -839,7 +840,8 @@ class DemStitcher(Component): if not self._keepDems: for dem in decompressedList: - os.remove(dem) + for d in glob.glob('*'+os.path.basename(dem.decode('UTF-8'))[:7]+'*'): + os.remove(d) if self._createXmlMetadata: self.createXmlMetadata(lat,lon,source,outname) if self._createRscMetadata: diff --git a/contrib/demUtils/demstitcher/DemStitcherND.py b/contrib/demUtils/demstitcher/DemStitcherND.py new file mode 100755 index 0000000..5597185 --- /dev/null +++ b/contrib/demUtils/demstitcher/DemStitcherND.py @@ -0,0 +1,271 @@ +#!/usr/bin/env python3 + + +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Copyright 2012 California Institute of Technology. ALL RIGHTS RESERVED. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# United States Government Sponsorship acknowledged. This software is subject to +# U.S. export control laws and regulations and has been classified as 'EAR99 NLR' +# (No [Export] License Required except when exporting to an embargoed country, +# end user, or in support of a prohibited end use). By downloading this software, +# the user agrees to comply with all applicable U.S. export laws and regulations. +# The user has the responsibility to obtain export licenses, or other export +# authority as may be required before exporting this software to any 'EAR99' +# embargoed foreign country or citizen of those countries. +# +# Author: Giangi Sacco +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + + + + + + +import isce +from ctypes import cdll +import os +import sys +import urllib.request, urllib.error, urllib.parse +from isce import logging +from iscesys.Component.Component import Component +from contrib.demUtils.DemStitcher import DemStitcher as DS +#Parameters definitions +URL1 = Component.Parameter('_url1', + public_name = 'URL1',default = 'http://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_HGT.001/2000.02.11', + type = str, + mandatory = False, + doc = "Url for the high resolution DEM. Used for NASADEM") +EXTRA_PREPEND1 = Component.Parameter('_extraPrepend1', + public_name = 'extra prepend 1',default = 'NASADEM_HGT', + type = str, + mandatory = False, + doc = "The actual file name might have some extra string compared to the conventional one." \ + + "This is for the high resolution files. Used for NASADEM") +HAS_EXTRAS = Component.Parameter('_hasExtras', + public_name = 'has extras',default = True, + type = bool, + mandatory = False, + doc = "Instead of having to provide the EXTRA_EXT or EXTRA_PREPEND empty when the extra extension " \ ++ "is not present, turn on this flag. Used for NASADEM") + +## This class provides a set of convenience method to retrieve and possibly combine different DEMs from the USGS server. +# \c NOTE: the latitudes and the longitudes that describe the DEMs refer to the bottom left corner of the image. +class DemStitcher(DS): + + + ## + # Given a latitude and longitude in degrees it returns the expected filename. + # @param lat \c int latitude in the range (-90,90). Actual data are restricted to (-60,60) or so. + # @param lon \c int longitude in the range [-180,180) or [0,360). + # @return \c string the filename for that location + + def createFilename(self,lat,lon,source = None): + + if lon > 180: + lon = -(360 - lon) + else: + lon = lon + ns,ew = self.convertCoordinateToString(lat,lon) + #make coords lower-case + ns=ns.lower() + ew=ew.lower() + if(self._hasExtras): + if(source and source == 1): + toPrepend = self._extraPrepend1 + "_" + else: + print('Unrecognized dem source',source) + raise Exception + + return toPrepend+ ns + ew + self._extension + self._zip + + else: + return ns + ew + self._extension + self._zip + + + def getUnzippedName(self,name,source = None): + name = name.replace(self._extraPrepend1+'_','') + return name.replace(self._zip,'.hgt') + + ## + # Given a list of filenames it fetches the corresponding + # compressed (zip format) DEMs. + # @param source \c int the type of DEM. source = 1 for 1 arcsec resolution data, + # source = 3 for 3 arcsec resolution data. + # @param listFile \c list of the filenames to be retrieved. + # @param downloadDir \c string the directory where the DEMs are downloaded. + # If the directory does not exists it will be created. If the argument is not + # provided then the files are downloaded in the location defined by the + # self._downloadDir that is defaulted to the current directory. + # @param region \c list \c strings regions where to look for the files. It must + # have the same length of \c listFile. If not provided the files are searched by + # scanning the content of each region. Use method getRegionList to get the list of + # possible regions for a given source. Set region only if sure that all the requested + # file are contained in it. + + def getDems(self,source,listFile,downloadDir = None,region = None): + if downloadDir is None: + downloadDir = self._downloadDir + else: + self._downloadDir = downloadDir + + if not (downloadDir) is None: + try: + os.makedirs(downloadDir) + except: + #dir already exists + pass + for fileNow in listFile: + url = self.getFullHttp(source) + opener = urllib.request.URLopener() + try: + if not os.path.exists(os.path.join(downloadDir,fileNow)): + if(self._un is None or self._pw is None): + #opener.retrieve(url + fileNow,os.path.join(downloadDir,fileNow)) + if os.path.exists(os.path.join(os.environ['HOME'],'.netrc')): + command = 'curl -n -L -c $HOME/.earthdatacookie -b $HOME/.earthdatacookie -k -f -O ' + os.path.join(url,fileNow) + else: + self.logger.error('Please create a .netrc file in your home directory containing\nmachine urs.earthdata.nasa.gov\n\tlogin yourusername\n\tpassword yourpassword') + sys.exit(1) + else: + command = 'curl -k -f -u ' + self._un + ':' + self._pw + ' -O ' + os.path.join(url,fileNow) + # curl with -O download in working dir, so save current, move to donwloadDir + # nd get back once download is finished + cwd = os.getcwd() + os.chdir(downloadDir) + print(command) + if os.system(command): + os.chdir(cwd) + raise Exception + os.chdir(cwd) + self._downloadReport[fileNow] = self._succeded + except Exception as e: + self.logger.warning('There was a problem in retrieving the file %s. Exception %s'%(os.path.join(url,fileNow),str(e))) + self._downloadReport[fileNow] = self._failed + + + #still need to call it since the initialization calls the _url so the setter of + #url does not get called + def _configure(self): + pass + + def _facilities(self): + super(DemStitcher,self)._facilities() + + def __setstate__(self,d): + self.__dict__.update(d) + self.logger = logging.getLogger('isce.contrib.demUtils.DemStitcherV3') + libName = os.path.join(os.path.dirname(__file__),self._loadLibName) + ##self._keepAfterFailed = False #if True keeps the downloaded files even if the stitching failed. + self._lib = cdll.LoadLibrary(libName) + return + + def main(self): + # prevent from deliting local files + if(self._useLocalDirectory): + self._keepAfterFailed = True + self._keepDems = True + # is a metadata file is created set the right type + if(self._meta == 'xml'): + self.setCreateXmlMetadata(True) + elif(self._meta == 'rsc'): + self.setCreateRscMetadata(True) + # check for the action to be performed + if(self._action == 'stitch'): + if(self._bbox): + lat = self._bbox[0:2] + lon = self._bbox[2:4] + if (self._outputFile is None): + self._outputFile = self.defaultName(self._bbox) + + if not(self.stitchDems(lat,lon,self._source,self._outputFile,self._downloadDir, \ + keep=self._keepDems)): + print('Could not create a stitched DEM. Some tiles are missing') + else: + if(self._correct): + width = self.getDemWidth(lon,self._source) + self.correct(os.path.join(self._downloadDir,self._outputFile), \ + self._source,width,min(lat[0],lat[1]),min(lon[0],lon[1])) + #self.correct(self._output,self._source,width,min(lat[0],lat[1]),min(lon[0],lon[1])) + else: + print('Error. The "bbox" attribute must be specified when the action is "stitch"') + raise ValueError + elif(self._action == 'download'): + if(self._bbox): + lat = self._bbox[0:2] + lon = self._bbox[2:4] + self.getDemsInBox(lat,lon,self._source,self._downloadDir) + #can make the bbox and pairs mutually esclusive if replace the if below with elif + if(self._pairs): + self.downloadFilesFromList(self._pairs[::2],self._pairs[1::2],self._source,self._downloadDir) + if(not (self._bbox or self._pairs)): + print('Error. Either the "bbox" attribute or the "pairs" attribute must be specified when --action download is used') + raise ValueError + + else: + print('Unrecognized action ',self._action) + return + + if(self._report): + for k,v in list(self._downloadReport.items()): + print(k,'=',v) + + #use this logic so the right http is returned + + def getFullHttp(self,source): + return self._url1 if source == 1 else self._url3 + ''' + parameter_list = ( + URL1, + EXTRA_PREPEND1, + HAS_EXTRAS, + USERNAME, + PASSWORD, + KEEP_AFTER_FAILED, + DIRECTORY, + ACTION, + CORRECT, + META, + SOURCE, + NO_FILLING, + FILLING_VALUE, + BBOX, + PAIRS, + KEEP_DEMS, + REPORT, + USE_LOCAL_DIRECTORY, + OUTPUT_FILE + ) + ''' + parameter_list = ( + URL1, + EXTRA_PREPEND1, + HAS_EXTRAS + ) + DS.parameter_list + + family = 'demstitcher' + + def __init__(self,family = '', name = ''): + + super(DemStitcher, self).__init__(family if family else self.__class__.family, name=name) + # logger not defined until baseclass is called + self._extension = '' + self._zip = '.zip' + + #to make it working with other urls, make sure that the second part of the url + #it's /srtm/version2_1/SRTM(1,3) + self._remove = ['.jpg','.xml'] + if not self.logger: + self.logger = logging.getLogger('isce.contrib.demUtils.DemStitcherV3') diff --git a/contrib/demUtils/demstitcher/SConscript b/contrib/demUtils/demstitcher/SConscript index cd12597..91ee6cb 100644 --- a/contrib/demUtils/demstitcher/SConscript +++ b/contrib/demUtils/demstitcher/SConscript @@ -50,6 +50,6 @@ helpList,installHelp = envdemStitcher['HELP_BUILDER'](envdemStitcher,'../__init_ envdemStitcher.Install(installHelp,helpList) envdemStitcher.Alias('install',installHelp) -listFiles = ['DemStitcher.py','DemStitcherV3.py'] +listFiles = ['DemStitcher.py','DemStitcherND.py','DemStitcherV3.py'] envdemStitcher.Install(install,listFiles) envdemStitcher.Alias('install',install) diff --git a/contrib/stack/topsStack/Stack.py b/contrib/stack/topsStack/Stack.py index fc8ecb6..54603af 100644 --- a/contrib/stack/topsStack/Stack.py +++ b/contrib/stack/topsStack/Stack.py @@ -223,6 +223,7 @@ class config(object): self.f.write('defomax : ' + self.defoMax + '\n') self.f.write('rlks : ' + self.rangeLooks + '\n') self.f.write('alks : ' + self.azimuthLooks + '\n') + self.f.write('rmfilter : ' + self.rmFilter + '\n') self.f.write('method : ' + self.unwMethod + '\n') def unwrapSnaphu(self, function): @@ -282,7 +283,7 @@ class run(object): os.makedirs(self.config_path, exist_ok=True) for slcdate in acquisitionDates: - configName = os.path.join(self.config_path,'config_'+slcdate) + configName = os.path.join(self.config_path,'config_unpack_'+slcdate) configObj = config(configName) configObj.configure(self) configObj.dirName = safe_dict[slcdate].safe_file @@ -377,11 +378,15 @@ class run(object): self.runf.write(self.text_cmd + 'subsetMaster.py -m ' + os.path.join(self.work_dir, 'master') + ' -g ' + os.path.join(self.work_dir, 'geom_master') + '\n') - def overlap_geo2rdr_resample(self, slaveList): + + def geo2rdr_offset(self, slaveList, fullBurst='False'): for slave in slaveList: master = self.master_date - configName = os.path.join(self.config_path ,'config_resamp_overlap_'+slave) + if fullBurst == 'True': + configName = os.path.join(self.config_path, 'config_fullBurst_geo2rdr_' + slave) + else: + configName = os.path.join(self.config_path, 'config_overlap_geo2rdr_'+slave) ########### configObj = config(configName) configObj.configure(self) @@ -389,16 +394,42 @@ class run(object): configObj.masterDir = os.path.join(self.work_dir, 'master') configObj.geom_master = os.path.join(self.work_dir, 'geom_master') configObj.coregSlaveDir = os.path.join(self.work_dir, 'coreg_slaves/'+slave) - configObj.overlapTrueOrFalse = 'True' + if fullBurst == 'True': + configObj.misreg_az = os.path.join(self.work_dir, 'misreg/azimuth/dates/' + slave + '.txt') + configObj.misreg_rng = os.path.join(self.work_dir, 'misreg/range/dates/' + slave + '.txt') + configObj.overlapTrueOrFalse = 'False' + else: + configObj.overlapTrueOrFalse = 'True' configObj.geo2rdr('[Function-1]') - ########### - configObj.interferogram_prefix = 'coarse' - configObj.masterDir = os.path.join(self.work_dir,'master') - configObj.resamp_withCarrier('[Function-2]') - ########### configObj.finalize() del configObj - self.runf.write(self.text_cmd + 'SentinelWrapper.py -c ' + configName + '\n') + self.runf.write(self.text_cmd + 'SentinelWrapper.py -c ' + configName + '\n') + + def resample_with_carrier(self, slaveList, fullBurst='False'): + for slave in slaveList: + master = self.master_date + if fullBurst == 'True': + configName = os.path.join(self.config_path, 'config_fullBurst_resample_' + slave) + else: + configName = os.path.join(self.config_path, 'config_overlap_resample_' + slave) + ########### + configObj = config(configName) + configObj.configure(self) + configObj.slaveDir = os.path.join(self.work_dir, 'slaves/' + slave) + configObj.masterDir = os.path.join(self.work_dir, 'master') + configObj.coregSlaveDir = os.path.join(self.work_dir, 'coreg_slaves/' + slave) + configObj.interferogram_prefix = 'coarse' + configObj.masterDir = os.path.join(self.work_dir, 'master') + if fullBurst == 'True': + configObj.misreg_az = os.path.join(self.work_dir, 'misreg/azimuth/dates/' + slave + '.txt') + configObj.misreg_rng = os.path.join(self.work_dir, 'misreg/range/dates/' + slave + '.txt') + configObj.overlapTrueOrFalse = 'False' + else: + configObj.overlapTrueOrFalse = 'True' + configObj.resamp_withCarrier('[Function-1]') + configObj.finalize() + del configObj + self.runf.write(self.text_cmd + 'SentinelWrapper.py -c ' + configName + '\n') def pairs_misregistration(self, dateList, safe_dict): # generating overlap interferograms, estimate azimuth misregistration for each pair: @@ -453,37 +484,13 @@ class run(object): self.runf.write(self.text_cmd + 'invertMisreg.py -i ' + os.path.join(self.work_dir,'misreg/azimuth/pairs/') + ' -o ' + os.path.join(self.work_dir,'misreg/azimuth/dates/') + '\n') self.runf.write(self.text_cmd + 'invertMisreg.py -i ' + os.path.join(self.work_dir,'misreg/range/pairs/') + ' -o ' + os.path.join(self.work_dir,'misreg/range/dates/') + '\n') - def geo2rdr_resample(self, slaveList): - # geometry offsets and resampling each full burst slave SLC - for slave in slaveList: - master = self.master_date - configName = os.path.join(self.config_path , 'config_resamp_' + slave) - ########### - configObj = config(configName) - configObj.configure(self) - configObj.slaveDir = os.path.join(self.work_dir, 'slaves/' + slave) - configObj.masterDir = os.path.join(self.work_dir, 'master') - configObj.geom_master = os.path.join(self.work_dir, 'geom_master') - configObj.coregSlaveDir = os.path.join(self.work_dir, 'coreg_slaves/' + slave) - configObj.misreg_az = os.path.join(self.work_dir, 'misreg/azimuth/dates/' + slave + '.txt') - configObj.misreg_rng = os.path.join(self.work_dir, 'misreg/range/dates/' + slave + '.txt') - configObj.overlapTrueOrFalse = 'False' - configObj.geo2rdr('[Function-1]') - ########### - configObj.interferogram_prefix = 'coarse' - configObj.masterDir = os.path.join(self.work_dir, 'master') - configObj.resamp_withCarrier('[Function-2]') - ########### - configObj.finalize() - del configObj - self.runf.write(self.text_cmd + 'SentinelWrapper.py -c ' + configName + '\n') - def extractStackValidRegion(self): masterDir = os.path.join(self.work_dir, 'master') coregSlaveDir = os.path.join(self.work_dir, 'coreg_slaves') self.runf.write(self.text_cmd + 'extractCommonValidRegion.py -m ' + masterDir + ' -s ' + coregSlaveDir + '\n') - def burstIgram_mergeBurst(self, dateList, safe_dict, pairs): + def generate_burstIgram(self, dateList, safe_dict, pairs): + for date in dateList: safe_dict[date].slc = os.path.join(self.work_dir, 'coreg_slaves/'+date) safe_dict[self.master_date].slc = os.path.join(self.work_dir , 'master') @@ -491,8 +498,7 @@ class run(object): master = pair[0] slave = pair[1] interferogramDir = os.path.join(self.work_dir, 'interferograms/' + master + '_' + slave) - mergedDir = os.path.join(self.work_dir, 'merged/interferograms/' + master + '_' + slave) - configName = os.path.join(self.config_path ,'config_igram_' + master + '_' + slave) + configName = os.path.join(self.config_path ,'config_generate_igram_' + master + '_' + slave) configObj = config(configName) configObj.configure(self) configObj.masterDir = safe_dict[master].slc @@ -502,9 +508,26 @@ class run(object): configObj.flatten = 'False' configObj.overlapTrueOrFalse = 'False' configObj.generateIgram('[Function-1]') + configObj.finalize() + del configObj - configObj.master = configObj.interferogramDir - configObj.dirName = configObj.master + self.runf.write(self.text_cmd + 'SentinelWrapper.py -c ' + configName + '\n') + + def igram_mergeBurst(self, dateList, safe_dict, pairs): + for date in dateList: + safe_dict[date].slc = os.path.join(self.work_dir, 'coreg_slaves/'+date) + safe_dict[self.master_date].slc = os.path.join(self.work_dir , 'master') + for pair in pairs: + master = pair[0] + slave = pair[1] + interferogramDir = os.path.join(self.work_dir, 'interferograms/' + master + '_' + slave) + mergedDir = os.path.join(self.work_dir, 'merged/interferograms/' + master + '_' + slave) + configName = os.path.join(self.config_path ,'config_merge_igram_' + master + '_' + slave) + configObj = config(configName) + configObj.configure(self) + configObj.interferogram_prefix = 'fine' + configObj.master = interferogramDir + configObj.dirName = interferogramDir configObj.namePattern = 'fine*int' configObj.mergedFile = mergedDir + '/' + configObj.interferogram_prefix + '.int' configObj.mergeBurstsMethod = 'top' @@ -513,7 +536,7 @@ class run(object): configObj.useVirtualFiles = 'True' configObj.multiLook = 'True' configObj.stack = os.path.join(self.work_dir, 'stack') - configObj.mergeBurst('[Function-2]') + configObj.mergeBurst('[Function-1]') configObj.finalize() del configObj @@ -657,6 +680,7 @@ class run(object): configObj.cohName = os.path.join(mergedDir,'filt_fine.cor') configObj.unwName = os.path.join(mergedDir,'filt_fine.unw') configObj.noMCF = noMCF + configObj.rmFilter = self.rmFilter configObj.master = os.path.join(self.work_dir,'master') configObj.defoMax = defoMax configObj.unwMethod = self.unwMethod diff --git a/contrib/stack/topsStack/stackSentinel.py b/contrib/stack/topsStack/stackSentinel.py index 0103500..c70b819 100755 --- a/contrib/stack/topsStack/stackSentinel.py +++ b/contrib/stack/topsStack/stackSentinel.py @@ -162,6 +162,9 @@ def createParser(): parser.add_argument('-useGPU', '--useGPU', dest='useGPU',action='store_true', default=False, help='Allow App to use GPU when available') + parser.add_argument('-rmFilter', '--rmFilter', dest='rmFilter', action='store_true', default=False, + help='Make an extra unwrap file in which filtering effect is removed') + return parser def cmdLineParse(iargs = None): @@ -457,10 +460,16 @@ def slcStack(inps, acquisitionDates, stackMasterDate, slaveDates, safe_dict, upd runObj.extractOverlaps() runObj.finalize() - i+=1 + i += 1 runObj = run() - runObj.configure(inps, 'run_' + str(i) + "_overlap_geo2rdr_resample") - runObj.overlap_geo2rdr_resample(slaveDates) + runObj.configure(inps, 'run_' + str(i) + "_overlap_geo2rdr") + runObj.geo2rdr_offset(slaveDates) + runObj.finalize() + + i += 1 + runObj = run() + runObj.configure(inps, 'run_' + str(i) + "_overlap_resample") + runObj.resample_with_carrier(slaveDates) runObj.finalize() i+=1 @@ -478,10 +487,16 @@ def slcStack(inps, acquisitionDates, stackMasterDate, slaveDates, safe_dict, upd runObj.timeseries_misregistration() runObj.finalize() - i+=1 + i += 1 runObj = run() - runObj.configure(inps, 'run_' + str(i) + "_geo2rdr_resample") - runObj.geo2rdr_resample(slaveDates) + runObj.configure(inps, 'run_' + str(i) + "_fullBurst_geo2rdr") + runObj.geo2rdr_offset(slaveDates, fullBurst='True') + runObj.finalize() + + i += 1 + runObj = run() + runObj.configure(inps, 'run_' + str(i) + "_fullBurst_resample") + runObj.resample_with_carrier(slaveDates, fullBurst='True') runObj.finalize() i+=1 @@ -546,8 +561,14 @@ def interferogramStack(inps, acquisitionDates, stackMasterDate, slaveDates, safe i+=1 runObj = run() - runObj.configure(inps, 'run_' + str(i) + "_merge_burst_igram") - runObj.burstIgram_mergeBurst(acquisitionDates, safe_dict, pairs) + runObj.configure(inps, 'run_' + str(i) + "_generate_burst_igram") + runObj.generate_burstIgram(acquisitionDates, safe_dict, pairs) + runObj.finalize() + + i += 1 + runObj = run() + runObj.configure(inps, 'run_' + str(i) + "_merge_burst_igram") + runObj.igram_mergeBurst(acquisitionDates, safe_dict, pairs) runObj.finalize() i+=1 diff --git a/contrib/stack/topsStack/unwrap.py b/contrib/stack/topsStack/unwrap.py index c297dc8..e6c9182 100755 --- a/contrib/stack/topsStack/unwrap.py +++ b/contrib/stack/topsStack/unwrap.py @@ -38,9 +38,11 @@ from isceobj.Constants import SPEED_OF_LIGHT import argparse import os import pickle -import sys +import gdal +import numpy as np #import shelve import s1a_isce_utils as ut +from isceobj.Util.ImageUtil import ImageLib as IML def createParser(): ''' @@ -70,6 +72,9 @@ def createParser(): parser.add_argument('-m', '--method', dest='method', type=str, default='icu', help='unwrapping method') + parser.add_argument('--rmfilter', action='store_true', default=False, + help='remove the effect of filtering from final unwrapped interferograms') + return parser @@ -276,14 +281,48 @@ def runUnwrapIcu(infile, outfile): #unwImage.finalizeImage() unwImage.renderHdr() + +def remove_filter(intfile, filtfile, unwfile): + + outunw = os.path.abspath(unwfile).split('filt_') + outunw = outunw[0] + outunw[1] + + ds_unw = gdal.Open(unwfile + ".vrt", gdal.GA_ReadOnly) + width = ds_unw.RasterXSize + length = ds_unw.RasterYSize + + unw_phase = np.memmap(unwfile, dtype='f', mode='r', shape=(2, length, width)) + filt_phase = np.memmap(filtfile, dtype='f', mode='r', shape=(length, width)) + int_phase = np.memmap(intfile, dtype='f', mode='r', shape=(length, width)) + + outfile = np.memmap(outunw, dtype='f', mode='w+', shape=(2, length, width)) + + for line in range(length): + integer_jumps = unw_phase[1, line, :] - np.angle(filt_phase[line, :]) + outfile[1, line, :] = integer_jumps + np.angle(int_phase[line, :]) + outfile[0, line, :] = unw_phase[0, line, :] + + unwImage = isceobj.Image.createImage() + unwImage.setFilename(outunw) + unwImage.setWidth(width) + unwImage.setLength(length) + unwImage.imageType = 'unw' + unwImage.bands = 2 + unwImage.scheme = 'BIL' + unwImage.dataType = 'FLOAT' + unwImage.setAccessMode('read') + unwImage.renderHdr() + + return + def main(iargs=None): ''' The main driver. ''' inps = cmdLineParse(iargs) - interferogramDir = os.path.dirname(inps.intfile) print ('unwrapping method : ' , inps.method) + if inps.method == 'snaphu': if inps.nomcf: fncall = runUnwrap @@ -298,6 +337,13 @@ def main(iargs=None): elif inps.method == 'icu': runUnwrapIcu(inps.intfile, inps.unwfile) + if inps.rmfilter: + filtfile = os.path.abspath(inps.intfile) + intfile = filtfile.split('filt_') + intfile = intfile[0] + intfile[1] + + remove_filter(intfile, filtfile, inps.unwfile) + if __name__ == '__main__': diff --git a/examples/input_files/alos2/alos2App.xml b/examples/input_files/alos2/alos2App.xml index 8674070..69c86af 100644 --- a/examples/input_files/alos2/alos2App.xml +++ b/examples/input_files/alos2/alos2App.xml @@ -120,14 +120,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/alos2burstApp.xml b/examples/input_files/alos2/alos2burstApp.xml index b95c7cb..18622ca 100644 --- a/examples/input_files/alos2/alos2burstApp.xml +++ b/examples/input_files/alos2/alos2burstApp.xml @@ -120,12 +120,13 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/scansar-scansar/1/alos2App.xml b/examples/input_files/alos2/example_input_files/scansar-scansar/1/alos2App.xml index df2ad17..c597ce2 100644 --- a/examples/input_files/alos2/example_input_files/scansar-scansar/1/alos2App.xml +++ b/examples/input_files/alos2/example_input_files/scansar-scansar/1/alos2App.xml @@ -107,14 +107,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/scansar-scansar/2/alos2App.xml b/examples/input_files/alos2/example_input_files/scansar-scansar/2/alos2App.xml index 01a4f6e..e5257ce 100644 --- a/examples/input_files/alos2/example_input_files/scansar-scansar/2/alos2App.xml +++ b/examples/input_files/alos2/example_input_files/scansar-scansar/2/alos2App.xml @@ -107,14 +107,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/scansar-scansar/3/alos2App.xml b/examples/input_files/alos2/example_input_files/scansar-scansar/3/alos2App.xml index a13212e..03fa5a7 100644 --- a/examples/input_files/alos2/example_input_files/scansar-scansar/3/alos2App.xml +++ b/examples/input_files/alos2/example_input_files/scansar-scansar/3/alos2App.xml @@ -107,14 +107,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/scansar-scansar/4/alos2App.xml b/examples/input_files/alos2/example_input_files/scansar-scansar/4/alos2App.xml index d0139df..6ddc679 100644 --- a/examples/input_files/alos2/example_input_files/scansar-scansar/4/alos2App.xml +++ b/examples/input_files/alos2/example_input_files/scansar-scansar/4/alos2App.xml @@ -107,14 +107,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/scansar-scansar_7s/alos2App.xml b/examples/input_files/alos2/example_input_files/scansar-scansar_7s/alos2App.xml index 5ddcf32..b412162 100644 --- a/examples/input_files/alos2/example_input_files/scansar-scansar_7s/alos2App.xml +++ b/examples/input_files/alos2/example_input_files/scansar-scansar_7s/alos2App.xml @@ -107,14 +107,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/scansar-scansar_burst/1/alos2burstApp.xml b/examples/input_files/alos2/example_input_files/scansar-scansar_burst/1/alos2burstApp.xml index bf16d48..0e7facf 100644 --- a/examples/input_files/alos2/example_input_files/scansar-scansar_burst/1/alos2burstApp.xml +++ b/examples/input_files/alos2/example_input_files/scansar-scansar_burst/1/alos2burstApp.xml @@ -106,12 +106,13 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/scansar-scansar_burst/2/alos2burstApp.xml b/examples/input_files/alos2/example_input_files/scansar-scansar_burst/2/alos2burstApp.xml index 55b8b39..73e7074 100644 --- a/examples/input_files/alos2/example_input_files/scansar-scansar_burst/2/alos2burstApp.xml +++ b/examples/input_files/alos2/example_input_files/scansar-scansar_burst/2/alos2burstApp.xml @@ -106,12 +106,13 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/scansar-scansar_burst/3/alos2burstApp.xml b/examples/input_files/alos2/example_input_files/scansar-scansar_burst/3/alos2burstApp.xml index d89f4f7..39789d1 100644 --- a/examples/input_files/alos2/example_input_files/scansar-scansar_burst/3/alos2burstApp.xml +++ b/examples/input_files/alos2/example_input_files/scansar-scansar_burst/3/alos2burstApp.xml @@ -106,12 +106,13 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/scansar-scansar_burst/4/alos2burstApp.xml b/examples/input_files/alos2/example_input_files/scansar-scansar_burst/4/alos2burstApp.xml index 6a43c82..c31747a 100644 --- a/examples/input_files/alos2/example_input_files/scansar-scansar_burst/4/alos2burstApp.xml +++ b/examples/input_files/alos2/example_input_files/scansar-scansar_burst/4/alos2burstApp.xml @@ -106,12 +106,13 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/scansar-stripmap/1/alos2App.xml b/examples/input_files/alos2/example_input_files/scansar-stripmap/1/alos2App.xml index e85866b..c483d1d 100644 --- a/examples/input_files/alos2/example_input_files/scansar-stripmap/1/alos2App.xml +++ b/examples/input_files/alos2/example_input_files/scansar-stripmap/1/alos2App.xml @@ -105,14 +105,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/scansar-stripmap/2/alos2App.xml b/examples/input_files/alos2/example_input_files/scansar-stripmap/2/alos2App.xml index d0b712b..b8723d2 100644 --- a/examples/input_files/alos2/example_input_files/scansar-stripmap/2/alos2App.xml +++ b/examples/input_files/alos2/example_input_files/scansar-stripmap/2/alos2App.xml @@ -105,14 +105,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/stripmap-stripmap/1/alos2App.xml b/examples/input_files/alos2/example_input_files/stripmap-stripmap/1/alos2App.xml index 3980014..16f7fcb 100644 --- a/examples/input_files/alos2/example_input_files/stripmap-stripmap/1/alos2App.xml +++ b/examples/input_files/alos2/example_input_files/stripmap-stripmap/1/alos2App.xml @@ -105,14 +105,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/stripmap-stripmap/2/alos2App.xml b/examples/input_files/alos2/example_input_files/stripmap-stripmap/2/alos2App.xml index af8e9df..9f1f628 100644 --- a/examples/input_files/alos2/example_input_files/stripmap-stripmap/2/alos2App.xml +++ b/examples/input_files/alos2/example_input_files/stripmap-stripmap/2/alos2App.xml @@ -107,14 +107,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/stripmap-stripmap/3/alos2App.xml b/examples/input_files/alos2/example_input_files/stripmap-stripmap/3/alos2App.xml index cc617e7..bd4d4c9 100644 --- a/examples/input_files/alos2/example_input_files/stripmap-stripmap/3/alos2App.xml +++ b/examples/input_files/alos2/example_input_files/stripmap-stripmap/3/alos2App.xml @@ -105,14 +105,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 diff --git a/examples/input_files/alos2/example_input_files/stripmap-stripmap/4/alos2App.xml b/examples/input_files/alos2/example_input_files/stripmap-stripmap/4/alos2App.xml index 0edca0d..66ca33b 100644 --- a/examples/input_files/alos2/example_input_files/stripmap-stripmap/4/alos2App.xml +++ b/examples/input_files/alos2/example_input_files/stripmap-stripmap/4/alos2App.xml @@ -105,14 +105,15 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450