diff --git a/contrib/stack/stripmapStack/README.md b/contrib/stack/stripmapStack/README.md new file mode 100644 index 0000000..4bf9260 --- /dev/null +++ b/contrib/stack/stripmapStack/README.md @@ -0,0 +1,85 @@ +## StripMap stack processor + +The detailed algorithms and workflow for stack processing of stripmap SAR data can be found here: + ++ Fattahi, H., M. Simons, and P. Agram (2017), InSAR Time-Series Estimation of the Ionospheric Phase Delay: An Extension of the Split Range-Spectrum Technique, IEEE Transactions on Geoscience and Remote Sensing, 55(10), 5984-5996, doi:[10.1109/TGRS.2017.2718566](https://ieeexplore.ieee.org/abstract/document/7987747/). + +----------------------------------- + +To use the stripmap stack processor, make sure to add the path of your `contrib/stack/stripmapStack` folder to your `$PATH` environment varibale. + +Currently supported workflows include a coregistered stack of SLC, interferograms, ionospheric delays. + +Here are some notes to get started with processing stacks of stripmap data with ISCE. + +#### 1. Create your project folder somewhere + +``` +mkdir MauleAlosDT111 +cd MauleAlosDT111 +``` + +#### 2. Prepare DEM +a) create a folder for DEM; +b) create a DEM using dem.py with SNWE of your study area in integer; +d) Keep only ".dem.wgs84", ".dem.wgs84.vrt" and ".dem.wgs84.xml" and remove unnecessary files; +d) fix the path of the file in the xml file of the DEM by using fixImageXml.py. + +``` +mkdir DEM; cd DEM +dem.py -a stitch -b -37 -31 -72 -69 -r -s 1 -c +rm demLat*.dem demLat*.dem.xml demLat*.dem.vrt +fixImageXml.py -f -i demLat*.dem.wgs84 +cd .. +``` + +#### 3. Download data + +##### 3.1 create a folder to download SAR data (i.e. ALOS-1 data from ASF) + +``` +mkdir download +cd download +``` + +##### 3.2 Download the data that that you want to process to the "downlowd" directory + +#### 4. Prepare SAR data + +Once all data have been downloaded, we need to unzip them and move them to different folders and getting ready for unpacking and then SLC generation. This can be done by running the following command in a directory above "download": + +``` +prepRawALOS.py -i download/ -o SLC +``` + +This command generates an empty SLC folder and a run file called: "run_unPackALOS". + +You could also run prepRawSensor.py which aims to recognize the sensor data automatically followed by running the sensor specific preparation script. For now we include support for ALOS and CSK raw data, but it is trivial to expand and include other sensors as unpacking routines are already included in the distribution. + +``` +prepRawSensor.py -i download/ -o SLC +``` + +#### 5. Execute the commands in "run_unPackALOS" file + +If you have a cluster that you can submit jobs, you can submit each line of command to a processor. The commands are independent and can be run in parallel. + +After successfully running the previous step, you should see acquisition dates in the SLC folder and the ".raw" files for each acquisition. + +Note: For ALOS-1, If there is an acquisition that does not include .raw file, this is most likely due to PRF change between frames and can not be currently handled by ISCE. You have to ignore those. + +#### 6. Run "stackStripmap.py" + +This will generate many config and run files that need to be executed. Here is an example: + +``` +stackStripMap.py -s SLC/ -d DEM/demLat*.dem.wgs84 -t 250 -b 1000 -a 14 -r 4 -u snaphu +``` + +This will produce: +a) baseline folder, which contains baseline information +b) pairs.png which is a baseline-time plot of the network of interferograms +c) configs: which contains the configuration parameter to run different InSAR processing steps +d) run_files: a folder that includes several run and job files that needs to be run in order + +#### 7. Execute the commands in run files (run_1*, run_2*, etc) in the "run_files" folder diff --git a/contrib/stack/stripmapStack/README.txt b/contrib/stack/stripmapStack/README.txt deleted file mode 100644 index e6ae915..0000000 --- a/contrib/stack/stripmapStack/README.txt +++ /dev/null @@ -1,64 +0,0 @@ - -The detailed algorithms for stack processing of stripmap data can be found here: - -H. Fattahi, M. Simons, and P. Agram, "InSAR Time-Series Estimation of the Ionospheric Phase Delay: An Extension of the Split Range-Spectrum Technique", IEEE Trans. Geosci. Remote Sens., vol. 55, no. 10, 5984-5996, 2017. (https://ieeexplore.ieee.org/abstract/document/7987747/) - - ------------------------------------ - -Notes on stripmap stack processor: - -Here are some notes to get started with processing stacks of stripmap data with ISCE. - - -1- create a folder somewhere for your project - -mkdir MauleT111 -cd MauleT111 - -2- create a DEM: - -dem.py -a stitch -b -37 -31 -72 -69 -r -s 1 -c - -3- Keep only ".dem.wgs84", ".dem.wgs84.vrt" and ".dem.wgs84.xml" and remove unnecessary files - -4- fix the path of the file in the xml file of the DEM by using this command: - -fixImageXml.py -f -i demLat_S37_S31_Lon_W072_W069.dem.wgs84 - -5- create a folder to download the ALOS-1 data from ASF: - -mkdir download -cd download - -6- Download the data that that you want to process to the downlowd directory. - -7- once all data have been downloaded, we need to unzip them and move them to different folders and getting ready for unpacking and then SLC generation. -This can be done by running the following command in a directory above "download": - -prepRawALOS.py -i download/ -o SLC - -This command generates an empty SLC folder and a run file called: "run_unPackALOS". -You could also run prepRawSensor.py which aims to recognize the sensor data automatically followed by running the sensor specific preparation script. For now we include support for ALOS and CSK raw data, but it is trivial to expand and include other sensors as unpacking routines are already included in the distribution. - -prepRawSensor.py -i download/ -o SLC - -8- execute the commands inside run_unPackALOS file. If you have a cluster that you can submit jobs, you can submit each line of command to a processor. The commands are independent and can be run in parallel. - -9- After successfully running the previous step, you should see acquisition dates in the SLC folder and the ".raw" files for each acquisition - -Note: For ALOS-1, If there is an acquisition that does not include .raw file, this is most likely due to PRF change between frames and can not be currently handled by ISCE. You have to ignore those. - -10- run stackStripmap.py which will generate many config and run files that need to be executed. Here is an example: - -stackStripMap.py -s SLC/ -d demLat_S37_S31_Lon_W072_W069.dem.wgs84 -t 250 -b 1000 -a 14 -r 4 -u snaphu - -This will produce: -a) baseline folder, which contains baseline information -b) pairs.png which is a baseline-time plot of the network of interferograms -c) configs: which contains the configuration parameter to run different InSAR processing steps -d) run_files: a folder that includes several run and job files that needs to be run in order - -11- execute the commands in run files (run_1, run_2, etc) in the run_files folder - - diff --git a/contrib/stack/stripmapStack/stackStripMap.py b/contrib/stack/stripmapStack/stackStripMap.py index 07c1f68..728ffe0 100755 --- a/contrib/stack/stripmapStack/stackStripMap.py +++ b/contrib/stack/stripmapStack/stackStripMap.py @@ -5,7 +5,7 @@ import os, imp, sys, glob import argparse import configparser -import datetime +import datetime import numpy as np import shelve import isce @@ -20,7 +20,7 @@ defoMax = '2' maxNodes = 72 def createParser(): - parser = argparse.ArgumentParser( description='Preparing the directory structure and config files for stack processing of Sentinel data') + parser = argparse.ArgumentParser( description='Preparing the directory structure and config files for stack processing of StripMap data') parser.add_argument('-s', '--slc_directory', dest='slcDir', type=str, required=True, help='Directory with all stripmap SLCs') @@ -31,7 +31,7 @@ def createParser(): help='Working directory ') parser.add_argument('-d', '--dem', dest='dem', type=str, required=True, - help='Directory with the DEM (.xml and .vrt files)') + help='DEM file (with .xml and .vrt files)') parser.add_argument('-m', '--master_date', dest='masterDate', type=str, default=None, help='Directory with master acquisition') @@ -43,47 +43,54 @@ def createParser(): help='Baseline threshold (max bperp in meters)') parser.add_argument('-a', '--azimuth_looks', dest='alks', type=str, default='10', - help='Number of looks in azimuth (automaticly computed as AspectR*looks when "S" or "sensor" is defined to give approximately square multi-look pixels)') + help='Number of looks in azimuth (automaticly computed as AspectR*looks when ' + '"S" or "sensor" is defined to give approximately square multi-look pixels)') parser.add_argument('-r', '--range_looks', dest='rlks', type=str, default='10', help='Number of looks in range') parser.add_argument('-S', '--sensor', dest='sensor', type=str, required=False, help='SAR sensor used to define square multi-look pixels') - parser.add_argument('-L', '--low_band_frequency', dest='fL', type=str, default=None, - help='low band frequency') - parser.add_argument('-H', '--high_band_frequency', dest='fH', type=str, default=None, - help='high band frequency') - parser.add_argument('-B', '--subband_bandwidth ', dest='bandWidth', type=str, default=None, - help='sub-band band width') - parser.add_argument('-u', '--unw_method', dest='unwMethod', type=str, default='snaphu' - , help='unwrapping method (icu, snaphu, or snaphu2stage)') + + parser.add_argument('-u', '--unw_method', dest='unwMethod', type=str, default='snaphu', + help='unwrapping method (icu, snaphu, or snaphu2stage)') parser.add_argument('-f','--filter_strength', dest='filtStrength', type=str, default=filtStrength, help='strength of Goldstein filter applied to the wrapped phase before spatial coherence estimation.' ' Default: {}'.format(filtStrength)) - parser.add_argument('--filter_sigma_x', dest='filterSigmaX', type=str, default='100' - , help='filter sigma for gaussian filtering the dispersive and nonDispersive phase') + iono = parser.add_argument_group('Ionosphere', 'Configurationas for ionospheric correction') + iono.add_argument('-L', '--low_band_frequency', dest='fL', type=str, default=None, + help='low band frequency') + iono.add_argument('-H', '--high_band_frequency', dest='fH', type=str, default=None, + help='high band frequency') + iono.add_argument('-B', '--subband_bandwidth ', dest='bandWidth', type=str, default=None, + help='sub-band band width') - parser.add_argument('--filter_sigma_y', dest='filterSigmaY', type=str, default='100.0', - help='sigma of the gaussian filter in Y direction, default=100') + iono.add_argument('--filter_sigma_x', dest='filterSigmaX', type=str, default='100', + help='filter sigma for gaussian filtering the dispersive and nonDispersive phase') - parser.add_argument('--filter_size_x', dest='filterSizeX', type=str, default='800.0', - help='size of the gaussian kernel in X direction, default = 800') + iono.add_argument('--filter_sigma_y', dest='filterSigmaY', type=str, default='100.0', + help='sigma of the gaussian filter in Y direction, default=100') - parser.add_argument('--filter_size_y', dest='filterSizeY', type=str, default='800.0', - help='size of the gaussian kernel in Y direction, default=800') + iono.add_argument('--filter_size_x', dest='filterSizeX', type=str, default='800.0', + help='size of the gaussian kernel in X direction, default = 800') - parser.add_argument('--filter_kernel_rotation', dest='filterKernelRotation', type=str, default='0.0', - help='rotation angle of the filter kernel in degrees (default = 0.0)') + iono.add_argument('--filter_size_y', dest='filterSizeY', type=str, default='800.0', + help='size of the gaussian kernel in Y direction, default=800') - parser.add_argument('-W', '--workflow', dest='workflow', type=str, default='slc' - , help='The InSAR processing workflow : (slc, interferogram, ionosphere)') + iono.add_argument('--filter_kernel_rotation', dest='filterKernelRotation', type=str, default='0.0', + help='rotation angle of the filter kernel in degrees (default = 0.0)') - parser.add_argument('-z', '--zero', dest='zerodop', action='store_true', default=False, help='Use zero doppler geometry for processing - Default : No') - parser.add_argument('--nofocus', dest='nofocus', action='store_true', default=False, help='If input data is already focused to SLCs - Default : do focus') - parser.add_argument('-c', '--text_cmd', dest='text_cmd', type=str, default='' - , help='text command to be added to the beginning of each line of the run files. Example : source ~/.bash_profile;') - parser.add_argument('-useGPU', '--useGPU', dest='useGPU',action='store_true', default=False, help='Allow App to use GPU when available') + parser.add_argument('-W', '--workflow', dest='workflow', type=str, default='slc', + help='The InSAR processing workflow : (slc, interferogram, ionosphere)') + + parser.add_argument('-z', '--zero', dest='zerodop', action='store_true', default=False, + help='Use zero doppler geometry for processing - Default : No') + parser.add_argument('--nofocus', dest='nofocus', action='store_true', default=False, + help='If input data is already focused to SLCs - Default : do focus') + parser.add_argument('-c', '--text_cmd', dest='text_cmd', type=str, default='', + help='text command to be added to the beginning of each line of the run files. Example : source ~/.bash_profile;') + parser.add_argument('-useGPU', '--useGPU', dest='useGPU',action='store_true', default=False, + help='Allow App to use GPU when available') parser.add_argument('--summary', dest='summary', action='store_true', default=False, help='Show summary only') return parser diff --git a/contrib/stack/stripmapStack/topo.py b/contrib/stack/stripmapStack/topo.py index b730af5..2b175a6 100755 --- a/contrib/stack/stripmapStack/topo.py +++ b/contrib/stack/stripmapStack/topo.py @@ -1,17 +1,17 @@ #!/usr/bin/env python3 + +import os import argparse +import shelve +import datetime +import shutil +import numpy as np import isce import isceobj -import numpy as np -import shelve -import os -import datetime -import shutil from isceobj.Constants import SPEED_OF_LIGHT from isceobj.Util.Poly2D import Poly2D from mroipac.looks.Looks import Looks - def createParser(): ''' Command line parser. @@ -48,7 +48,7 @@ class Dummy(object): def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): - + from isceobj.Planet.Planet import Planet from zerodop.GPUtopozero.GPUtopozero import PyTopozero from isceobj import Constants as CN @@ -84,14 +84,14 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): omethod = 2 # LEGENDRE INTERPOLATION else: omethod = 0 # HERMITE INTERPOLATION - + # tracking doppler specifications if nativedop and (dop is not None): try: coeffs = dop._coeffs except: coeffs = dop - + polyDoppler = Poly2D() polyDoppler.setWidth(width) polyDoppler.setLength(length) @@ -109,10 +109,10 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): polyDoppler.createPoly2D() - # dem + # dem demImage.setCaster('read','FLOAT') demImage.createImage() - + # slant range file slantRangeImage = Poly2D() slantRangeImage.setWidth(width) @@ -130,12 +130,12 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): dataType = 'DOUBLE' latImage.initImage(latFilename,accessMode,width,dataType) latImage.createImage() - + # lon file lonImage = isceobj.createImage() lonImage.initImage(lonFilename,accessMode,width,dataType) lonImage.createImage() - + # LOS file losImage = isceobj.createImage() dataType = 'FLOAT' @@ -144,7 +144,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): losImage.initImage(losFilename,accessMode,width,dataType,bands=bands,scheme=scheme) losImage.setCaster('write','DOUBLE') losImage.createImage() - + # height file heightImage = isceobj.createImage() dataType = 'DOUBLE' @@ -158,7 +158,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): incImage.initImage(incFilename,accessMode,width,dataType,bands=bands,scheme=scheme) incImage.createImage() incImagePtr = incImage.getImagePointer() - + maskImage = isceobj.createImage() dataType = 'BYTE' bands = 1 @@ -168,7 +168,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): else: incImagePtr = 0 maskImagePtr = 0 - + # initalize planet elp = Planet(pname='Earth').ellipsoid @@ -214,14 +214,14 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): topo.set_orbitBasis(1) # Is this ever different? topo.createOrbit() # Initializes the empty orbit to the right allocated size count = 0 - + for sv in info.orbit.stateVectors.list: td = DTU.seconds_since_midnight(sv.getTime()) pos = sv.getPosition() vel = sv.getVelocity() topo.set_orbitVector(count,td,pos[0],pos[1],pos[2],vel[0],vel[1],vel[2]) count += 1 - + # run topo topo.runTopo() @@ -244,13 +244,13 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): # los file descr = '''Two channel Line-Of-Sight geometry image (all angles in degrees). Represents vector drawn from target to platform. - Channel 1: Incidence angle measured from vertical at target (always +ve). + Channel 1: Incidence angle measured from vertical at target (always +ve). Channel 2: Azimuth angle measured from North in Anti-clockwise direction.''' losImage.setImageType('bil') losImage.addDescription(descr) losImage.finalizeImage() losImage.renderHdr() - + # dem/ height file demImage.finalizeImage() @@ -259,7 +259,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): descr = '''Two channel angle file. Channel 1: Angle between ray to target and the vertical at the sensor Channel 2: Local incidence angle accounting for DEM slope at target''' - + incImage.addDescription(descr) incImage.finalizeImage() incImage.renderHdr() @@ -268,7 +268,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): maskImage.addDescription(descr) maskImage.finalizeImage() maskImage.renderHdr() - + if slantRangeImage: try: slantRangeImage.finalizeImage() @@ -276,7 +276,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): pass -def runTopoCPU(info, demImage, dop=None, +def runTopoCPU(info, demImage, dop=None, nativedop=False, legendre=False): from zerodop.topozero import createTopozero from isceobj.Planet.Planet import Planet @@ -298,7 +298,7 @@ def runTopoCPU(info, demImage, dop=None, topo.numberRangeLooks = info.numberRangeLooks topo.numberAzimuthLooks = info.numberAzimuthLooks topo.lookSide = info.lookSide - topo.sensingStart = info.sensingStart + datetime.timedelta(seconds = ((info.numberAzimuthLooks - 1) /2) / info.prf) + topo.sensingStart = info.sensingStart + datetime.timedelta(seconds = ((info.numberAzimuthLooks - 1) /2) / info.prf) topo.rangeFirstSample = info.rangeFirstSample + ((info.numberRangeLooks - 1)/2) * info.slantRangePixelSpacing topo.demInterpolationMethod='BIQUINTIC' @@ -331,9 +331,10 @@ def runTopoCPU(info, demImage, dop=None, topo.topo() return + def runSimamp(outdir, hname='z.rdr'): from iscesys.StdOEL.StdOELPy import create_writer - + #####Run simamp stdWriter = create_writer("log","",True,filename='sim.log') objShade = isceobj.createSimamplitude() @@ -461,9 +462,9 @@ def extractInfo(frame, inps): info.prf = frame.PRF info.radarWavelength = frame.radarWavelegth info.orbit = frame.getOrbit() - - info.width = frame.getNumberOfSamples() - info.length = frame.getNumberOfLines() + + info.width = frame.getNumberOfSamples() + info.length = frame.getNumberOfLines() info.sensingStop = frame.getSensingStop() info.outdir = inps.outdir @@ -472,7 +473,7 @@ def extractInfo(frame, inps): def main(iargs=None): - + inps = cmdLineParse(iargs) # see if the user compiled isce with GPU enabled @@ -502,11 +503,9 @@ def main(iargs=None): doppler = db['doppler'] except: doppler = frame._dopplerVsPixel - db.close() - ####Setup dem demImage = isceobj.createDemImage() demImage.load(inps.dem + '.xml') @@ -522,7 +521,6 @@ def main(iargs=None): info.incFilename = os.path.join(info.outdir, 'incLocal.rdr') info.maskFilename = os.path.join(info.outdir, 'shadowMask.rdr') - runTopo(info,demImage,dop=doppler,nativedop=inps.nativedop, legendre=inps.legendre) runSimamp(os.path.dirname(info.heightFilename),os.path.basename(info.heightFilename)) @@ -540,4 +538,3 @@ if __name__ == '__main__': Main driver. ''' main() -