From 776f31286818e0a171ebb72c425bc70a0259d332 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 18 Apr 2023 09:48:04 -0700 Subject: [PATCH 01/33] updating iono workflow --- isce2_topsapp/__main__.py | 13 +- isce2_topsapp/iono_proc.py | 652 ++++++++++++++++++++++++++++--------- 2 files changed, 514 insertions(+), 151 deletions(-) diff --git a/isce2_topsapp/__main__.py b/isce2_topsapp/__main__.py index 9b5243c2..c03902ad 100644 --- a/isce2_topsapp/__main__.py +++ b/isce2_topsapp/__main__.py @@ -193,15 +193,16 @@ def gunw_slc(): ) # Run ionospheric correction + # MG: correct burst jumps, e.g. needed for Arabian + # processing. TODO: We need a trigger function for + # this option (it adds almost double time to iono) + # example: look at filt_toposphase.unw or .flat + # and analyze if there are any burst jumps, + # if yes, set this option True if args.estimate_ionosphere_delay: iono_processing( - reference_slc_zips=loc_data["ref_paths"], - secondary_slc_zips=loc_data["sec_paths"], - orbit_directory=loc_data["orbit_directory"], - extent=loc_data["processing_extent"], - dem_for_proc=loc_data["full_res_dem_path"], - dem_for_geoc=loc_data["low_res_dem_path"], mask_filename=loc_data["water_mask"], + correct_burst_jumps=False, ) ref_properties = loc_data["reference_properties"] diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index cef5b0be..fa0848c5 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -1,173 +1,169 @@ +# Author: Marin Govorcin +# Copyright 2023 +# California Institute of Technology + import os import site -import subprocess from pathlib import Path from typing import Union +from skimage import morphology +import scipy.signal as ss +import multiprocessing import numpy as np +from isce.applications import topsApp +from isce.components.isceobj.TopsProc import runIon from isce.components import isceobj from isce.components.isceobj.TopsProc.runMergeBursts import ( - interpolateDifferentNumberOfLooks, + interpolateDifferentNumberOfLooks, multilook, + mergeBursts2, mergeBox ) -from jinja2 import Template + from numpy.typing import NDArray from osgeo import gdal from tqdm import tqdm +''' +TODO: add unwrapping changes to ionosphere swath by swath +''' # List of parameters for ionospheric correction: -# https://github.com/isce-framework/isce2/blob/8af43d2b9f291e2370e9112a3ead1449d0f2a386/examples/input_files/topsApp.xml#L182 +# isce2/examples/input_files/topsApp.xml#L182 -IONO_STEPS = ["subband", "rawion", "grd2ion", "filt_gaussian"] -TEMPLATE_DIR = Path(__file__).parent / "templates" +GEOCODE_LIST_ION = ["merged/topophase.ion"] -GEOCODE_LIST_BASE = ["merged/topophase.ion"] +omp_env = os.getenv('OMP_NUM_THREADS') +# set default number of threads to 4 if env does not exist +if not omp_env: + omp_env = str(4) def iono_processing( *, - reference_slc_zips: list, - secondary_slc_zips: list, - orbit_directory: str, - extent: list, - dem_for_proc: str, - dem_for_geoc: str, - azimuth_looks: int = 7, - range_looks: int = 19, - swaths: list = None, - do_esd: bool = False, - esd_coherence_threshold: float = 0.7, - mask_filename: str = None, -): - swaths = swaths or [1, 2, 3] - # for [ymin, ymax, xmin, xmax] - extent_isce = [extent[k] for k in [1, 3, 0, 2]] + topsapp_xml_filename: str = 'topsApp.xml', + mask_filename: str = '', + correct_burst_jumps: bool = False, + num_threads: str = '4') -> None: + ''' + NOTE: If water mask is not used, the code will return to its + default processing with addition of using briding of unwrapped + phase components, and modifed adaptive gaussian filtering. + Outlier removal and masking using connected component 0 will + be skipped + ''' + + # Update the number of threads + if num_threads == 'all': + # Use all possible threads + num_threads = str(multiprocessing.cpu_count()) + + os.environ["OMP_NUM_THREADS"] = str(num_threads) # Update PATH with ISCE2 applications + # Need for isce2/bin/imageMath.py in runIon.unwrap function isce_app_path = Path(f"{site.getsitepackages()[0]}" "/isce/applications/") os.environ["PATH"] += ":" + str(isce_app_path) - tops_app_cmd = f"{isce_app_path}/topsApp.py" - - # Define topsApp input xml parameters - ion_kwargs = { - "orbit_directory": orbit_directory, - "output_reference_directory": "reference", - "output_secondary_directory": "secondary", - "ref_zip_file": reference_slc_zips, - "sec_zip_file": secondary_slc_zips, - "region_of_interest": extent_isce, - "demFilename": dem_for_proc, - "geocodeDemFilename": dem_for_geoc, - "filter_strength": 0.5, - "do_unwrap": True, - "use_virtual_files": True, - "do_esd": do_esd, - "esd_coherence_threshold": esd_coherence_threshold, - # IONO PARAMETERS - "estimate_ionosphere_delay": True, - "ion_burstProperties": False, - "ion_polyFit": True, - "ion_correctAzimuthshift": 1, - "azimuth_looks": azimuth_looks, - "range_looks": range_looks, - "swaths": swaths, - "geocode_list": GEOCODE_LIST_BASE, - } - - with open(TEMPLATE_DIR / "topsapp_iono_template.xml", "r") as file: - template = Template(file.read()) - - # Step-1 Run ionospheric correction step subband - iono_xml = template.render( - **ion_kwargs, ion_startStep=IONO_STEPS[0], ion_stopStep=IONO_STEPS[0] - ) + # Use Connected component 0 to mask unwrapped interferogram + # to improve masking after using water mask, e.g. remove noisy + # pixels along the coastline + conncomp_flag = True if mask_filename else False - with open("ionoApp.xml", "w") as file: - file.write(iono_xml) + # Use outlier removal only is water mask exist + # if False return to defualt topsProc iono processing + outlier_removal = True if mask_filename else False - step_cmd = f"{tops_app_cmd} ionoApp.xml --dostep=ion" - result = subprocess.run(step_cmd, shell=True) - if result.returncode != 0: - raise ValueError("TopsApp failed at step: ion-subband") + # Outlier removal settings for gaussian filtering + if outlier_removal: + coh_threshold = 0.8 + sigma_rule = 1 # remove value above mean +- 1*std + else: + coh_threshold = 0.0 + sigma_rule = 0 - # Step-2 Mask upper and lower band burst wrapped interferograms - if mask_filename: - # Project mask to burst image coordinate space - workdir = Path(".").resolve() - mask_iono_ifg_bursts(workdir, mask_filename) + # Load input file + topsapp = topsApp.TopsInSAR(name="topsApp", cmdline=[topsapp_xml_filename]) + topsapp.configure() - # Step-3 Compute ionospheric correction - iono_xml = template.render( - **ion_kwargs, ion_startStep=IONO_STEPS[1], ion_stopStep=IONO_STEPS[-1] - ) + # Load topsApp PICKLE object from previous topsApp step: fineresamp + topsapp_dir = Path('.').resolve() + topsapp.loadPickleObj(str(topsapp_dir / 'PICKLE/fineresamp')) - with open("ionoApp.xml", "w") as file: - file.write(iono_xml) + # IONOSPHERE + # Run iono setup + ionParam = runIon.setup(topsapp) - step_cmd = f"{tops_app_cmd} ionoApp.xml --dostep=ion" - result = subprocess.run(step_cmd, shell=True) - if result.returncode != 0: - raise ValueError("TopsApp failed at step: ion-rawion2esd") + # Run iono step subband + # Generate lower and upper band fin interferograms + # NOTE: doing resampling ^ speed with n_threads + # Max threads used in resampling is 8, avoid using + # all threads for resampling as somehow this causes slowdown + runIon.subband(topsapp, ionParam) - # Step-4 mergeBursts - # Create merged/topophase.ion file - merge_bursts(range_looks=range_looks, azimuth_looks=azimuth_looks) + # Mask fine interferograms + if mask_filename: + # Project mask to burst image rdr coordinate space + mask_iono_ifg_bursts(topsapp_dir, mask_filename) - # Step-5 Geocode ionospheric correction outputs - step_cmd = f"{tops_app_cmd} ionoApp.xml --dostep=geocode" - result = subprocess.run(step_cmd, shell=True) - if result.returncode != 0: - raise ValueError("TopsApp failed at step: geocode (ion)") + if ionParam.calIonWithMerged: + # Run iono step rawion : merge + runIon.merge(topsapp, ionParam) + # Run iono step rawion : unwrap + unwrap(topsapp, ionParam, + use_bridging=True, use_conncomp=conncomp_flag) -def merge_bursts( - range_looks: int = 19, - azimuth_looks: int = 7, - ion_rangeLooks: int = 200, - ion_azimuthLooks: int = 50, - mergedir: Union[str, Path] = "./merged", -) -> None: - """ - Reference: - https://github.com/isce-framework/isce2/blob/8af43d2b9f291e2370e9112a3ead1449d0f2a386/components/isceobj/TopsProc/runMergeBursts.py#L851 + # Run iono step rawion : ionosphere + runIon.ionosphere(topsapp, ionParam) + else: + # This mode is used for cross + # Sentinel-1A/B interferogram + # NOTE: need to re-write this step to include + # changes with unwrapping, aka bridging + # and using conncomp=0 to mask noise + runIon.ionSwathBySwath(topsapp, ionParam) + + # Run iono step grd2ion + # Resample ionosphere from ground to iono layer + # NOTE: this is not necessary as filtering is using + # not projected raw ionosphere (on the ground) + # runIon.grd2ion(topsapp, ionParam) + + # Run iono step; filt_gaussian + filt_gaussian(topsapp, ionParam, + coh_threshold=coh_threshold, sigma_rule=sigma_rule) + + # Step when using azimuth shift correction + # between bursts + if correct_burst_jumps: + # ionosphere shift + runIon.ionosphere_shift(topsapp, ionParam) + + # resample from ionospheric layer to ground layer, + # get ionosphere for each burst + runIon.ion2grd(topsapp, ionParam) + + # esd + runIon.esd(topsapp, ionParam) + + # Create merged/topophase.ion file + # using ion/ion_burst/ files + merge_multilook_bursts(topsapp, input_dir='ion/ion_burst', + output_filename='topophase.ion') - """ - mergedIfgname = "topophase.flat" - mergedIonname = "topophase.ion" - ######################################### + else: + # Create merged/topophase.ion file + merge_bursts(input_file="ion/ion_cal/filt.ion", + output_filename="topophase.ion") - ionFilt = "ion/ion_cal/filt.ion" - img = isceobj.createImage() - img.load(ionFilt + ".xml") - ionFiltImage = ( - np.fromfile(ionFilt, dtype=np.float32).reshape( - img.length * 2, img.width) - )[1: img.length * 2: 2, :] - img = isceobj.createImage() - img.load(os.path.join(mergedir, mergedIfgname + ".xml")) + # Geocode ionospheric correction outputs + # This step can be faster ^ with num_of_threads + topsapp.runGeocode(GEOCODE_LIST_ION, topsapp.do_unwrap, + topsapp.geocode_bbox) - # interpolate original - ionFiltImageOut = interpolateDifferentNumberOfLooks( - ionFiltImage, - img.length, - img.width, - range_looks, - azimuth_looks, - ion_rangeLooks, - ion_azimuthLooks, - ) - ionFiltOut = os.path.join(mergedir, mergedIonname) - ionFiltImageOut.astype(np.float32).tofile(ionFiltOut) - - image = isceobj.createImage() - image.setDataType("FLOAT") - image.setFilename(ionFiltOut) - image.extraFilename = ionFiltOut + ".vrt" - image.setWidth(img.width) - image.setLength(img.length) - image.renderHdr() + # Return number of threads to default + os.environ["OMP_NUM_THREADS"] = omp_env def mask_iono_ifg_bursts(tops_dir: Path, @@ -183,7 +179,7 @@ def get_burst(x): # Get all geometry files geom_files = list(tops_dir.glob("geom_reference/IW*/lat*.rdr")) - # NOTE: THis can be easily parallized, skipped for now + # NOTE: This can be easily parallized, skip for now with tqdm(total=len(geom_files)) as pbar: for lat in geom_files: pbar.set_description( @@ -229,7 +225,109 @@ def get_burst(x): pbar.update() +def merge_bursts( + input_file: str = "ion/ion_cal/filt.ion", + output_filename: str = "topophase.ion", + range_looks: int = 19, + azimuth_looks: int = 7, + ion_rangeLooks: int = 200, + ion_azimuthLooks: int = 50, + mergedir: Union[str, Path] = "./merged", +) -> None: + """ + Reference: + isce2/components/isceobj/TopsProc/runMergeBursts.py#L851 + + NOTE: change interpolation method from 'cubic' to 'nearest' in + interpolateDifferentNumberOfLooks for raw ionosphere merge. + Cubic creates artifacts around the edges. + + + """ + mergedIfgname = "topophase.flat" + ######################################### + + img = isceobj.createImage() + img.load(input_file + ".xml") + + ionFiltImage = ( + np.fromfile(input_file, dtype=np.float32).reshape( + img.length * 2, img.width) + )[1: img.length * 2: 2, :] + img = isceobj.createImage() + img.load(os.path.join(mergedir, mergedIfgname + ".xml")) + + # interpolate original + ionFiltImageOut = interpolateDifferentNumberOfLooks( + ionFiltImage, + img.length, + img.width, + range_looks, + azimuth_looks, + ion_rangeLooks, + ion_azimuthLooks, + ) + ionFiltOut = os.path.join(mergedir, output_filename) + ionFiltImageOut.astype(np.float32).tofile(ionFiltOut) + + image = isceobj.createImage() + image.setDataType("FLOAT") + image.setFilename(ionFiltOut) + image.extraFilename = ionFiltOut + ".vrt" + image.setWidth(img.width) + image.setLength(img.length) + image.renderHdr() + + +def merge_multilook_bursts(self, + input_dir: str = 'ion/ion_burst', + output_filename: str = 'topophase.ion', + mergedir: Union[str, Path] = "./merged") -> None: + """ + Reference: + isce2/components/isceobj/TopsProc/runMergeBursts.py#L776 + + """ + + if (self.numberRangeLooks == 1) and (self.numberAzimuthLooks == 1): + suffix = '' + else: + suffix = '.full' + + # Get frames (subswaths) + frames = [] + burstIndex = [] + swathList = self._insar.getValidSwathList(self.swaths) + for swath in swathList: + minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1) + if minBurst == maxBurst: + print('Skipping processing of swath {0}'.format(swath)) + continue + ifg = self._insar.loadProduct(os.path.join( + self._insar.fineIfgDirname, 'IW{0}.xml'.format(swath))) + frames.append(ifg) + burstIndex.append([int(swath), minBurst, maxBurst]) + + # Determine merged size + box = mergeBox(frames) + + # Merge iono products + mergeBursts2(frames, + os.path.join(input_dir, 'IW%d', 'burst_%02d.ion'), + burstIndex, box, + os.path.join(mergedir, output_filename+suffix), + virtual=True, validOnly=True) + + # Create merged/topophase.ion file + multilook(os.path.join(mergedir, output_filename+suffix), + outname=os.path.join(mergedir, output_filename), + alks=self.numberAzimuthLooks, + rlks=self.numberRangeLooks) + + # UTILITIES FOR MASKING + + def raster_geo2radar( rasterFilename: Union[str, Path], latFilename: Union[str, Path], @@ -259,10 +357,10 @@ def raster_geo2radar( """ # Open mask file - mask_ds = gdal.Open(rasterFilename + ".vrt") + mask_ds = gdal.Open(str(rasterFilename) + ".vrt") # Open lon and lat file - lon_ds = gdal.Open(lonFilename + ".vrt") - lat_ds = gdal.Open(latFilename + ".vrt") + lon_ds = gdal.Open(str(lonFilename) + ".vrt") + lat_ds = gdal.Open(str(latFilename) + ".vrt") # Get lon and lat arrays lons = lon_ds.ReadAsArray() @@ -286,10 +384,6 @@ def raster_geo2radar( mask_radarcoord[inboundIndex] = mask_ds.ReadAsArray()[ lineIdx[inboundIndex], sampleIdx[inboundIndex] ] - # close gdal instances - mask_ds = None - lon_ds = None - lat_ds = None # Save in isce format if saveFlag: @@ -297,18 +391,23 @@ def raster_geo2radar( image = isceobj.createImage() image.setDataType("BYTE") image.setFilename(outputFilename) - image.extraFilename = outputFilename + ".vrt" + image.extraFilename = str(outputFilename) + ".vrt" image.setWidth(mask_radarcoord.shape[1]) image.setLength(mask_radarcoord.shape[0]) image.renderHdr() + # close gdal instances + mask_ds = None + lon_ds = None + lat_ds = None + return mask_radarcoord def mask_interferogram( ifgFilename: Union[str, Path], maskArray: NDArray, - outFilename: Union[str, Path] = None, + outFilename: Union[str, Path] = '', ) -> None: """ This routine uses mask np.array in rdr to mask wrapped interferogram @@ -323,12 +422,14 @@ def mask_interferogram( """ # Read interferogram - int_ds = gdal.Open(ifgFilename + ".vrt") + int_ds = gdal.Open(str(ifgFilename) + ".vrt") int_array = int_ds.ReadAsArray() if not np.array_equal(int_array.shape, maskArray.shape): raise ValueError('Mask array dimensions do not match' - ' the interferogram dimensions!') + ' the interferogram dimensions!' + f'mask: {maskArray.shape} vs' + f'ifg: {int_array.shape}') # Mask interferogram int_array.imag[np.bool_(maskArray)] = 0 @@ -337,7 +438,7 @@ def mask_interferogram( # Write masked interferogram if outFilename: driver = gdal.GetDriverByName("ISCE") - outdata = driver.CreateCopy(ifgFilename + "_msk", int_ds) + outdata = driver.CreateCopy(str(ifgFilename) + "_msk", int_ds) outdata.GetRasterBand(1).WriteArray(int_array) outdata.FlushCache() # saves to disk!! outdata = None @@ -345,8 +446,8 @@ def mask_interferogram( # create isce aux files image = isceobj.createIntImage() image.setDataType("cfloat") - image.setFilename(outFilename + "_msk") - image.extraFilename = outFilename + "_msk" + ".vrt" + image.setFilename(str(outFilename) + "_msk") + image.extraFilename = str(outFilename) + "_msk" + ".vrt" image.setWidth(int_ds.RasterXSize) image.setLength(int_ds.RasterYSize) image.setAccessMode("READ") @@ -356,3 +457,264 @@ def mask_interferogram( # Overwrite existing interferogram int_ds = None int_array.astype(np.complex64).tofile(ifgFilename) + + +def brige_components(unwrapped_ifg: str, connected_components: str) -> None: + """ + This routine preforms "bridging' of unwrapped phase connected components + Each component is shifted with its median value + + unwrapped_ifg : str + path to unwrapped interferogram + connected_components : str + path to connected components (labels) + """ + print(f'Do bridging {unwrapped_ifg}') + ifg = gdal.Open(unwrapped_ifg, gdal.GA_Update) + ifg_conn = gdal.Open(connected_components) + + # Load the interferogram and the connected component labels + interferogram = ifg.GetRasterBand(2).ReadAsArray() + labels = ifg_conn.GetRasterBand(1).ReadAsArray() + + # Loop over each connected component in the interferogram + for i in range(1, np.max(labels)+1): + + # Create a binary mask for the current component + mask = np.zeros_like(interferogram, dtype=bool) + mask[labels == i] = True + + # Apply a binary closing operation to the mask + # to bridge any unwrapping errors + mask = morphology.binary_closing(mask) + + # Calculate the median phase value for the current component + median_phase = np.median(interferogram[mask]) + + # Subtract the median phase from the current component + # to remove any phase jumps + interferogram[mask] -= median_phase + + # Write + ifg.GetRasterBand(2).WriteArray(interferogram) + ifg.FlushCache() + ifg = None + + +# MODIFIED FUNCTIONS FROM ISCE2/TOPSPROC, TRIED TO KEEP THEM +# IN AN ORIGINAL FORM + +def unwrap(self, ionParam, use_bridging=True, use_conncomp=True) -> None: + ''' + unwrap lower and upper band interferograms + ref: isce2/components/isceobj/TopsProc/runIon.py#L915 + + M.G. April 2023 : Added option to use connected component 0 + to mask out unreliable pixels in unwrapped interferograms + Added option to use bridging between diconnected unwrapped + phase components + ''' + + # Get number of looks for second multilooking + nrange0 = ionParam.numberRangeLooks0 + nazimuth0 = ionParam.numberAzimuthLooks0 + + # Set number of looks for range0 and azimuth 0 + # same as for range and azimuth to skip multilook unw at first + ionParam.numberRangeLooks0 = ionParam.numberRangeLooks + ionParam.numberAzimuthLooks0 = ionParam.numberAzimuthLooks + + # Run Snaphu unwrapping + runIon.unwrap(self, ionParam) + + # Bridge and mask using connected components + if use_bridging | use_conncomp: + lower_unw_ifg = Path(ionParam.ionDirname, ionParam.lowerDirname, + ionParam.mergedDirname, + self._insar.unwrappedIntFilename) + upper_unw_ifg = Path(ionParam.ionDirname, ionParam.upperDirname, + ionParam.mergedDirname, + self._insar.unwrappedIntFilename) + + for ion_ifg in [lower_unw_ifg, upper_unw_ifg]: + # Shifted unwrapped componentes with their median value + if use_bridging: + brige_components(str(ion_ifg), str(ion_ifg) + '.conncomp') + # Use connected component 0 to mask unwrapped interferograms + if use_conncomp: + runIon.maskUnwrap(str(ion_ifg), str(ion_ifg) + '.conncomp') + + # Multilook + ionParam.numberRangeLooks0 = nrange0 + ionParam.numberAzimuthLooks0 = nazimuth0 + runIon.multilook_unw(self, ionParam, ionParam.mergedDirname) + + +def filt_gaussian(self, ionParam, + coh_threshold=0.5, sigma_rule=2) -> None: + ''' + This function filters image using gaussian filter + + REF: isce2/components/isceobj/TopsProc/runIon.py#L1906 + + M.G. April 2023 : Added option to mask out pixel using coherence + threshold and sigma rule for outlier removal before gaussian + filtering + ''' + + print('filtering ionosphere') + # Using not projected ionosphere + ionfile = Path(ionParam.ionDirname, ionParam.ioncalDirname, + ionParam.ionRawNoProj) + corfile = Path(ionParam.ionDirname, ionParam.ioncalDirname, + ionParam.ionCorNoProj) + + outfile = Path(ionParam.ionDirname, ionParam.ioncalDirname, + ionParam.ionFilt) + + # Read files using gdal + iono_ds = gdal.Open(str(ionfile)) + cor_ds = gdal.Open(str(corfile)) + + ion = iono_ds.GetRasterBand(2).ReadAsArray() + amp = iono_ds.GetRasterBand(1).ReadAsArray() + # Coherence is also stored as 2 band float file + cor = cor_ds.GetRasterBand(2).ReadAsArray() + + length = iono_ds.RasterYSize + width = iono_ds.RasterXSize + + # NOTE from runIon: AFTER COHERENCE IS RESAMPLED + # AT grd2ion, THERE ARE SOME WIRED VALUES + cor[np.nonzero(cor < 0)] = 0.0 + cor[np.nonzero(cor > 1)] = 0.0 + + # Get surface fitted to iono correction, polynomial fitting + if ionParam.ionFit: + ion_fit = runIon.weight_fitting(ion, cor, width, length, + 1, 1, 1, 1, 2, 0.85) + ion -= ion_fit * (ion != 0) + else: + ion_fit = np.zeros(ion.shape) + + # MG. Remove pixels with low coherence + # NOTE output is more smoothed, as + # filtering is less affected with noisy pixels + + if coh_threshold != 0.0: + ion[cor < coh_threshold] = 0 + + # Remove outliers using std + if sigma_rule != 0: + mean = np.mean(np.ma.masked_equal(ion, 0)) + std = np.std(np.ma.masked_equal(ion, 0)) + + ion[ion > mean + sigma_rule*std] = 0 + ion[ion < mean - sigma_rule*std] = 0 + + # Run modifed adaptive filtering MG April 2023 + filt = adaptive_gaussian(ion, cor**14, + ionParam.ionFilteringWinsizeMax, + ionParam.ionFilteringWinsizeMin) + + # Add estimated surface if exists + filt += ion_fit * (filt != 0) + + # Save + ion = np.zeros((length*2, width), dtype=np.float32) + ion[0:length*2:2, :] = amp + ion[1:length*2:2, :] = filt + ion.astype(np.float32).tofile(str(outfile)) + img = isceobj.createImage() + img.load(str(ionfile) + '.xml') + img.filename = str(outfile) + img.extraFilename = str(outfile) + '.vrt' + img.renderHdr() + + # Close gdal files + iono_ds = None + cor_ds = None + + +def adaptive_gaussian(ionos, wgt, size_max, size_min): + ''' + This program performs Gaussian filtering with adaptive window size. + + REF: isce2/components/isceobj/TopsProc/runIon.py#L1846 + + ionos: ionosphere + wgt: weight + size_max: maximum window size + size_min: minimum window size + + MG. NOTE: masking out with zeros is messing with statistics std + and mean, so it is important to mask zeros (aka nans) + before calculation + + ''' + + length = (ionos.shape)[0] + width = (ionos.shape)[1] + flag = (ionos != 0) * (wgt != 0) + ionos *= flag + wgt *= flag + + size_num = 100 + size = np.linspace(size_min, size_max, num=size_num, endpoint=True) + std = np.zeros((length, width, size_num)) + flt = np.zeros((length, width, size_num)) + out = np.zeros((length, width, 1)) + + # calculate filterd image and standard deviation + # sigma of window size: size_max + sigma = size_max / 2.0 + for i in range(size_num): + size2 = np.int32(np.around(size[i])) + if size2 % 2 == 0: + size2 += 1 + if (i+1) % 10 == 0: + print('min win: %4d, max win: %4d, current win: %4d' % ( + np.int32(np.around(size_min)), + np.int32(np.around(size_max)), + size2)) + + g2d = runIon.gaussian(size2, sigma*size2/size_max, scale=1.0) + scale = ss.fftconvolve(wgt, g2d, mode='same') + flt[:, :, i] = ss.fftconvolve( + ionos*wgt, g2d, mode='same') / (scale + (scale == 0)) + + # variance of resulting filtered sample + scale = scale**2 + var = ss.fftconvolve(wgt, g2d**2, mode='same') / (scale + (scale == 0)) + # in case there is a large area without data where scale is very small, + # which leads to wired values in variance + var[~flag] = 0 # MG mask the variance where there are 0s in input data + var[np.nonzero(var < 0)] = 0 + std[:, :, i] = np.sqrt(var) + + std_mv = np.mean(std[np.nonzero(std != 0)], dtype=np.float64) + diff_max = np.amax(np.absolute(std - std_mv)) + std_mv + 1 + std[np.nonzero(std == 0)] = diff_max + + # Find indexes with minimum values along the axis 2 + ixs = np.argmin(np.absolute(std - std_mv), axis=2) + + index = np.nonzero(np.ones((length, width))) + \ + (ixs.reshape(length*width), ) + + out = flt[index] + out = out.reshape((length, width)) + + # MG: use mean value to replace masked value, zeros introduce + # artifacts at the edges in the final output + out[~flag] = np.mean(np.ma.masked_array(out, mask=~flag)) + + # remove artifacts due to varying wgt + size_smt = size_min + if size_smt % 2 == 0: + size_smt += 1 + g2d = runIon.gaussian(size_smt, size_smt/2.0, scale=1.0) + scale = ss.fftconvolve((out != 0), g2d, mode='same') + out2 = ss.fftconvolve(out, g2d, mode='same') / (scale + (scale == 0)) + + return out2 From 40fc06c61799b6cff941309331ac75db24167be4 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 18 Apr 2023 10:14:54 -0700 Subject: [PATCH 02/33] update the import order --- isce2_topsapp/iono_proc.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index fa0848c5..2caf1c87 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -4,21 +4,24 @@ import os import site +import multiprocessing from pathlib import Path from typing import Union -from skimage import morphology + import scipy.signal as ss -import multiprocessing +from skimage import morphology + import numpy as np +from isce.components import isceobj from isce.applications import topsApp from isce.components.isceobj.TopsProc import runIon -from isce.components import isceobj + from isce.components.isceobj.TopsProc.runMergeBursts import ( - interpolateDifferentNumberOfLooks, multilook, - mergeBursts2, mergeBox -) + interpolateDifferentNumberOfLooks, mergeBox, + mergeBursts2, multilook, +) from numpy.typing import NDArray from osgeo import gdal from tqdm import tqdm From 6646ef281c270310ec12e99f788cb6030951afcc Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 18 Apr 2023 10:18:37 -0700 Subject: [PATCH 03/33] import order 2 --- isce2_topsapp/iono_proc.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index 2caf1c87..8fa10bac 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -3,20 +3,19 @@ # California Institute of Technology import os -import site import multiprocessing +import site from pathlib import Path from typing import Union import scipy.signal as ss +import numpy as np from skimage import morphology - -import numpy as np -from isce.components import isceobj from isce.applications import topsApp -from isce.components.isceobj.TopsProc import runIon +from isce.components import isceobj +from isce.components.isceobj.TopsProc import runIon from isce.components.isceobj.TopsProc.runMergeBursts import ( interpolateDifferentNumberOfLooks, mergeBox, mergeBursts2, multilook, From 9f1d8d0a7ee3b152a660965725c3b1de5f464742 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 18 Apr 2023 10:24:11 -0700 Subject: [PATCH 04/33] import order 3 --- isce2_topsapp/iono_proc.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index 8fa10bac..7a2a13e8 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -2,18 +2,18 @@ # Copyright 2023 # California Institute of Technology -import os import multiprocessing +import os import site from pathlib import Path from typing import Union -import scipy.signal as ss import numpy as np -from skimage import morphology +import scipy.signal as ss from isce.applications import topsApp from isce.components import isceobj +from skimage import morphology from isce.components.isceobj.TopsProc import runIon from isce.components.isceobj.TopsProc.runMergeBursts import ( @@ -21,6 +21,7 @@ mergeBursts2, multilook, ) + from numpy.typing import NDArray from osgeo import gdal from tqdm import tqdm From 66fec393c77fe95e2bece2cb5a4bedb674ddd18c Mon Sep 17 00:00:00 2001 From: Marin Govorcin Date: Tue, 18 Apr 2023 10:30:20 -0700 Subject: [PATCH 05/33] Update CHANGELOG.md --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index bf0d7a65..0ca148dc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,9 +11,11 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ### Fixed * For Solid Earth Tide computation, derive coordinates and spacing from geotrans as opposed to latitude/longitude metadata arrays * Include topsapp_iono template. +* Fixed ionosphere computation over water, includes masking conncomp zero, phase bridging, and modified adaptive gaussian filtering ### Added * localize_data within __main__.py added option to use/not use water mask for ionosphere processing +* Added option to estimate burst phase jumps in ionosphere computation ## [0.2.3] From 2003812d3fa7ad523977ef23a5d930b4f0a168aa Mon Sep 17 00:00:00 2001 From: Marin Govorcin Date: Tue, 18 Apr 2023 10:39:05 -0700 Subject: [PATCH 06/33] Update CHANGELOG.md --- CHANGELOG.md | 2 -- 1 file changed, 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0ca148dc..bf0d7a65 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,11 +11,9 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ### Fixed * For Solid Earth Tide computation, derive coordinates and spacing from geotrans as opposed to latitude/longitude metadata arrays * Include topsapp_iono template. -* Fixed ionosphere computation over water, includes masking conncomp zero, phase bridging, and modified adaptive gaussian filtering ### Added * localize_data within __main__.py added option to use/not use water mask for ionosphere processing -* Added option to estimate burst phase jumps in ionosphere computation ## [0.2.3] From cd677829cc3cf973c5085e677292a2ff23726e55 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 18 Apr 2023 10:34:06 -0700 Subject: [PATCH 07/33] fix import order 3 --- isce2_topsapp/iono_proc.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index 7a2a13e8..53c6aaf1 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -8,13 +8,8 @@ from pathlib import Path from typing import Union -import numpy as np -import scipy.signal as ss - from isce.applications import topsApp from isce.components import isceobj -from skimage import morphology - from isce.components.isceobj.TopsProc import runIon from isce.components.isceobj.TopsProc.runMergeBursts import ( interpolateDifferentNumberOfLooks, mergeBox, @@ -22,8 +17,15 @@ ) +import numpy as np from numpy.typing import NDArray + from osgeo import gdal + +import scipy.signal as ss + +from skimage import morphology + from tqdm import tqdm ''' From 8fde8f474ec512032adb917c747bba65155daa25 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 18 Apr 2023 10:37:50 -0700 Subject: [PATCH 08/33] updated changelog --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index bf0d7a65..0ca148dc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,9 +11,11 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ### Fixed * For Solid Earth Tide computation, derive coordinates and spacing from geotrans as opposed to latitude/longitude metadata arrays * Include topsapp_iono template. +* Fixed ionosphere computation over water, includes masking conncomp zero, phase bridging, and modified adaptive gaussian filtering ### Added * localize_data within __main__.py added option to use/not use water mask for ionosphere processing +* Added option to estimate burst phase jumps in ionosphere computation ## [0.2.3] From f95c3da77602240f911b05553c794ff3139fabbc Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 18 Apr 2023 10:39:11 -0700 Subject: [PATCH 09/33] update import order 4 --- isce2_topsapp/iono_proc.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index 53c6aaf1..dddd6fd8 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -10,6 +10,7 @@ from isce.applications import topsApp from isce.components import isceobj +import numpy as np from isce.components.isceobj.TopsProc import runIon from isce.components.isceobj.TopsProc.runMergeBursts import ( interpolateDifferentNumberOfLooks, mergeBox, @@ -17,15 +18,10 @@ ) -import numpy as np from numpy.typing import NDArray - -from osgeo import gdal - import scipy.signal as ss - +from osgeo import gdal from skimage import morphology - from tqdm import tqdm ''' From 37625f51b5bb48f8e644773a0054837e090775ee Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 18 Apr 2023 10:40:35 -0700 Subject: [PATCH 10/33] import order 5 --- isce2_topsapp/iono_proc.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index dddd6fd8..d52503b7 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -9,17 +9,16 @@ from typing import Union from isce.applications import topsApp -from isce.components import isceobj import numpy as np +from isce.components import isceobj from isce.components.isceobj.TopsProc import runIon from isce.components.isceobj.TopsProc.runMergeBursts import ( interpolateDifferentNumberOfLooks, mergeBox, mergeBursts2, multilook, ) - -from numpy.typing import NDArray import scipy.signal as ss +from numpy.typing import NDArray from osgeo import gdal from skimage import morphology from tqdm import tqdm From f0c151a90cee18804ef35b2bf35bb4941ca000d6 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 18 Apr 2023 10:42:36 -0700 Subject: [PATCH 11/33] import order 6 --- isce2_topsapp/iono_proc.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index d52503b7..995b0046 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -8,16 +8,16 @@ from pathlib import Path from typing import Union -from isce.applications import topsApp import numpy as np +from isce.applications import topsApp from isce.components import isceobj from isce.components.isceobj.TopsProc import runIon +import scipy.signal as ss from isce.components.isceobj.TopsProc.runMergeBursts import ( interpolateDifferentNumberOfLooks, mergeBox, mergeBursts2, multilook, ) -import scipy.signal as ss from numpy.typing import NDArray from osgeo import gdal from skimage import morphology From 2097b6574de773f90b2eaab1ce9f04b69a897ad1 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 18 Apr 2023 10:44:36 -0700 Subject: [PATCH 12/33] import order 7 --- isce2_topsapp/iono_proc.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index 995b0046..0bd4dd23 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -11,13 +11,11 @@ import numpy as np from isce.applications import topsApp from isce.components import isceobj -from isce.components.isceobj.TopsProc import runIon import scipy.signal as ss +from isce.components.isceobj.TopsProc import runIon from isce.components.isceobj.TopsProc.runMergeBursts import ( interpolateDifferentNumberOfLooks, mergeBox, - mergeBursts2, multilook, - -) + mergeBursts2, multilook) from numpy.typing import NDArray from osgeo import gdal from skimage import morphology From ab6d3559b6339d2772d3248b74c5f51ea01f8e9c Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 18 Apr 2023 10:45:34 -0700 Subject: [PATCH 13/33] import order 9 --- isce2_topsapp/iono_proc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index 0bd4dd23..c5db8c82 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -9,9 +9,9 @@ from typing import Union import numpy as np +import scipy.signal as ss from isce.applications import topsApp from isce.components import isceobj -import scipy.signal as ss from isce.components.isceobj.TopsProc import runIon from isce.components.isceobj.TopsProc.runMergeBursts import ( interpolateDifferentNumberOfLooks, mergeBox, From a09041ab5fdf5217ed1ff107339daa6d3ad83773 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Wed, 19 Apr 2023 14:39:22 -0700 Subject: [PATCH 14/33] added reference and applied comments --- isce2_topsapp/iono_proc.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index c5db8c82..d716413c 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -6,7 +6,7 @@ import os import site from pathlib import Path -from typing import Union +from typing import Union, Type import numpy as np import scipy.signal as ss @@ -22,6 +22,9 @@ from tqdm import tqdm ''' +Reference: Liang et al. (2019): Ionospheric Correction of InSAR Time Series +Analysis of C-band Sentinel-1 TOPS Data, doi:0.1109/TGRS.2019.2908494 + TODO: add unwrapping changes to ionosphere swath by swath ''' @@ -277,7 +280,7 @@ def merge_bursts( image.renderHdr() -def merge_multilook_bursts(self, +def merge_multilook_bursts(self: Type[topsApp.TopsInSAR], input_dir: str = 'ion/ion_burst', output_filename: str = 'topophase.ion', mergedir: Union[str, Path] = "./merged") -> None: @@ -480,11 +483,10 @@ def brige_components(unwrapped_ifg: str, connected_components: str) -> None: # Create a binary mask for the current component mask = np.zeros_like(interferogram, dtype=bool) - mask[labels == i] = True + mask = (labels == i) # Apply a binary closing operation to the mask - # to bridge any unwrapping errors - mask = morphology.binary_closing(mask) + mask = morphology.binary_closing(mask, structure=None, iterations=1) # Calculate the median phase value for the current component median_phase = np.median(interferogram[mask]) @@ -502,7 +504,8 @@ def brige_components(unwrapped_ifg: str, connected_components: str) -> None: # MODIFIED FUNCTIONS FROM ISCE2/TOPSPROC, TRIED TO KEEP THEM # IN AN ORIGINAL FORM -def unwrap(self, ionParam, use_bridging=True, use_conncomp=True) -> None: +def unwrap(self: Type[topsApp.TopsInSAR], ionParam: Type[runIon.dummy], + use_bridging: bool = True, use_conncomp: bool = True) -> None: ''' unwrap lower and upper band interferograms ref: isce2/components/isceobj/TopsProc/runIon.py#L915 @@ -548,8 +551,8 @@ def unwrap(self, ionParam, use_bridging=True, use_conncomp=True) -> None: runIon.multilook_unw(self, ionParam, ionParam.mergedDirname) -def filt_gaussian(self, ionParam, - coh_threshold=0.5, sigma_rule=2) -> None: +def filt_gaussian(self: Type[topsApp.TopsInSAR], ionParam: Type[runIon.dummy], + coh_threshold: float = 0.5, sigma_rule: int = 2) -> None: ''' This function filters image using gaussian filter @@ -634,7 +637,8 @@ def filt_gaussian(self, ionParam, cor_ds = None -def adaptive_gaussian(ionos, wgt, size_max, size_min): +def adaptive_gaussian(ionos: NDArray, wgt: NDArray, + size_max: int = 200, size_min: int = 100) -> NDArray: ''' This program performs Gaussian filtering with adaptive window size. From b8be8582d439d785886fddc65392d08658e7e260 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Wed, 19 Apr 2023 14:41:51 -0700 Subject: [PATCH 15/33] import order 10 --- isce2_topsapp/iono_proc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index d716413c..931c88e7 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -6,7 +6,7 @@ import os import site from pathlib import Path -from typing import Union, Type +from typing import Type, Union import numpy as np import scipy.signal as ss From 4812032ade56f5f3a1c89b881084ef293ce04bfa Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Wed, 19 Apr 2023 14:44:14 -0700 Subject: [PATCH 16/33] updated CHANGELOG --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3580126a..c353ac13 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,7 +10,7 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ### Fixed * For Solid Earth Tide computation, derive coordinates and spacing from geotrans as opposed to latitude/longitude metadata arrays -* Include topsapp_iono template. +* Include topsapp_iono template. * Increases DEM buffer to .4 from .1 to ensure the extent of at least two bursts (~40 km) are added when retrieving DEM (because estimated footprint can differ from what ISCE2 generates for a GUNW extent) * Fixed ionosphere computation over water, includes masking conncomp zero, phase bridging, and modified adaptive gaussian filtering From 2fb49d1707c085fe53bcf8d74ac766c8a12e97bb Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 24 Apr 2023 09:06:52 -0700 Subject: [PATCH 17/33] add_crossS1A_B_support --- isce2_topsapp/Untitled-1.ipynb | 587 +++++++++++++++++++++++++++++++++ 1 file changed, 587 insertions(+) create mode 100644 isce2_topsapp/Untitled-1.ipynb diff --git a/isce2_topsapp/Untitled-1.ipynb b/isce2_topsapp/Untitled-1.ipynb new file mode 100644 index 00000000..5e6f6f25 --- /dev/null +++ b/isce2_topsapp/Untitled-1.ipynb @@ -0,0 +1,587 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create sensing grid\n", + "\n", + "from mintpy import iono_tec" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import site\n", + "from pathlib import Path\n", + "from isce2_topsapp import iono_proc\n", + "from isce.applications import topsApp\n", + "# Update PATH with ISCE2 applications\n", + "isce_app_path = Path(f\"{site.getsitepackages()[0]}\" \"/isce/applications/\")\n", + "os.environ[\"PATH\"] += \":\" + str(isce_app_path)\n", + "\n", + "tops_app_cmd = f\"{isce_app_path}/topsApp.py\"\n", + "os.chdir('/u/trappist-r0/govorcin/ARIA-tools/IONO/isce/docker_A_B')\n", + "#Load input file\n", + "topsapp = topsApp.TopsInSAR(name=\"topsApp\", cmdline=['topsApp.xml'])\n", + "topsapp.configure()\n", + "\n", + "#Load pickle objekt from 1 step before ion: fine resampale\n", + "topsapp.loadPickleObj('/u/trappist-r0/govorcin/ARIA-tools/IONO/isce/docker_A_B/PICKLE/fineresamp')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ionParam = iono_proc.runIon.setup(topsapp)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "iono_proc.ionSwathBySwath(topsapp, ionParam)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import os\n", + "import datetime as dt\n", + "\n", + "import geopandas as gpd\n", + "import numpy as np\n", + "import pandas as pd\n", + "from matplotlib import pyplot as plt, colorbar, ticker, colors\n", + "from scipy import interpolate, stats\n", + "from shapely.geometry import Point\n", + "\n", + "from mintpy import iono_tec\n", + "from mintpy.objects import sensor, ionex\n", + "from mintpy.simulation import iono\n", + "from mintpy.utils import ptime, readfile, writefile\n", + "plt.rcParams.update({'font.size': 12})\n", + "\n", + "#proj_dir = os.path.expanduser('.')\n", + "#work_dir = os.path.join(proj_dir, 'TEC')\n", + "#os.chdir(work_dir)\n", + "#print('Go to directory:', work_dir)\n", + "\n", + "# aux info\n", + "tec_dir = os.path.expanduser('~/data/aux/IONEX')\n", + "#tf_file = os.path.join(work_dir, 'tframe_left_look.gpkg')\n", + "\n", + "inc_angle = 42\n", + "inc_angle_iono = iono.incidence_angle_ground2iono(inc_angle)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tec_files = []\n", + "for year in range(2014, 2023):\n", + " date_list = ptime.get_date_range(f'{year}0101', f'{year}1231')\n", + " tec_files += iono_tec.download_ionex_files(date_list, tec_dir, sol_code='jpl')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "from numpy.typing import NDArray\n", + "from isce.components.isceobj.Constants import SPEED_OF_LIGHT\n", + "from isce.applications import topsApp\n", + "\n", + "def get_sensingtime_grid(self, lats, lons, hgts) -> NDArray:\n", + "\n", + " #Get frames and orbit\n", + " swathList = self._insar.getValidSwathList(topsapp.swaths)\n", + " frames = []\n", + " for swath in swathList:\n", + " referenceProduct = self._insar.loadProduct( os.path.join(self._insar.fineCoregDirname, 'IW{0}.xml'.format(swath)))\n", + " frames.append(referenceProduct)\n", + "\n", + " # Load orbits\n", + " orb = self._insar.getMergedOrbit(frames)\n", + "\n", + " # Create sensing grid array\n", + " length, width = lat_ds.RasterYSize, lat_ds.RasterXSize \n", + "\n", + " sensing_grid = np.zeros((length, width, 2), dtype=np.float64)\n", + "\n", + " for iy in range(0,length,1):\n", + " for ix in range(0,width,1):\n", + " if lats[iy,ix] !=0:\n", + " pixel_sensing_az_rg = orb.geo2rdr(llh=(lats[iy,ix],\n", + " lons[iy,ix],\n", + " hgts[iy,ix]))\n", + " \n", + " sensing_grid[iy,ix,0] = pixel_sensing_az_rg[0].timestamp()\n", + " sensing_grid[iy,ix,1] = (2*pixel_sensing_az_rg[1]) / SPEED_OF_LIGHT\n", + " \n", + " return sensing_grid" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "import pandas as pd\n", + "import numpy as np\n", + "\n", + "def get_nc_groups(filename_nc : str):\n", + " from osgeo import gdal\n", + " #open nc file\n", + " ds = gdal.Open(filename_nc)\n", + " metadata = ds.GetMetadata()\n", + "\n", + " # Get Groups\n", + " groups_df = pd.DataFrame(columns=['LAYER','PATH'])\n", + " for layer in ds.GetSubDatasets():\n", + " layer_df = pd.DataFrame(data={'LAYER':[layer[0].split(':')[-1].split('/')[-1]],\n", + " 'PATH':[layer[0]]})\n", + " groups_df = pd.concat([groups_df, layer_df], ignore_index=True)\n", + " #close\n", + " ds = None\n", + "\n", + " return groups_df, metadata\n", + "\n", + "def read_nc_group(group_path:str):\n", + " from osgeo import gdal\n", + " # group_path : 'NETCDF:\"$path/file.nc\":group_layer\n", + " # example: 'NETCDF:\"$PATH/S1-GUNW-D-R-021-tops-20230210_20230129-033529-00035E_00034N-PP-9780-v2_0_6.nc\":/science/grids/data/unwrappedPhase'\n", + "\n", + " # open group\n", + " ds = gdal.Open(group_path, gdal.GA_ReadOnly)\n", + "\n", + " # Get group/raster array\n", + " data = ds.ReadAsArray()\n", + " nodata = ds.GetRasterBand(1).GetNoDataValue()\n", + " metadata = ds.GetMetadata()\n", + "\n", + " # Get raster geographical information\n", + " trans = ds.GetGeoTransform()\n", + " xsize = ds.RasterXSize\n", + " ysize = ds.RasterYSize\n", + " snwe = [trans[3] + ysize * trans[5], trans[3],\n", + " trans[0], trans[0] + xsize * trans[1]]\n", + " lon_spacing = trans[1]\n", + " lat_spacing = trans[5]\n", + " \n", + " projection = ds.GetProjection()\n", + "\n", + " # wrap raster info in dict\n", + " raster_dict = {'NODATA' : nodata,\n", + " 'LENGTH' : ysize,\n", + " 'WIDTH' : xsize,\n", + " 'SNWE' : snwe,\n", + " 'LON_SPACING' : lon_spacing,\n", + " 'LAT_SPACING' : lat_spacing,\n", + " 'PROJECTION' : projection}\n", + "\n", + " # close\n", + " ds = None\n", + "\n", + " return data, raster_dict, metadata\n", + "\n", + "# get image extent from attr['snwe'] for plotting\n", + "def snwe_to_extent(snwe):\n", + " extent = [snwe[2], snwe[3], snwe[0], snwe[1]]\n", + " \n", + " return extent" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "os.chdir('/u/trappist-r0/govorcin/ARIA-tools/IONO/isce/docker_test2')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "os.listdir('S1-GUNW-A-R-064-tops-20211120_20211108-014951-00119W_00030N-PP-aed8-v2_0_6')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "get_nc_groups('S1-GUNW-A-R-064-tops-20211120_20211108-014951-00119W_00030N-PP-aed8-v2_0_6/S1-GUNW-A-R-064-tops-20211120_20211108-014951-00119W_00030N-PP-aed8-v2_0_6.nc')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from osgeo import gdal\n", + "ds = gdal.Open('S1-GUNW-A-R-064-tops-20211120_20211108-014951-00119W_00030N-PP-aed8-v2_0_6/S1-GUNW-A-R-064-tops-20211120_20211108-014951-00119W_00030N-PP-aed8-v2_0_6.nc')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for layer in ds.GetSubDatasets():\n", + " print(layer)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "metadata = ds.GetMetadata()\n", + "metadata" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import h5py" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "file = h5py.File('merged/metadata.h5', 'r')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dateime(midtight) + timedelta(second=10000)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data = file['cube']\n", + "print(data.keys())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data['lats'], data['lons'], data['heights']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data['lats'].shape[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data['lats'][:]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sensing_grid = np.zeros((data['lats'].shape[0],\n", + " data['lons'].shape[0], \n", + " data['heights'].shape[0]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "plt.imshow(sensing_grid[:,:,0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nz_idx = np.nonzero(data['lats'][:])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "os.path.join(os.getcwd(),'PICKLE/geocode')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import h5py\n", + "import os\n", + "from isce.components.isceobj.Constants import SPEED_OF_LIGHT\n", + "from datetime import timedelta\n", + "from isce.applications import topsApp\n", + "\n", + "file = h5py.File('merged/metadata.h5', 'r')\n", + "data = file['cube']\n", + "\n", + "\n", + "#Load input file\n", + "topsapp = topsApp.TopsInSAR(name=\"topsApp\", cmdline=['topsApp.xml'])\n", + "topsapp.configure()\n", + "\n", + "#Load pickle object\n", + "topsapp.loadPickleObj(os.path.join(os.getcwd(),'PICKLE/geocode'))\n", + "\n", + "#Get frames and orbit\n", + "swathList = topsapp._insar.getValidSwathList(topsapp.swaths)\n", + "frames = []\n", + "for swath in swathList:\n", + " referenceProduct = topsapp._insar.loadProduct(os.path.join(topsapp._insar.fineCoregDirname, 'IW{0}.xml'.format(swath)))\n", + " frames.append(referenceProduct)\n", + "\n", + "# Load orbits\n", + "orb = topsapp._insar.getMergedOrbit(frames)\n", + "\n", + "sensing_grid = np.zeros((data['lats'].shape[0],\n", + " data['lons'].shape[0], \n", + " data['heights'].shape[0]))\n", + "\n", + "for iy, lat in enumerate(data['lats'][:]):\n", + " for ix, lon in enumerate(data['lons'][:]):\n", + " for iz, hgt in enumerate(data['heights'][:]):\n", + " pixel_sensing_az_rg = orb.geo2rdr(llh=(lat,lon,hgt)) # approx, should alwys start from 1 swath\n", + " sensing = pixel_sensing_az_rg[0] + timedelta(seconds=(pixel_sensing_az_rg[1]) / SPEED_OF_LIGHT) # azimuth_sensing_time +( ix * PRF + range/speed of light)\n", + " sensing_grid[iy,ix,iz] = sensing.timestamp()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(lat,lon,hgt) 1/PRF" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "orb.geo2rdr(llh=(30,-115,1050))[0] - orb.geo2rdr(llh=(30,-115,1000))[0] " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pixel_sensing_az_rg[1] / SPEED_OF_LIGHT" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "frames[1].bursts.burst1.firstValidSample * (1/frames[1].bursts.burst1.prf)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "frames[1].bursts.burst1.lastValidSample * (1/frames[1].bursts.burst1.prf) + pixel_sensing_az_rg[1] / SPEED_OF_LIGHT" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "np.broadcast_to(data['lats'][:], (72,47))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from datetime import datetime" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "datetime.fromtimestamp(sensing_grid[:,:,3])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from datetime import datetime\n", + "[datetime.fromtimestamp(i) for i in sensing_grid[0,:,3]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.imshow(sensing_grid[:,:,])\n", + "plt.colorbar()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "datetime.fromtimestamp()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "lats = data['lats'][:]\n", + "lons = data['lons'][:]\n", + "heights = data['heights'][:]\n", + "\n", + "# Reshape the arrays to create a 3D grid\n", + "ny, nx, nz = (lats.shape[0], \n", + " lons.shape[0], \n", + " heights.shape[0])\n", + "lats = lats.reshape(ny, nx, 1)\n", + "lons = lons.reshape(ny, nx, 1)\n", + "heights = heights.reshape(1, 1, nz)\n", + "\n", + "\n", + "pixels_sensing_az_rg = orb.geo2rdr(llh=(lats, lons, heights))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "coords = np.stack((data['lats'][:], data['lons'][:], data['heights'][:]), axis=-1)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "topsapp_env", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.16" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} From a4dab420f6e4bd27ac31487b171f5e23f14f5f6a Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 24 Apr 2023 09:08:14 -0700 Subject: [PATCH 18/33] add_cross_S1_B_support --- isce2_topsapp/iono_proc.py | 410 ++++++++++++++++++++++++++++++++++--- 1 file changed, 381 insertions(+), 29 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index 931c88e7..9ea0b7ff 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -8,6 +8,8 @@ from pathlib import Path from typing import Type, Union +import shutil +import datetime import numpy as np import scipy.signal as ss from isce.applications import topsApp @@ -15,6 +17,7 @@ from isce.components.isceobj.TopsProc import runIon from isce.components.isceobj.TopsProc.runMergeBursts import ( interpolateDifferentNumberOfLooks, mergeBox, + adjustValidWithLooks, mergeBurstsVirtual, mergeBursts2, multilook) from numpy.typing import NDArray from osgeo import gdal @@ -25,13 +28,10 @@ Reference: Liang et al. (2019): Ionospheric Correction of InSAR Time Series Analysis of C-band Sentinel-1 TOPS Data, doi:0.1109/TGRS.2019.2908494 -TODO: add unwrapping changes to ionosphere swath by swath +List of parameters for ionospheric correction: + isce2/examples/input_files/topsApp.xml#L182 ''' -# List of parameters for ionospheric correction: -# isce2/examples/input_files/topsApp.xml#L182 - - GEOCODE_LIST_ION = ["merged/topophase.ion"] omp_env = os.getenv('OMP_NUM_THREADS') @@ -120,10 +120,9 @@ def iono_processing( else: # This mode is used for cross # Sentinel-1A/B interferogram - # NOTE: need to re-write this step to include - # changes with unwrapping, aka bridging - # and using conncomp=0 to mask noise - runIon.ionSwathBySwath(topsapp, ionParam) + # runIon.ionSwathBySwath(topsapp, ionParam) + ionSwathBySwath(topsapp, ionParam, + use_bridging=True, conncomp_flag=conncomp_flag) # Run iono step grd2ion # Resample ionosphere from ground to iono layer @@ -196,7 +195,7 @@ def get_burst(x): output_dir = tops_dir / f"mask/{swath}" output_dir.mkdir(parents=True, exist_ok=True) output_file = output_dir / f"msk_{burst}.rdr" - # Projct mask to radar coordinate space + # Project mask to radar coordinate space raster_geo2radar(mask_filename, str(lat), lon, str(output_file)) pbar.update() @@ -243,7 +242,6 @@ def merge_bursts( interpolateDifferentNumberOfLooks for raw ionosphere merge. Cubic creates artifacts around the edges. - """ mergedIfgname = "topophase.flat" ######################################### @@ -287,7 +285,6 @@ def merge_multilook_bursts(self: Type[topsApp.TopsInSAR], """ Reference: isce2/components/isceobj/TopsProc/runMergeBursts.py#L776 - """ if (self.numberRangeLooks == 1) and (self.numberAzimuthLooks == 1): @@ -327,8 +324,6 @@ def merge_multilook_bursts(self: Type[topsApp.TopsInSAR], # UTILITIES FOR MASKING - - def raster_geo2radar( rasterFilename: Union[str, Path], latFilename: Union[str, Path], @@ -368,14 +363,12 @@ def raster_geo2radar( lats = lat_ds.ReadAsArray() # Translate raster from geographical to radar coordinate space - lineIdx = np.int32( + lineIdx = np.array( (lats - mask_ds.GetGeoTransform()[3]) / - mask_ds.GetGeoTransform()[5] + 0.5 - ) - sampleIdx = np.int32( + mask_ds.GetGeoTransform()[5] + 0.5, dtype=np.int32) + sampleIdx = np.array( (lons - mask_ds.GetGeoTransform()[0]) / - mask_ds.GetGeoTransform()[1] + 0.5 - ) + mask_ds.GetGeoTransform()[1] + 0.5, dtype=np.int32) inboundIndex = np.logical_and( np.logical_and(lineIdx >= 0, lineIdx <= mask_ds.RasterYSize - 1), np.logical_and(sampleIdx >= 0, sampleIdx <= mask_ds.RasterXSize - 1), @@ -460,10 +453,12 @@ def mask_interferogram( int_array.astype(np.complex64).tofile(ifgFilename) -def brige_components(unwrapped_ifg: str, connected_components: str) -> None: +def brige_components(unwrapped_ifg: str, + connected_components: str) -> None: """ - This routine preforms "bridging' of unwrapped phase connected components - Each component is shifted with its median value + This routine preforms "bridging' of unwrapped phase + connected components. Each component is shifted with + its median value unwrapped_ifg : str path to unwrapped interferogram @@ -475,7 +470,7 @@ def brige_components(unwrapped_ifg: str, connected_components: str) -> None: ifg_conn = gdal.Open(connected_components) # Load the interferogram and the connected component labels - interferogram = ifg.GetRasterBand(2).ReadAsArray() + interferogram = ifg.GetRasterBand(2).ReadAsArray() # unw Phase labels = ifg_conn.GetRasterBand(1).ReadAsArray() # Loop over each connected component in the interferogram @@ -486,7 +481,7 @@ def brige_components(unwrapped_ifg: str, connected_components: str) -> None: mask = (labels == i) # Apply a binary closing operation to the mask - mask = morphology.binary_closing(mask, structure=None, iterations=1) + mask = morphology.binary_closing(mask) # Calculate the median phase value for the current component median_phase = np.median(interferogram[mask]) @@ -504,8 +499,10 @@ def brige_components(unwrapped_ifg: str, connected_components: str) -> None: # MODIFIED FUNCTIONS FROM ISCE2/TOPSPROC, TRIED TO KEEP THEM # IN AN ORIGINAL FORM -def unwrap(self: Type[topsApp.TopsInSAR], ionParam: Type[runIon.dummy], - use_bridging: bool = True, use_conncomp: bool = True) -> None: +def unwrap(self: Type[topsApp.TopsInSAR], + ionParam: Type[runIon.dummy], + use_bridging: bool = True, + use_conncomp: bool = True) -> None: ''' unwrap lower and upper band interferograms ref: isce2/components/isceobj/TopsProc/runIon.py#L915 @@ -551,8 +548,10 @@ def unwrap(self: Type[topsApp.TopsInSAR], ionParam: Type[runIon.dummy], runIon.multilook_unw(self, ionParam, ionParam.mergedDirname) -def filt_gaussian(self: Type[topsApp.TopsInSAR], ionParam: Type[runIon.dummy], - coh_threshold: float = 0.5, sigma_rule: int = 2) -> None: +def filt_gaussian(self: Type[topsApp.TopsInSAR], + ionParam: Type[runIon.dummy], + coh_threshold: float = 0.5, + sigma_rule: int = 2) -> None: ''' This function filters image using gaussian filter @@ -591,6 +590,11 @@ def filt_gaussian(self: Type[topsApp.TopsInSAR], ionParam: Type[runIon.dummy], cor[np.nonzero(cor > 1)] = 0.0 # Get surface fitted to iono correction, polynomial fitting + # weight_fitting(iono_phase, coherence, width, length, + # num_range_looks-input, num_azimuth_looks-input + # num_range_looks-output, num_azimuth_looks-output + # polynomial_order, coherence_threshold) + if ionParam.ionFit: ion_fit = runIon.weight_fitting(ion, cor, width, length, 1, 1, 1, 1, 2, 0.85) @@ -720,3 +724,351 @@ def adaptive_gaussian(ionos: NDArray, wgt: NDArray, out2 = ss.fftconvolve(out, g2d, mode='same') / (scale + (scale == 0)) return out2 + + +def ionSwathBySwath(self: Type[topsApp.TopsInSAR], + ionParam: Type[runIon.dummy], + use_bridging: bool = True, + conncomp_flag: bool = True) -> None: + ''' + This routine merge, unwrap and compute ionosphere swath by swath, + and then adjust phase difference between adjacent swaths caused + by relative range timing error between adjacent swaths. + + This routine includes the following steps in the merged-swath + processing: + + merge(self, ionParam) + unwrap(self, ionParam) + ionosphere(self, ionParam) + + MG: did little clean up of the code, replaced step 2 & 3 with + functions + ''' + + ######################################### + # SET PARAMETERS HERE + numberRangeLooks = ionParam.numberRangeLooks + numberAzimuthLooks = ionParam.numberAzimuthLooks + numberRangeLooks0 = ionParam.numberRangeLooks0 + numberAzimuthLooks0 = ionParam.numberAzimuthLooks0 + mergedDirname = ionParam.mergedDirname # MG + ioncalDirname = ionParam.ioncalDirname # MG + + # THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self) + corThresholdSwathAdj = 0.85 + + ######################################### + + print('computing ionosphere swath by swath') + # if ionParam.calIonWithMerged == False: + warningInfo = (f'{datetime.datetime.now()} calculating ionosphere' + f'swath by swath, there may be slight phase error' + f'between subswaths\n') + with open(os.path.join(ionParam.ionDirname, ionParam.warning), 'a') as f: + f.write(warningInfo) + + # get bursts + numValidSwaths = 0 + swathList = self._insar.getValidSwathList(self.swaths) + for swath in swathList: + minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1) + if minBurst == maxBurst: + # print('Skipping processing of swath {0}'.format(swath)) + continue + numValidSwaths += 1 + + if numValidSwaths <= 1: + raise Exception( + 'There are less than one subswaths, no need to use swath-by-swath' + 'method to compute ionosphere!') + else: + xmlDirname = os.path.join( + ionParam.ionDirname, ionParam.lowerDirname, ionParam.fineIfgDirname) + + merge_kwargs = {'numberRangeLooks': ionParam.numberRangeLooks, + 'numberAzimuthLooks': ionParam.numberAzimuthLooks} + + (box, burstValidBox, + burstValidBox2, frames) = runIon.getMergeBox(self, xmlDirname, + **merge_kwargs) + + # compute ionosphere swath by swath + corList = [] + ampList = [] + ionosList = [] + nswath = len(swathList) + ii = -1 + for i in range(nswath): + swath = swathList[i] + # MG Change merged dir for each swath + ionParam.mergedDirname = mergedDirname + '_IW{0}'.format(swath) + ionParam.ioncalDirname = ioncalDirname + '_IW{0}'.format(swath) + + minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1) + if minBurst == maxBurst: + print('Skipping processing of swath {0}'.format(swath)) + continue + else: + ii += 1 + + ######################################################## + # STEP 1. MERGE THE BURSTS OF A SWATH + ######################################################## + + dirs = [ionParam.lowerDirname, ionParam.upperDirname] + for dirx in dirs: + outputFilename = self._insar.mergedIfgname + outputDirname = os.path.join( + ionParam.ionDirname, dirx, ionParam.mergedDirname) + os.makedirs(outputDirname, exist_ok=True) + suffix = '.full' + if (numberRangeLooks0 == 1) and (numberAzimuthLooks0 == 1): + suffix = '' + + # merge + burstPattern = 'burst_%02d.int' + burstDirname = os.path.join( + ionParam.ionDirname, dirx, ionParam.fineIfgDirname) + + ifg = self._insar.loadProduct(os.path.join( + burstDirname, 'IW{0}.xml'.format(swath))) + bst = [os.path.join(burstDirname, f'IW{swath}', + burstPattern % (x+1)) for x in range(minBurst, + maxBurst)] + + # doing adjustment before use + adjustment_rvalid = np.int32(np.around(numberRangeLooks/8.0)) + adjustValidWithLooks([ifg], box, + numberAzimuthLooks, numberRangeLooks, + edge=0, avalid='strict', + rvalid=adjustment_rvalid) + mergeBurstsVirtual([ifg], [bst], box, os.path.join( + outputDirname, outputFilename+suffix)) + + # take looks + if suffix not in ['', None]: + multilook(os.path.join(outputDirname, outputFilename+suffix), + os.path.join(outputDirname, outputFilename), + numberAzimuthLooks0, + numberRangeLooks0) + else: + print('skipping multilooking') + + # The orginal coherence calculated by topsApp.py is not good at all, + # use the following coherence instead + lowerintfile = os.path.join(ionParam.ionDirname, + ionParam.lowerDirname, + ionParam.mergedDirname, + self._insar.mergedIfgname) + upperintfile = os.path.join(ionParam.ionDirname, + ionParam.upperDirname, + ionParam.mergedDirname, + self._insar.mergedIfgname) + corfile = os.path.join(ionParam.ionDirname, + ionParam.lowerDirname, + ionParam.mergedDirname, + self._insar.correlationFilename) + + img = isceobj.createImage() + img.load(lowerintfile + '.xml') + width = img.width + length = img.length + lowerint = np.fromfile( + lowerintfile, dtype=np.complex64).reshape(length, width) + upperint = np.fromfile( + upperintfile, dtype=np.complex64).reshape(length, width) + + ######################################################## + # slight filtering to improve the estimation accuarcy + # of swath difference + if 1 and shutil.which('psfilt1') != None: + cmd1 = 'mv {} tmp'.format(lowerintfile) + cmd2 = 'psfilt1 tmp {} {} .3 32 8'.format(lowerintfile, width) + cmd3 = 'rm tmp' + cmd4 = 'mv {} tmp'.format(upperintfile) + cmd5 = 'psfilt1 tmp {} {} .3 32 8'.format(upperintfile, width) + cmd6 = 'rm tmp' + + runIon.runCmd(cmd1) + runIon.runCmd(cmd2) + runIon.runCmd(cmd3) + runIon.runCmd(cmd4) + runIon.runCmd(cmd5) + runIon.runCmd(cmd6) + ############################################################### + + # compute coherence only using interferogram + # here I use differential interferogram of lower and upper band + # interferograms so that coherence is not affected by fringes + cord = runIon.cal_coherence( + lowerint*np.conjugate(upperint), win=3, edge=4) + cor = np.zeros((length*2, width), dtype=np.float32) + cor[0:length*2:2, + :] = np.sqrt((np.absolute(lowerint)+np.absolute(upperint))/2.0) + cor[1:length*2:2, :] = cord + cor.astype(np.float32).tofile(corfile) + + # img = isceobj.Image.createUnwImage() + img = isceobj.createOffsetImage() + img.setFilename(corfile) + img.extraFilename = corfile + '.vrt' + img.setWidth(width) + img.setLength(length) + img.renderHdr() + + # MG: replace unwrap, and compute iono code as + # the only change is the iono merged directory + + # STEP 2. UNWRAP SWATH INTERFEROGRAM + unwrap(self, ionParam, + use_bridging=use_bridging, + use_conncomp=conncomp_flag) + + # STEP 3. COMPUTE IONOSPHERE + runIon.ionosphere(self, ionParam) + + # Load the results of ionosphere computation + outDir = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname) + outFilename = os.path.join(outDir, ionParam.ionRawNoProj) + corFilename = os.path.join(outDir, ionParam.ionCorNoProj) + + img = isceobj.createImage() + img.load(outFilename + '.xml') + width = img.width + length = img.length + + ionos = (np.fromfile(outFilename, dtype=np.float32).reshape( + length*2, width))[1:length*2:2, :] + amp = (np.fromfile(outFilename, dtype=np.float32).reshape( + length*2, width))[0:length*2:2, :] + cor = (np.fromfile(corFilename, dtype=np.float32).reshape( + length*2, width))[1:length*2:2, :] + + ionosList.append(ionos) + corList.append(cor) + ampList.append(amp) + + # MG Get the last lower unwrapped + # interferogram path + lowerUnwfile = os.path.join(ionParam.ionDirname, + ionParam.lowerDirname, + ionParam.mergedDirname, + self._insar.unwrappedIntFilename) + + # MG: use image size from lower unwrapped + # interferogram swath + img = isceobj.createImage() + img.load(lowerUnwfile + '.xml') + width = img.width + length = img.length + + # MG: reset ionParm mergedDirname + ionParam.mergedDirname = mergedDirname + ionParam.ioncalDirname = ioncalDirname + + # do adjustment between ajacent swaths + if numValidSwaths == 3: + adjustList = [ionosList[0], ionosList[2]] + else: + adjustList = [ionosList[0]] + for adjdata in adjustList: + masked_cor = corList[1] > corThresholdSwathAdj + index = np.nonzero( + (adjdata != 0) * (ionosList[1] != 0) * masked_cor) + if index[0].size < 5: + print(f'WARNING: too few samples available for adjustment' + f'between swaths: {index[0].size} with coherence ' + f'threshold: {corThresholdSwathAdj}') + print(' no adjustment made!') + print( + ' to do ajustment, please consider using' + ' lower coherence threshold') + else: + print(f'number of samples available for adjustment in the' + f' overlap area: {index[0].size}') + + # use weighted mean instead + wgt = corList[1][index]**14 + diff = np.sum((ionosList[1] - adjdata)[index] * wgt / + np.sum(wgt, dtype=np.float64), dtype=np.float64) + + index2 = np.nonzero(adjdata != 0) + adjdata[index2] = adjdata[index2] + diff + + # get merged ionosphere + ampMerged = np.zeros((length, width), dtype=np.float32) + corMerged = np.zeros((length, width), dtype=np.float32) + ionosMerged = np.zeros((length, width), dtype=np.float32) + for i in range(numValidSwaths): + nBurst = len(burstValidBox[i]) + for j in range(nBurst): + + # index after multi-looking in merged image, index + # starts from 1 + first_line = np.int32( + np.around((burstValidBox[i][j][0] - 1) / + numberAzimuthLooks + 1)) + last_line = np.int32( + np.around(burstValidBox[i][j][1] / + numberAzimuthLooks)) + first_sample = np.int32( + np.around((burstValidBox[i][j][2] - 1) / + numberRangeLooks + 1)) + last_sample = np.int32( + np.around(burstValidBox[i][j][3] / + numberRangeLooks)) + + corMerged[first_line-1:last_line-1+1, + first_sample-1:last_sample-1+1] = \ + corList[i][first_line-1:last_line-1 + + 1, first_sample-1:last_sample-1+1] + + ampMerged[first_line-1:last_line-1+1, + first_sample-1:last_sample-1+1] = \ + ampList[i][first_line-1:last_line-1 + + 1, first_sample-1:last_sample-1+1] + + ionosMerged[first_line-1:last_line-1+1, + first_sample-1:last_sample-1+1] = \ + ionosList[i][first_line-1:last_line - + 1+1, first_sample-1:last_sample-1+1] + + # remove an empirical ramp + # MG this part can be commented out as ramp is not applied to iono + # leave it here for now + if ionParam.rampRemovel != 0: + warningInfo = (f'{datetime.datetime.now()} calculating ionosphere for ' + f'cross S-1A/B interferogram, an empirical ramp is ' + f'removed from estimated ionosphere\n') + + warning_file = os.path.join(ionParam.ionDirname, ionParam.warning) + with open(warning_file, 'a') as f: + f.write(warningInfo) + + abramp = runIon.cal_cross_ab_ramp( + swathList, box[1], numberRangeLooks, ionParam.passDirection) + if ionParam.rampRemovel == -1: + abramp *= -1.0 + # currently do not apply this + # ionosMerged -= abramp[None, :] + + # dump ionosphere + outDir = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname) + os.makedirs(outDir, exist_ok=True) + outFilename = os.path.join(outDir, ionParam.ionRawNoProj) + ion = np.zeros((length*2, width), dtype=np.float32) + ion[0:length*2:2, :] = ampMerged + ion[1:length*2:2, :] = ionosMerged + ion.astype(np.float32).tofile(outFilename) + img.filename = outFilename + img.extraFilename = outFilename + '.vrt' + img.renderHdr() + + # dump coherence + outFilename = os.path.join(outDir, ionParam.ionCorNoProj) + ion[1:length*2:2, :] = corMerged + ion.astype(np.float32).tofile(outFilename) + img.filename = outFilename + img.extraFilename = outFilename + '.vrt' + img.renderHdr() From cd9b5c49214b3dd77969f6a66d98df0d3e98e38b Mon Sep 17 00:00:00 2001 From: Marin Govorcin Date: Mon, 24 Apr 2023 09:12:58 -0700 Subject: [PATCH 19/33] Delete Untitled-1.ipynb --- isce2_topsapp/Untitled-1.ipynb | 587 --------------------------------- 1 file changed, 587 deletions(-) delete mode 100644 isce2_topsapp/Untitled-1.ipynb diff --git a/isce2_topsapp/Untitled-1.ipynb b/isce2_topsapp/Untitled-1.ipynb deleted file mode 100644 index 5e6f6f25..00000000 --- a/isce2_topsapp/Untitled-1.ipynb +++ /dev/null @@ -1,587 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Create sensing grid\n", - "\n", - "from mintpy import iono_tec" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "import site\n", - "from pathlib import Path\n", - "from isce2_topsapp import iono_proc\n", - "from isce.applications import topsApp\n", - "# Update PATH with ISCE2 applications\n", - "isce_app_path = Path(f\"{site.getsitepackages()[0]}\" \"/isce/applications/\")\n", - "os.environ[\"PATH\"] += \":\" + str(isce_app_path)\n", - "\n", - "tops_app_cmd = f\"{isce_app_path}/topsApp.py\"\n", - "os.chdir('/u/trappist-r0/govorcin/ARIA-tools/IONO/isce/docker_A_B')\n", - "#Load input file\n", - "topsapp = topsApp.TopsInSAR(name=\"topsApp\", cmdline=['topsApp.xml'])\n", - "topsapp.configure()\n", - "\n", - "#Load pickle objekt from 1 step before ion: fine resampale\n", - "topsapp.loadPickleObj('/u/trappist-r0/govorcin/ARIA-tools/IONO/isce/docker_A_B/PICKLE/fineresamp')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ionParam = iono_proc.runIon.setup(topsapp)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "iono_proc.ionSwathBySwath(topsapp, ionParam)\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib inline\n", - "import os\n", - "import datetime as dt\n", - "\n", - "import geopandas as gpd\n", - "import numpy as np\n", - "import pandas as pd\n", - "from matplotlib import pyplot as plt, colorbar, ticker, colors\n", - "from scipy import interpolate, stats\n", - "from shapely.geometry import Point\n", - "\n", - "from mintpy import iono_tec\n", - "from mintpy.objects import sensor, ionex\n", - "from mintpy.simulation import iono\n", - "from mintpy.utils import ptime, readfile, writefile\n", - "plt.rcParams.update({'font.size': 12})\n", - "\n", - "#proj_dir = os.path.expanduser('.')\n", - "#work_dir = os.path.join(proj_dir, 'TEC')\n", - "#os.chdir(work_dir)\n", - "#print('Go to directory:', work_dir)\n", - "\n", - "# aux info\n", - "tec_dir = os.path.expanduser('~/data/aux/IONEX')\n", - "#tf_file = os.path.join(work_dir, 'tframe_left_look.gpkg')\n", - "\n", - "inc_angle = 42\n", - "inc_angle_iono = iono.incidence_angle_ground2iono(inc_angle)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tec_files = []\n", - "for year in range(2014, 2023):\n", - " date_list = ptime.get_date_range(f'{year}0101', f'{year}1231')\n", - " tec_files += iono_tec.download_ionex_files(date_list, tec_dir, sol_code='jpl')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "\n", - "from numpy.typing import NDArray\n", - "from isce.components.isceobj.Constants import SPEED_OF_LIGHT\n", - "from isce.applications import topsApp\n", - "\n", - "def get_sensingtime_grid(self, lats, lons, hgts) -> NDArray:\n", - "\n", - " #Get frames and orbit\n", - " swathList = self._insar.getValidSwathList(topsapp.swaths)\n", - " frames = []\n", - " for swath in swathList:\n", - " referenceProduct = self._insar.loadProduct( os.path.join(self._insar.fineCoregDirname, 'IW{0}.xml'.format(swath)))\n", - " frames.append(referenceProduct)\n", - "\n", - " # Load orbits\n", - " orb = self._insar.getMergedOrbit(frames)\n", - "\n", - " # Create sensing grid array\n", - " length, width = lat_ds.RasterYSize, lat_ds.RasterXSize \n", - "\n", - " sensing_grid = np.zeros((length, width, 2), dtype=np.float64)\n", - "\n", - " for iy in range(0,length,1):\n", - " for ix in range(0,width,1):\n", - " if lats[iy,ix] !=0:\n", - " pixel_sensing_az_rg = orb.geo2rdr(llh=(lats[iy,ix],\n", - " lons[iy,ix],\n", - " hgts[iy,ix]))\n", - " \n", - " sensing_grid[iy,ix,0] = pixel_sensing_az_rg[0].timestamp()\n", - " sensing_grid[iy,ix,1] = (2*pixel_sensing_az_rg[1]) / SPEED_OF_LIGHT\n", - " \n", - " return sensing_grid" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from pathlib import Path\n", - "import pandas as pd\n", - "import numpy as np\n", - "\n", - "def get_nc_groups(filename_nc : str):\n", - " from osgeo import gdal\n", - " #open nc file\n", - " ds = gdal.Open(filename_nc)\n", - " metadata = ds.GetMetadata()\n", - "\n", - " # Get Groups\n", - " groups_df = pd.DataFrame(columns=['LAYER','PATH'])\n", - " for layer in ds.GetSubDatasets():\n", - " layer_df = pd.DataFrame(data={'LAYER':[layer[0].split(':')[-1].split('/')[-1]],\n", - " 'PATH':[layer[0]]})\n", - " groups_df = pd.concat([groups_df, layer_df], ignore_index=True)\n", - " #close\n", - " ds = None\n", - "\n", - " return groups_df, metadata\n", - "\n", - "def read_nc_group(group_path:str):\n", - " from osgeo import gdal\n", - " # group_path : 'NETCDF:\"$path/file.nc\":group_layer\n", - " # example: 'NETCDF:\"$PATH/S1-GUNW-D-R-021-tops-20230210_20230129-033529-00035E_00034N-PP-9780-v2_0_6.nc\":/science/grids/data/unwrappedPhase'\n", - "\n", - " # open group\n", - " ds = gdal.Open(group_path, gdal.GA_ReadOnly)\n", - "\n", - " # Get group/raster array\n", - " data = ds.ReadAsArray()\n", - " nodata = ds.GetRasterBand(1).GetNoDataValue()\n", - " metadata = ds.GetMetadata()\n", - "\n", - " # Get raster geographical information\n", - " trans = ds.GetGeoTransform()\n", - " xsize = ds.RasterXSize\n", - " ysize = ds.RasterYSize\n", - " snwe = [trans[3] + ysize * trans[5], trans[3],\n", - " trans[0], trans[0] + xsize * trans[1]]\n", - " lon_spacing = trans[1]\n", - " lat_spacing = trans[5]\n", - " \n", - " projection = ds.GetProjection()\n", - "\n", - " # wrap raster info in dict\n", - " raster_dict = {'NODATA' : nodata,\n", - " 'LENGTH' : ysize,\n", - " 'WIDTH' : xsize,\n", - " 'SNWE' : snwe,\n", - " 'LON_SPACING' : lon_spacing,\n", - " 'LAT_SPACING' : lat_spacing,\n", - " 'PROJECTION' : projection}\n", - "\n", - " # close\n", - " ds = None\n", - "\n", - " return data, raster_dict, metadata\n", - "\n", - "# get image extent from attr['snwe'] for plotting\n", - "def snwe_to_extent(snwe):\n", - " extent = [snwe[2], snwe[3], snwe[0], snwe[1]]\n", - " \n", - " return extent" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "os.chdir('/u/trappist-r0/govorcin/ARIA-tools/IONO/isce/docker_test2')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "os.listdir('S1-GUNW-A-R-064-tops-20211120_20211108-014951-00119W_00030N-PP-aed8-v2_0_6')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "get_nc_groups('S1-GUNW-A-R-064-tops-20211120_20211108-014951-00119W_00030N-PP-aed8-v2_0_6/S1-GUNW-A-R-064-tops-20211120_20211108-014951-00119W_00030N-PP-aed8-v2_0_6.nc')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from osgeo import gdal\n", - "ds = gdal.Open('S1-GUNW-A-R-064-tops-20211120_20211108-014951-00119W_00030N-PP-aed8-v2_0_6/S1-GUNW-A-R-064-tops-20211120_20211108-014951-00119W_00030N-PP-aed8-v2_0_6.nc')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for layer in ds.GetSubDatasets():\n", - " print(layer)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "metadata = ds.GetMetadata()\n", - "metadata" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import h5py" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "file = h5py.File('merged/metadata.h5', 'r')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dateime(midtight) + timedelta(second=10000)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "data = file['cube']\n", - "print(data.keys())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "data['lats'], data['lons'], data['heights']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "data['lats'].shape[0]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "data['lats'][:]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sensing_grid = np.zeros((data['lats'].shape[0],\n", - " data['lons'].shape[0], \n", - " data['heights'].shape[0]))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from matplotlib import pyplot as plt\n", - "plt.imshow(sensing_grid[:,:,0])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "nz_idx = np.nonzero(data['lats'][:])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "os.path.join(os.getcwd(),'PICKLE/geocode')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import h5py\n", - "import os\n", - "from isce.components.isceobj.Constants import SPEED_OF_LIGHT\n", - "from datetime import timedelta\n", - "from isce.applications import topsApp\n", - "\n", - "file = h5py.File('merged/metadata.h5', 'r')\n", - "data = file['cube']\n", - "\n", - "\n", - "#Load input file\n", - "topsapp = topsApp.TopsInSAR(name=\"topsApp\", cmdline=['topsApp.xml'])\n", - "topsapp.configure()\n", - "\n", - "#Load pickle object\n", - "topsapp.loadPickleObj(os.path.join(os.getcwd(),'PICKLE/geocode'))\n", - "\n", - "#Get frames and orbit\n", - "swathList = topsapp._insar.getValidSwathList(topsapp.swaths)\n", - "frames = []\n", - "for swath in swathList:\n", - " referenceProduct = topsapp._insar.loadProduct(os.path.join(topsapp._insar.fineCoregDirname, 'IW{0}.xml'.format(swath)))\n", - " frames.append(referenceProduct)\n", - "\n", - "# Load orbits\n", - "orb = topsapp._insar.getMergedOrbit(frames)\n", - "\n", - "sensing_grid = np.zeros((data['lats'].shape[0],\n", - " data['lons'].shape[0], \n", - " data['heights'].shape[0]))\n", - "\n", - "for iy, lat in enumerate(data['lats'][:]):\n", - " for ix, lon in enumerate(data['lons'][:]):\n", - " for iz, hgt in enumerate(data['heights'][:]):\n", - " pixel_sensing_az_rg = orb.geo2rdr(llh=(lat,lon,hgt)) # approx, should alwys start from 1 swath\n", - " sensing = pixel_sensing_az_rg[0] + timedelta(seconds=(pixel_sensing_az_rg[1]) / SPEED_OF_LIGHT) # azimuth_sensing_time +( ix * PRF + range/speed of light)\n", - " sensing_grid[iy,ix,iz] = sensing.timestamp()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "(lat,lon,hgt) 1/PRF" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "orb.geo2rdr(llh=(30,-115,1050))[0] - orb.geo2rdr(llh=(30,-115,1000))[0] " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "pixel_sensing_az_rg[1] / SPEED_OF_LIGHT" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "frames[1].bursts.burst1.firstValidSample * (1/frames[1].bursts.burst1.prf)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "frames[1].bursts.burst1.lastValidSample * (1/frames[1].bursts.burst1.prf) + pixel_sensing_az_rg[1] / SPEED_OF_LIGHT" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "np.broadcast_to(data['lats'][:], (72,47))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from datetime import datetime" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "datetime.fromtimestamp(sensing_grid[:,:,3])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from datetime import datetime\n", - "[datetime.fromtimestamp(i) for i in sensing_grid[0,:,3]]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.imshow(sensing_grid[:,:,])\n", - "plt.colorbar()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "datetime.fromtimestamp()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "lats = data['lats'][:]\n", - "lons = data['lons'][:]\n", - "heights = data['heights'][:]\n", - "\n", - "# Reshape the arrays to create a 3D grid\n", - "ny, nx, nz = (lats.shape[0], \n", - " lons.shape[0], \n", - " heights.shape[0])\n", - "lats = lats.reshape(ny, nx, 1)\n", - "lons = lons.reshape(ny, nx, 1)\n", - "heights = heights.reshape(1, 1, nz)\n", - "\n", - "\n", - "pixels_sensing_az_rg = orb.geo2rdr(llh=(lats, lons, heights))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "coords = np.stack((data['lats'][:], data['lons'][:], data['heights'][:]), axis=-1)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "topsapp_env", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.16" - }, - "orig_nbformat": 4 - }, - "nbformat": 4, - "nbformat_minor": 2 -} From b490cd7c2635c2e26606d1d9221e4256f8615199 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 24 Apr 2023 09:13:20 -0700 Subject: [PATCH 20/33] import order --- isce2_topsapp/iono_proc.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index 9ea0b7ff..d75f0f0a 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -6,19 +6,18 @@ import os import site from pathlib import Path +import datetime +import shutil from typing import Type, Union -import shutil -import datetime import numpy as np import scipy.signal as ss from isce.applications import topsApp from isce.components import isceobj from isce.components.isceobj.TopsProc import runIon from isce.components.isceobj.TopsProc.runMergeBursts import ( - interpolateDifferentNumberOfLooks, mergeBox, - adjustValidWithLooks, mergeBurstsVirtual, - mergeBursts2, multilook) + adjustValidWithLooks, interpolateDifferentNumberOfLooks, + mergeBox, mergeBursts2, mergeBurstsVirtual, multilook) from numpy.typing import NDArray from osgeo import gdal from skimage import morphology @@ -456,8 +455,8 @@ def mask_interferogram( def brige_components(unwrapped_ifg: str, connected_components: str) -> None: """ - This routine preforms "bridging' of unwrapped phase - connected components. Each component is shifted with + This routine preforms "bridging' of unwrapped phase + connected components. Each component is shifted with its median value unwrapped_ifg : str @@ -882,7 +881,7 @@ def ionSwathBySwath(self: Type[topsApp.TopsInSAR], ######################################################## # slight filtering to improve the estimation accuarcy # of swath difference - if 1 and shutil.which('psfilt1') != None: + if 1 and shutil.which('psfilt1') is not None: cmd1 = 'mv {} tmp'.format(lowerintfile) cmd2 = 'psfilt1 tmp {} {} .3 32 8'.format(lowerintfile, width) cmd3 = 'rm tmp' From 1745ab277885136d132f35f3137bc235c957ab27 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 24 Apr 2023 09:24:25 -0700 Subject: [PATCH 21/33] update import order 2 --- isce2_topsapp/iono_proc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index d75f0f0a..ad61d3b6 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -5,8 +5,8 @@ import multiprocessing import os import site -from pathlib import Path import datetime +from pathlib import Path import shutil from typing import Type, Union From 4f833c41eafc2f367da7a6360f479f3ff2d5558e Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 24 Apr 2023 09:25:30 -0700 Subject: [PATCH 22/33] update order 3 --- isce2_topsapp/iono_proc.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index ad61d3b6..4cb1f050 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -4,10 +4,10 @@ import multiprocessing import os -import site import datetime -from pathlib import Path +import site import shutil +from pathlib import Path from typing import Type, Union import numpy as np From b2a2ff983d4551184970b81c668da2948a5bb9ad Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 24 Apr 2023 09:28:43 -0700 Subject: [PATCH 23/33] update_order3 --- isce2_topsapp/iono_proc.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index 4cb1f050..fd79682d 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -3,10 +3,10 @@ # California Institute of Technology import multiprocessing -import os import datetime -import site +import os import shutil +import site from pathlib import Path from typing import Type, Union From 9b7684e3cca6afa791a8fa2dc2f31e5259672013 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 24 Apr 2023 09:37:57 -0700 Subject: [PATCH 24/33] update_import_order4 --- isce2_topsapp/iono_proc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index fd79682d..a80d45c6 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -2,8 +2,8 @@ # Copyright 2023 # California Institute of Technology -import multiprocessing import datetime +import multiprocessing import os import shutil import site From 6566a9530f5e6b5a421af1dbb21e6357e97ee3a1 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 25 Apr 2023 11:30:34 -0700 Subject: [PATCH 25/33] turn off ESD when using iono computation --- isce2_topsapp/__main__.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/isce2_topsapp/__main__.py b/isce2_topsapp/__main__.py index c03902ad..60924744 100644 --- a/isce2_topsapp/__main__.py +++ b/isce2_topsapp/__main__.py @@ -178,6 +178,10 @@ def gunw_slc(): json.dump(loc_data, open("loc_data.json", "w"), indent=2, cls=MetadataEncoder) + # Turn-off ESD when using ionospheric computation + if args.estimate_ionosphere_delay: + args.esd_coherence_threshold = -1 + topsapp_processing( reference_slc_zips=loc_data["ref_paths"], secondary_slc_zips=loc_data["sec_paths"], From e28266927fcebb5cf80189b3aa1968037ef672e7 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 8 May 2023 16:00:58 -0700 Subject: [PATCH 26/33] added derampinh before bridging --- isce2_topsapp/iono_proc.py | 64 +++++++++++++++++++++++++++++++++++++- 1 file changed, 63 insertions(+), 1 deletion(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index a80d45c6..f34e031a 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -452,6 +452,61 @@ def mask_interferogram( int_array.astype(np.complex64).tofile(ifgFilename) +def deramp(data: NDArray, mask_in: NDArray= None, + ramp_type: str= 'linear') -> NDArray: + ''' + Taken from: https://github.com/insarlab/MintPy + /src/mintpy/objects/ramp.py#L23 + + This function preforms deramping of unwrapped phase + + data : NDarray + 2D array data to be derampped + mask_in : NDarray + 2D np.ndarray, mask of pixels used for ramp estimation + ramp_type : str + name of ramp to be estimated. "linear" or 'quadratic" + + ''' + + dshape = data.shape + length, width = dshape[-2:] + + # for 2d only + data = data.reshape(-1, 1) + dmean = np.array(data).flatten() + + if mask_in is None: + mask_in = np.ones((length, width), dtype=np.float32) + mask = (mask_in != 0).flatten() + del mask_in + + # 2. ignore pixels with NaN and/or zero data value + mask *= ~np.isnan(dmean) + mask *= (dmean != 0.) + del dmean + + # design matrix + xx, yy = np.meshgrid(np.arange(0, width), + np.arange(0, length)) + xx = np.array(xx, dtype=np.float32).reshape(-1, 1) + yy = np.array(yy, dtype=np.float32).reshape(-1, 1) + ones = np.ones(xx.shape, dtype=np.float32) + if ramp_type == 'linear': + G = np.hstack((yy, xx, ones)) + elif ramp_type == 'quadratic': + G = np.hstack((yy**2, xx**2, yy*xx, yy, xx, ones)) + else: + raise ValueError(f'un-recognized ramp type: {ramp_type}') + + # estimate ramp + X = np.dot(np.linalg.pinv(G[mask, :], rcond=1e-15), data[mask, :]) + ramp = np.dot(G, X) + ramp = np.array(ramp, dtype=data.dtype) + + return ramp.reshape(dshape) + + def brige_components(unwrapped_ifg: str, connected_components: str) -> None: """ @@ -482,8 +537,15 @@ def brige_components(unwrapped_ifg: str, # Apply a binary closing operation to the mask mask = morphology.binary_closing(mask) + # Deramp the data, to remove any trend + # before estimating median + print('Deramp') + ramp = deramp(data=interferogram, mask_in=mask, + ramp_type='linear') + derampped_interferogram = interferogram - ramp + # Calculate the median phase value for the current component - median_phase = np.median(interferogram[mask]) + median_phase = np.median(derampped_interferogram[mask]) # Subtract the median phase from the current component # to remove any phase jumps From b52736f42f8e373983ee91cebedd8794f1b1c032 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 8 May 2023 16:02:27 -0700 Subject: [PATCH 27/33] fix falke8 comments --- isce2_topsapp/iono_proc.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index f34e031a..5d235b40 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -452,8 +452,7 @@ def mask_interferogram( int_array.astype(np.complex64).tofile(ifgFilename) -def deramp(data: NDArray, mask_in: NDArray= None, - ramp_type: str= 'linear') -> NDArray: +def deramp(data: NDArray, mask_in: NDArray = None, ramp_type: str = 'linear') -> NDArray: ''' Taken from: https://github.com/insarlab/MintPy /src/mintpy/objects/ramp.py#L23 From 8c78434241695e7ed21f30b0b8b980b4948e77fb Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 8 May 2023 16:03:34 -0700 Subject: [PATCH 28/33] fix #128 --- isce2_topsapp/iono_proc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/isce2_topsapp/iono_proc.py b/isce2_topsapp/iono_proc.py index 5d235b40..e841aa63 100644 --- a/isce2_topsapp/iono_proc.py +++ b/isce2_topsapp/iono_proc.py @@ -487,7 +487,7 @@ def deramp(data: NDArray, mask_in: NDArray = None, ramp_type: str = 'linear') -> # design matrix xx, yy = np.meshgrid(np.arange(0, width), - np.arange(0, length)) + np.arange(0, length)) xx = np.array(xx, dtype=np.float32).reshape(-1, 1) yy = np.array(yy, dtype=np.float32).reshape(-1, 1) ones = np.ones(xx.shape, dtype=np.float32) From 928f2ea818e2b7fdbe04d476df29100cf2465308 Mon Sep 17 00:00:00 2001 From: Charlie Marshak Date: Mon, 22 May 2023 12:22:31 -0700 Subject: [PATCH 29/33] fix slc tests --- isce2_topsapp/localize_slc.py | 14 +++++++++++++- tests/test_localize_slc.py | 12 +++++++++--- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/isce2_topsapp/localize_slc.py b/isce2_topsapp/localize_slc.py index 3adf2ec7..ad0de520 100644 --- a/isce2_topsapp/localize_slc.py +++ b/isce2_topsapp/localize_slc.py @@ -1,10 +1,13 @@ +import io import netrc from concurrent.futures import ThreadPoolExecutor +from functools import lru_cache from pathlib import Path from warnings import warn import asf_search as asf import geopandas as gpd +import requests from dateparser import parse from shapely.geometry import GeometryCollection, Polygon, shape from shapely.ops import unary_union @@ -120,9 +123,18 @@ def check_track_numbers(slc_properties: list): return False +@lru_cache(maxsize=None) +def get_world_df() -> gpd.GeoDataFrame: + natural_earth_url = ('https://www.naturalearthdata.com/' + 'http//www.naturalearthdata.com/download/10m/physical/ne_10m_land.zip') + resp = requests.get(natural_earth_url, headers={"User-Agent": "XY"}) + df_world = gpd.read_file(io.BytesIO(resp.content)) + return df_world + + def get_percent_water_from_ne_land(ifg_geo: Polygon): """Gets percent_water using Natural Earth Low Res Mask""" - df_world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres')) + df_world = get_world_df() world_geo = df_world.geometry.unary_union land_overlap = world_geo.intersection(ifg_geo) return (1 - land_overlap.area / ifg_geo.area) * 100 diff --git a/tests/test_localize_slc.py b/tests/test_localize_slc.py index be50b0ae..b1cfa562 100644 --- a/tests/test_localize_slc.py +++ b/tests/test_localize_slc.py @@ -67,20 +67,26 @@ def test_bad_tracks_with_same_flight_direction(): def test_warnings_over_water(): + wtr_msg = 'If there are not enough bursts over land - ISCE2 will fail.' + # Intersection over water ref_ids = ['S1B_IW_SLC__1SDV_20211210T153506_20211210T153533_029965_0393C6_D840'] sec_ids = ['S1A_IW_SLC__1SDV_20211122T153535_20211122T153602_040686_04D3DB_520A'] - with pytest.warns(RuntimeWarning): + with pytest.warns(RuntimeWarning) as records: download_slcs(ref_ids, sec_ids, frame_id=-1, dry_run=True) + assert any([wtr_msg in str(r.message) for r in records]) # Tibet (no water) ref_ids = ['S1A_IW_SLC__1SDV_20170817T120001_20170817T120028_017963_01E230_A23A'] sec_ids = ['S1A_IW_SLC__1SSV_20160717T115946_20160717T120014_012188_012E84_F684', 'S1A_IW_SLC__1SSV_20160717T120012_20160717T120039_012188_012E84_4198'] - with warnings.catch_warnings(): - warnings.simplefilter("error") + + # Make sure water warning is not thrown + with pytest.warns(Warning) as records: download_slcs(ref_ids, sec_ids, frame_id=-1, dry_run=True) + if records: + assert all([wtr_msg not in str(r.message) for r in records]) def test_bad_date_order(): From 30bcc9635ebf8afeb815be526ddd62f12ce2336c Mon Sep 17 00:00:00 2001 From: Charlie Marshak Date: Mon, 22 May 2023 13:34:57 -0700 Subject: [PATCH 30/33] fix warning catch --- tests/test_localize_slc.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_localize_slc.py b/tests/test_localize_slc.py index b1cfa562..5d55a063 100644 --- a/tests/test_localize_slc.py +++ b/tests/test_localize_slc.py @@ -73,9 +73,9 @@ def test_warnings_over_water(): ref_ids = ['S1B_IW_SLC__1SDV_20211210T153506_20211210T153533_029965_0393C6_D840'] sec_ids = ['S1A_IW_SLC__1SDV_20211122T153535_20211122T153602_040686_04D3DB_520A'] - with pytest.warns(RuntimeWarning) as records: + with pytest.warns(RuntimeWarning) as warning_records: download_slcs(ref_ids, sec_ids, frame_id=-1, dry_run=True) - assert any([wtr_msg in str(r.message) for r in records]) + assert any([wtr_msg in str(r.message) for r in warning_records]) # Tibet (no water) ref_ids = ['S1A_IW_SLC__1SDV_20170817T120001_20170817T120028_017963_01E230_A23A'] @@ -83,10 +83,10 @@ def test_warnings_over_water(): 'S1A_IW_SLC__1SSV_20160717T120012_20160717T120039_012188_012E84_4198'] # Make sure water warning is not thrown - with pytest.warns(Warning) as records: + with warnings.catch_warnings(record=True) as warning_records: download_slcs(ref_ids, sec_ids, frame_id=-1, dry_run=True) - if records: - assert all([wtr_msg not in str(r.message) for r in records]) + if warning_records: + assert all([wtr_msg not in str(r.message) for r in warning_records]) def test_bad_date_order(): From 6b8d5a0854116f014d2bfd708cbc48e800fb55ef Mon Sep 17 00:00:00 2001 From: Charlie Marshak Date: Mon, 22 May 2023 13:45:41 -0700 Subject: [PATCH 31/33] changelog --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index e0ca5745..6d5b255e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,8 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). * For Solid Earth Tide computation, derive coordinates and spacing from geotrans as opposed to latitude/longitude metadata arrays * Include topsapp_iono template. * Increases DEM buffer to .4 from .1 to ensure the extent of at least two bursts (~40 km) are added when retrieving DEM (because estimated footprint can differ from what ISCE2 generates for a GUNW extent) +* Catch warnings in tests and match messages to ensure package warnings do not fail test suite +* Read low resolution Natural Earth land masses from public url due to removal from geopandas package. ### Added * localize_data within __main__.py added option to use/not use water mask for ionosphere processing From f16dbf500e0bcf1f6b15025c79e4295e376b07c8 Mon Sep 17 00:00:00 2001 From: Charlie Marshak Date: Mon, 22 May 2023 13:45:52 -0700 Subject: [PATCH 32/33] put records outside of warning context --- tests/test_localize_slc.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_localize_slc.py b/tests/test_localize_slc.py index 5d55a063..e1cad719 100644 --- a/tests/test_localize_slc.py +++ b/tests/test_localize_slc.py @@ -85,8 +85,8 @@ def test_warnings_over_water(): # Make sure water warning is not thrown with warnings.catch_warnings(record=True) as warning_records: download_slcs(ref_ids, sec_ids, frame_id=-1, dry_run=True) - if warning_records: - assert all([wtr_msg not in str(r.message) for r in warning_records]) + if warning_records: + assert all([wtr_msg not in str(r.message) for r in warning_records]) def test_bad_date_order(): From bbf24961bfe3225e53a14e301842812e31a0b92e Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 23 May 2023 14:46:46 -0700 Subject: [PATCH 33/33] update changelog --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5b9b9452..67dd74d8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,7 +14,7 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). * Increases DEM buffer to .4 from .1 to ensure the extent of at least two bursts (~40 km) are added when retrieving DEM (because estimated footprint can differ from what ISCE2 generates for a GUNW extent) * Catch warnings in tests and match messages to ensure package warnings do not fail test suite * Read low resolution Natural Earth land masses from public url due to removal from geopandas package. - +* For ionosphere computation over water, includes masking conncomp zero, phase bridging, and modified adaptive gaussian filtering ### Added * localize_data within __main__.py added option to use/not use water mask for ionosphere processing