From 959a278e477c42fc8dfd138a7b86aaca24ea1615 Mon Sep 17 00:00:00 2001 From: piyushrpt Date: Mon, 11 Nov 2019 10:42:11 -0800 Subject: [PATCH 1/5] Look for autorift and install if found --- contrib/SConscript | 3 +++ 1 file changed, 3 insertions(+) diff --git a/contrib/SConscript b/contrib/SConscript index fb91670..ecc0baf 100644 --- a/contrib/SConscript +++ b/contrib/SConscript @@ -78,3 +78,6 @@ SConscript(rfi) SConscript('PyCuAmpcor/SConscript') SConscript('splitSpectrum/SConscript') SConscript('alos2proc/SConscript') + +if os.path.exists('geo_autoRIFT'): + SConscript('geo_autoRIFT/SConscript') From 37a510fd2d359dba2cf7a3ce1ad6c5a4b1871932 Mon Sep 17 00:00:00 2001 From: vbrancat Date: Tue, 12 Nov 2019 10:28:40 -0800 Subject: [PATCH 2/5] add rubbersheetRange and rubbersheetAzimuth in Factories.py --- components/isceobj/StripmapProc/Factories.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/components/isceobj/StripmapProc/Factories.py b/components/isceobj/StripmapProc/Factories.py index 7727f48..27497d8 100644 --- a/components/isceobj/StripmapProc/Factories.py +++ b/components/isceobj/StripmapProc/Factories.py @@ -112,7 +112,8 @@ createResampleSlc = _factory("runResampleSlc") createResampleSubbandSlc = _factory("runResampleSubbandSlc") createRefineSlaveTiming = _factory("runRefineSlaveTiming") createDenseOffsets = _factory("runDenseOffsets") -createRubbersheet = _factory("runRubbersheet") +createRubbersheetAzimuth = _factory("runRubbersheetAzimuth") # Modified by V. Brancato (10.07.2019) +createRubbersheetRange = _factory("runRubbersheetRange") # Modified by V. Brancato (10.07.2019) createInterferogram = _factory("runInterferogram") createCoherence = _factory("runCoherence") createFilter = _factory("runFilter") From 532cfa085b49a7c22a4723f3aebe38eb2db3fd18 Mon Sep 17 00:00:00 2001 From: vbrancat Date: Wed, 13 Nov 2019 10:33:46 -0800 Subject: [PATCH 3/5] Fixed some bugs in rubbersheeting and runInterferogram --- .../isceobj/StripmapProc/runInterferogram.py | 26 ++++--- .../StripmapProc/runResampleSubbandSlc.py | 2 + .../isceobj/StripmapProc/runRubbersheet.py | 1 + .../StripmapProc/runRubbersheetAzimuth.py | 55 ++++++++++----- .../StripmapProc/runRubbersheetRange.py | 68 ++++++++++++------- 5 files changed, 101 insertions(+), 51 deletions(-) diff --git a/components/isceobj/StripmapProc/runInterferogram.py b/components/isceobj/StripmapProc/runInterferogram.py index 71e3ddd..fbbde4f 100644 --- a/components/isceobj/StripmapProc/runInterferogram.py +++ b/components/isceobj/StripmapProc/runInterferogram.py @@ -31,14 +31,14 @@ def write_xml(fileName,width,length,bands,dataType,scheme): return None -def compute_FlatEarth(self,width,length): +def compute_FlatEarth(self,ifgFilename,width,length,radarWavelength): from imageMath import IML import logging # If rubbersheeting has been performed add back the range sheet offsets info = self._insar.loadProduct(self._insar.slaveSlcCropProduct) - radarWavelength = info.getInstrument().getRadarWavelength() + #radarWavelength = info.getInstrument().getRadarWavelength() rangePixelSize = info.getInstrument().getRangePixelSize() fact = 4 * np.pi* rangePixelSize / radarWavelength @@ -55,7 +55,7 @@ def compute_FlatEarth(self,width,length): rng2 = np.zeros((length,width)) # Open the interferogram - ifgFilename= os.path.join(self.insar.ifgDirname, self.insar.ifgFilename) + #ifgFilename= os.path.join(self.insar.ifgDirname, self.insar.ifgFilename) intf = np.memmap(ifgFilename+'.full',dtype=np.complex64,mode='r+',shape=(length,width)) for ll in range(length): @@ -126,7 +126,8 @@ def computeCoherence(slc1name, slc2name, corname, virtual=True): return # Modified by V. Brancato on 10.09.2019 (added self) -def generateIgram(self,imageSlc1, imageSlc2, resampName, azLooks, rgLooks): +# Modified by V. Brancato on 11.13.2019 (added radar wavelength for low and high band flattening +def generateIgram(self,imageSlc1, imageSlc2, resampName, azLooks, rgLooks,radarWavelength): objSlc1 = isceobj.createSlcImage() IU.copyAttributes(imageSlc1, objSlc1) objSlc1.setAccessMode('read') @@ -192,7 +193,7 @@ def generateIgram(self,imageSlc1, imageSlc2, resampName, azLooks, rgLooks): objCrossmul.crossmul(objSlc1, objSlc2, objInt, objAmp) # Remove Flat-Earth component - compute_FlatEarth(self,intWidth,lines) + compute_FlatEarth(self,resampInt,intWidth,lines,radarWavelength) # Perform Multilook multilook(resampInt+'.full', outname=resampInt, alks=azLooks, rlks=rgLooks) #takeLooks(objAmp,azLooks,rgLooks) @@ -206,7 +207,7 @@ def generateIgram(self,imageSlc1, imageSlc2, resampName, azLooks, rgLooks): return imageInt, imageAmp -def subBandIgram(self, masterSlc, slaveSlc, subBandDir): +def subBandIgram(self, masterSlc, slaveSlc, subBandDir,radarWavelength): img1 = isceobj.createImage() img1.load(masterSlc + '.xml') @@ -226,7 +227,7 @@ def subBandIgram(self, masterSlc, slaveSlc, subBandDir): interferogramName = os.path.join(ifgDir , self.insar.ifgFilename) - generateIgram(self,img1, img2, interferogramName, azLooks, rgLooks) + generateIgram(self,img1, img2, interferogramName, azLooks, rgLooks,radarWavelength) return interferogramName @@ -259,9 +260,9 @@ def runSubBandInterferograms(self): slaveHighBandSlc = os.path.join(coregDir , os.path.basename(slaveSlc)) ########## - interferogramName = subBandIgram(self, masterLowBandSlc, slaveLowBandSlc, self.insar.lowBandSlcDirname) + interferogramName = subBandIgram(self, masterLowBandSlc, slaveLowBandSlc, self.insar.lowBandSlcDirname,self.insar.lowBandRadarWavelength) - interferogramName = subBandIgram(self, masterHighBandSlc, slaveHighBandSlc, self.insar.highBandSlcDirname) + interferogramName = subBandIgram(self, masterHighBandSlc, slaveHighBandSlc, self.insar.highBandSlcDirname,self.insar.highBandRadarWavelength) def runFullBandInterferogram(self): logger.info("Generating interferogram") @@ -295,8 +296,11 @@ def runFullBandInterferogram(self): os.makedirs(ifgDir) interferogramName = os.path.join(ifgDir , self.insar.ifgFilename) - - generateIgram(self,img1, img2, interferogramName, azLooks, rgLooks) + + info = self._insar.loadProduct(self._insar.slaveSlcCropProduct) + radarWavelength = info.getInstrument().getRadarWavelength() + + generateIgram(self,img1, img2, interferogramName, azLooks, rgLooks,radarWavelength) ###Compute coherence diff --git a/components/isceobj/StripmapProc/runResampleSubbandSlc.py b/components/isceobj/StripmapProc/runResampleSubbandSlc.py index 2cd7a12..157764b 100644 --- a/components/isceobj/StripmapProc/runResampleSubbandSlc.py +++ b/components/isceobj/StripmapProc/runResampleSubbandSlc.py @@ -61,9 +61,11 @@ def resampleSlc(self,masterFrame, slaveFrame, imageSlc2, radarWavelength, coregD if not self.doRubbersheetingRange: print('Rubber sheeting in range is turned off, flattening the interferogram during resampling') flatten = True + print(flatten) else: print('Rubber sheeting in range is turned on, flattening the interferogram during interferogram formation') flatten=False + print(flatten) # end of Modification rObj.flatten = flatten diff --git a/components/isceobj/StripmapProc/runRubbersheet.py b/components/isceobj/StripmapProc/runRubbersheet.py index 3e10381..b16a58e 100644 --- a/components/isceobj/StripmapProc/runRubbersheet.py +++ b/components/isceobj/StripmapProc/runRubbersheet.py @@ -168,6 +168,7 @@ def runRubbersheet(self): # filtAzOffsetFile to it. resampleOffset(filtAzOffsetFile, geometryAzimuthOffset, sheetOffset) + print("I'm here") return None diff --git a/components/isceobj/StripmapProc/runRubbersheetAzimuth.py b/components/isceobj/StripmapProc/runRubbersheetAzimuth.py index 5035de3..75583f1 100755 --- a/components/isceobj/StripmapProc/runRubbersheetAzimuth.py +++ b/components/isceobj/StripmapProc/runRubbersheetAzimuth.py @@ -2,6 +2,9 @@ # Author: Heresh Fattahi # Copyright 2017 # +# Modified by V. Brancato +# Included offset filtering with no SNR +# import isce import isceobj @@ -20,28 +23,36 @@ def mask_filterNoSNR(denseOffsetFile,filterSize,outName): off_rg = ds.GetRasterBand(2).ReadAsArray() ds = None - off_az_copy = np.copy(off_az) - off_rg_copy = np.copy(off_rg) + # Remove missing values from ampcor + off_rg[np.where(off_rg < -9999)]=0 + off_az[np.where(off_az < -9999)]=0 - # Compute MAD - off_azm = ndimage.median_filter(off_az,filterSize) - off_rgm = ndimage.median_filter(off_rg,filterSize) - metric_rg = np.abs(off_rgm-off_rg) - metric_az = np.abs(off_azm-off_az) + # Store the offsets in a complex variable + off = off_rg + 1j*off_az - idx = np.where((metric_rg > 3) | (metric_az > 3)) + # Mask the azimuth offsets based on the MAD + mask = off_masking(off,filterSize,thre=3) + xoff_masked = np.ma.array(off.real,mask=mask) + yoff_masked = np.ma.array(off.imag,mask=mask) - # Remove missing values in the offset map - off_az_copy[np.where(off_az_copy<-9999)]=np.nan - off_az_copy[idx]=np.nan + # Delete unused variables + mask = None + off = None + # Remove residual noisy spots with a median filter on the azimuth offmap + yoff_masked.mask = yoff_masked.mask | \ + (ndimage.median_filter(xoff_masked.filled(fill_value=0),3) == 0) | \ + (ndimage.median_filter(yoff_masked.filled(fill_value=0),3) == 0) - # Fill the offsets with smoothed values - off_az_filled = fill_with_smoothed(off_az_copy,filterSize) + # Fill the data by iteratively using smoothed values + data = yoff_masked.data + data[yoff_masked.mask]=np.nan - # Median filter the offsets + off_az_filled = fill_with_smoothed(data,filterSize) + + # Apply median filter to smooth the azimuth offset map off_az_filled = ndimage.median_filter(off_az_filled,filterSize) # Save the filtered offsets @@ -63,9 +74,19 @@ def mask_filterNoSNR(denseOffsetFile,filterSize,outName): img.scheme = 'BIP' img.renderHdr() - return None + + return +def off_masking(off,filterSize,thre=2): + # Define the mask to fill the offsets + vram = ndimage.median_filter(off.real, filterSize) + vazm = ndimage.median_filter(off.imag, filterSize) + + mask = (np.abs(off.real-vram) > thre) | (np.abs(off.imag-vazm) > thre) | (off.imag == 0) | (off.real == 0) + + return mask + def fill(data, invalid=None): """ Replace the value of invalid 'data' cells (indicated by 'invalid') @@ -139,11 +160,11 @@ def fill_with_smoothed(off,filterSize): loop = 0 cnt2=1 - while (loop<100): + while (cnt2!=0 & loop<100): loop += 1 idx2= np.isnan(off_2filt) cnt2 = np.sum(np.count_nonzero(np.isnan(off_2filt))) - + print(cnt2) if cnt2 != 0: off_filt= convolve(off_2filt,kernel,boundary='extend',nan_treatment='interpolate') off_2filt[idx2]=off_filt[idx2] diff --git a/components/isceobj/StripmapProc/runRubbersheetRange.py b/components/isceobj/StripmapProc/runRubbersheetRange.py index 42b4cef..2f1ab3f 100755 --- a/components/isceobj/StripmapProc/runRubbersheetRange.py +++ b/components/isceobj/StripmapProc/runRubbersheetRange.py @@ -10,9 +10,10 @@ import isce import isceobj from osgeo import gdal from scipy import ndimage -from astropy.convolution import convolve import numpy as np import os +from astropy.convolution import convolve + def mask_filterNoSNR(denseOffsetFile,filterSize,outName): # Masking the offsets with a data-based approach @@ -23,27 +24,36 @@ def mask_filterNoSNR(denseOffsetFile,filterSize,outName): off_rg = ds.GetRasterBand(2).ReadAsArray() ds = None - off_az_copy = np.copy(off_az) - off_rg_copy = np.copy(off_rg) - - # Compute MAD - off_azm = ndimage.median_filter(off_az,filterSize) - off_rgm = ndimage.median_filter(off_rg,filterSize) - - metric_rg = np.abs(off_rgm-off_rg) - metric_az = np.abs(off_azm-off_az) - - - idx = np.where((metric_rg > 3) | (metric_az > 3)) - - # Remove missing data in the offset map - off_rg_copy[np.where(off_rg_copy<-9999)]=np.nan - off_rg_copy[idx]=np.nan + # Remove values reported as missing data (no value data from ampcor) + off_rg[np.where(off_rg < -9999)]=0 + off_az[np.where(off_az < -9999)]=0 + + # Store the offsets in a complex variable + off = off_rg + 1j*off_az - # Fill the offsets with smoothed values - off_rg_filled = fill_with_smoothed(off_rg_copy,filterSize) + # Mask the offset based on MAD + mask = off_masking(off,filterSize,thre=3) - # Median filter the offsets + xoff_masked = np.ma.array(off.real,mask=mask) + yoff_masked = np.ma.array(off.imag,mask=mask) + + # Delete not used variables + mask = None + off = None + + # Remove residual noisy spots with a median filter on the range offmap + xoff_masked.mask = xoff_masked.mask | \ + (ndimage.median_filter(xoff_masked.filled(fill_value=0),3) == 0) | \ + (ndimage.median_filter(yoff_masked.filled(fill_value=0),3) == 0) + + # Fill the range offset map iteratively with smoothed values + + data = xoff_masked.data + data[xoff_masked.mask]=np.nan + + off_rg_filled = fill_with_smoothed(data,filterSize) + + # Apply the median filter on the offset off_rg_filled = ndimage.median_filter(off_rg_filled,filterSize) # Save the filtered offsets @@ -64,8 +74,17 @@ def mask_filterNoSNR(denseOffsetFile,filterSize,outName): img.dataType = 'FLOAT' img.scheme = 'BIP' img.renderHdr() + + return + +def off_masking(off,filterSize,thre=2): + vram = ndimage.median_filter(off.real, filterSize) + vazm = ndimage.median_filter(off.imag, filterSize) + + mask = (np.abs(off.real-vram) > thre) | (np.abs(off.imag-vazm) > thre) | (off.imag == 0) | (off.real == 0) + + return mask - return None def fill(data, invalid=None): """ @@ -95,11 +114,11 @@ def fill_with_smoothed(off,filterSize): loop = 0 cnt2=1 - while (loop<100): + while (cnt2 !=0 & loop<100): loop += 1 idx2= np.isnan(off_2filt) cnt2 = np.sum(np.count_nonzero(np.isnan(off_2filt))) - + print(cnt2) if cnt2 != 0: off_filt= convolve(off_2filt,kernel,boundary='extend',nan_treatment='interpolate') off_2filt[idx2]=off_filt[idx2] @@ -255,3 +274,6 @@ def runRubbersheetRange(self): + + + From 5078c3ec8c5abbe04b626eff2ae44e30386046bf Mon Sep 17 00:00:00 2001 From: vbrancat Date: Fri, 15 Nov 2019 20:01:56 -0800 Subject: [PATCH 4/5] Add RubbersheetRange.py and runRubbersheetAzimuth.py in SConscript --- components/isceobj/StripmapProc/SConscript | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/isceobj/StripmapProc/SConscript b/components/isceobj/StripmapProc/SConscript index 2926124..89adcb2 100644 --- a/components/isceobj/StripmapProc/SConscript +++ b/components/isceobj/StripmapProc/SConscript @@ -49,7 +49,7 @@ listFiles = ['StripmapProc.py', 'runPreprocessor.py', 'runSplitSpectrum.py', 'Factories.py' , 'runDenseOffsets.py', 'runResampleSlc.py' , 'runUnwrapGrass.py', '__init__.py' , 'runDispersive.py' , 'runResampleSubbandSlc.py', 'runUnwrapIcu.py', 'runFilter.py' , 'runROI.py' , 'runUnwrapSnaphu.py', 'runCrop.py', - 'runGeo2rdr.py', 'runRubbersheet.py', '__StripmapProc.py' , 'runInterferogram.py', + 'runGeo2rdr.py', 'runRubbersheetRange.py', 'runRubbersheetAzimuth.py', '__StripmapProc.py' , 'runInterferogram.py', 'runVerifyDEM.py', 'runGeocode.py', 'Sensor.py' ] From cf66fb07898863bdf81fe3d59c0591fb667f5a73 Mon Sep 17 00:00:00 2001 From: vbrancat Date: Mon, 18 Nov 2019 08:29:50 -0800 Subject: [PATCH 5/5] added astropy to conda installation --- .circleci/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index a3aded6..98bb5bd 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -23,7 +23,7 @@ jobs: pwd mkdir config build install . /opt/conda/bin/activate root - conda install --yes cython gdal h5py libgdal pytest numpy fftw scipy basemap scons opencv hdf4 hdf5 netcdf4 libgcc libstdcxx-ng cmake + conda install --yes cython gdal h5py libgdal pytest numpy fftw scipy basemap scons opencv hdf4 hdf5 netcdf4 libgcc libstdcxx-ng cmake astropy yum install -y uuid-devel x11-devel motif-devel jq gcc-gfortran ln -s /opt/conda/bin/cython /opt/conda/bin/cython3 cd /opt/conda/lib