fix merge

.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

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")

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'
 ]
 

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

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

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
 
 

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)
+    # Remove missing values from ampcor 
-    off_rg_copy = np.copy(off_rg)
+    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)
+    # Store the offsets in a complex variable
-    metric_az = np.abs(off_azm-off_az)
+    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
+    # Delete unused variables
-    off_az_copy[np.where(off_az_copy<-9999)]=np.nan
+    mask = None
-    off_az_copy[idx]=np.nan
+    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
+    # Fill the data by iteratively using smoothed values    
-    off_az_filled = fill_with_smoothed(off_az_copy,filterSize)
+    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]

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)
+    # Remove values reported as missing data (no value data from ampcor)
-    off_rg_copy = np.copy(off_rg)
+    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)
+    # Store the offsets in a complex variable
-    off_rgm = ndimage.median_filter(off_rg,filterSize)
+    off = off_rg + 1j*off_az
-    
-    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
 
-    # Fill the offsets with smoothed values
+    # Mask the offset based on MAD
-    off_rg_filled = fill_with_smoothed(off_rg_copy,filterSize)
+    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):
 
 
 
+
+
+

contrib/SConscript

@@ -79,5 +79,5 @@ SConscript('PyCuAmpcor/SConscript')
 SConscript('splitSpectrum/SConscript')
 SConscript('alos2proc/SConscript')
 
-if os.path.exists('geoAutorift'):
+if os.path.exists('geo_autoRIFT'):
-    SConscript('geoAutorift/SConscript')
+    SConscript('geo_autoRIFT/SConscript')