[Help]: PS processing

[SteffanDavies]

I am having issues with sbas.ps_parallel, some indications on workflow would be helpful.

Is ps_parallel supposed to be run after intf_parallel and merge_parallel? Doing so causes {final date}_F{merged subswaths}.PRM missing file error.

FileNotFoundError: [Errno 2] No such file or directory: 'raw_asc/S1_20230530_ALL_F23.PRM'

Removing this file from dates and passing custom dates causes another error:

sbas.ps_parallel(dates=sbas.df.index[:-1])
Exception: "ValueError('mmap length is greater than file size')"

[AlexeyPechnikov]

It is not ready for public usage. I share the examples on Patreon for early access for sponsors only. For now, PS-SBAS works but we need ability to export the original resolution (~15x4 meter) grids for PS validation and so on. Also, there are two PS-SBAS approaches already available in PyGMTSAR and a lot of testing required to provide the recommendations how and when use them.

[AlexeyPechnikov]

@SteffanDavies See the example of PS+SBAS processing: https://www.linkedin.com/feed/update/urn:li:activity:7089507310389121024/

[SteffanDavies]

@SteffanDavies See the example of PS+SBAS processing: https://www.linkedin.com/feed/update/urn:li:activity:7089507310389121024/

I am currently testing this and in your example sbas.ps_parallel() is called before merging, however doing so will cause:

AssertionError: ERROR: multiple subswaths [2 3] found, merge them first using SBAS.merge_parallel()

So I suppose this only supports 1 subswath at the moment?

[AlexeyPechnikov]

It should work for all the subswaths. Are you sure SBAS.ps_parallel() does not work for multiple subswaths?

[SteffanDavies]

It should work for all the subswaths. Are you sure SBAS.ps_parallel() does not work for multiple subswaths?

If I change to only 1 subswath it works fine, otherwise it gives the above error.

Also ps_parallel() leads to huge unmanaged memory used even after it has finished. My 64 GB RAM and 16GB Swap were full after it finished.

[SteffanDavies]

More issues, topo_ra error and empty interferograms:

Dem is fine:

topo_ra error seems dask related?

Empty intf and corr (only plot 1 but all the same):

[AlexeyPechnikov]

For some reason, you have no active Dask client for the processing and most of PyGMTSAR operations are not available.

[SteffanDavies]

Dask is running, I've never had this error before.

[AlexeyPechnikov]

The Dask issue can be related to disconnected client or maybe Dask version change (when Dask or distributed libraries upgraded and inconsistent to the running cluster).

ps_parallel() leads to huge unmanaged memory used even after it has finished.

Please explain how much scenes do you have to analyse?

[SteffanDavies]

I reverted to PyGMTSAR pygmtsar-2023.5.3 and it is working. Something has broken in latest code.

[SteffanDavies]

The Dask issue can be related to disconnected client or maybe Dask version change (when Dask or distributed libraries upgraded and inconsistent to the running cluster).

ps_parallel() leads to huge unmanaged memory used even after it has finished.

Please explain how much scenes do you have to analyse?

181 scenes.

[AlexeyPechnikov]

But this is the recent PyPI version for now. GitHub version is in the alpha stage and I change some function calls to allow a single-pixel accuracy for PS geocoding and more. For SBAS analysis, we don’t need ~4m output accuracy, but it’s critical to verify detected PS.

[SteffanDavies]

But this is the recent PyPI version for now. GitHub version is in the alpha stage and I change some function calls to allow a single-pixel accuracy for PS geocoding and more. For SBAS analysis, we don’t need ~4m output accuracy, but it’s critical to verify detected PS.

I think the problem is related to your new dask delayed code. I was testing the new PS-SBAS to compare differences over railway.

[AlexeyPechnikov]

Yes, actually, I've fixed lots of issues in the github version last days. And you are right, ps_parallel does not work for multiple subswaths (it's easy to fix). About the memory consumption I need to test some ways — while it looks simple for sbas.ps_parallel() when the results should be stored on dist, it's not so trivial for sbas.ps_parallel(interactive=True) call when we expect to get the lazy result immediately.

[SteffanDavies]

Yes, actually, I've fixed lots of issues in the github version last days. And you are right, ps_parallel does not work for multiple subswaths (it's easy to fix). About the memory consumption I need to test some ways — while it looks simple for sbas.ps_parallel() when the results should be stored on dist, it's not so trivial for sbas.ps_parallel(interactive=True) call when we expect to get the lazy result immediately.

Glad I could help. Looking forward to new PS processing.

[AlexeyPechnikov]

You'd try functions https://github.com/mobigroup/gmtsar/blob/pygmtsar/pygmtsar/pygmtsar/SBAS_ps.py on top of the same commit as the notebook above (the recent commits have changed SBAS and STL processing). It works well for subswaths while the old intf_parallel "weight" supports only a single subswath grid. 200 scenes ADI calculation requires ~30 minutes for a single subswath on Apple iMac M1 or Apple Air M2.

[AlexeyPechnikov]

Some posts about known corner reflectors detection on Sentinel-1 amplitudes and InSAR coherence: https://www.linkedin.com/feed/update/urn:li:activity:7091290249380691968/ https://www.linkedin.com/feed/update/urn:li:share:7091295824273408001/

The normed dispersion is large for high amplitude pixels. The intensities-based Amplitude Dispersion Index (ADI) looks as

when the amplitude-based ADI is not usable:

[AlexeyPechnikov]

The weighted interferometry results validation using known corner reflectors: https://www.linkedin.com/feed/update/urn:li:activity:7092021658172981248/

[AlexeyPechnikov]

Checked and it works for the known corner reflector: https://www.linkedin.com/feed/update/urn:li:share:7098256457971732480/

[AlexeyPechnikov]

A set of screencasts available on the YouTube channel: https://m.youtube.com/channel/UCSEeXKAn9f_bDiTjT6l87Lg See the related example notebooks on Patreon https://www.patreon.com/pechnikov and on LinkedIn in PDF.

[AlexeyPechnikov]

@SteffanDavies GMTSAR subswath merging is not single-pixel accurate, and it introduces inconsistencies in PSI processing. Although we may not detect the shift in output products, the azimuthal accuracy of ±15m makes it impossible to perform single-pixel time-tracing. For the common multi-looking output of 60m or 90m, this is not an issue.

[SteffanDavies]

@SteffanDavies GMTSAR subswath merging is not single-pixel accurate, and it introduces inconsistencies in PSI processing. Although we may not detect the shift in output products, the azimuthal accuracy of ±15m makes it impossible to perform single-pixel time-tracing. For the common multi-looking output of 60m or 90m, this is not an issue.

The offset seems semi-consistent, wouldn't it be possible to do some kind of pixel shift in azimuth?

[AlexeyPechnikov]

GMTSAR calculates the offsets independently for every scene, and the accuracy is ±1 pixel for my test case. We cannot predict the shift, but it can be resolved using reference scene vertical shifts for all the scenes. For more details, please refer to my LinkedIn post: https://www.linkedin.com/feed/update/urn:li:activity:7110177621832843264/

[AlexeyPechnikov]

@SteffanDavies You can check the recent commit 030325a1983fad9394d564a88da78e0f78918edc in PyGMTSAR2 for performance evaluation on your hardware using the attached notebook: phasediff.ipynb.zip Please share your processed notebook in PDF (printable) format for reviewing the progress bars.

It should be significantly faster when using a processing chunk size of 2048 and with NetCDF compression disabled. The old configuration was tuned for a chunk size of 512 or 1024, and NetCDF compression was always enabled. Potentially, using Dask chunks that are four times larger allows Dask to work effectively on hosts with 4x4 times more cores compared to my 8-core development host, making it effective on 128+ core hosts.

P.S. The DEM is intentionally incomplete (check the upper left corner on the plot below) to verify that PyGMTSAR2 has no issues in this case, unlike GMTSAR where a single NODATA pixel on DEM crashes the processing completely.

[SteffanDavies]

@SteffanDavies You can check the recent commit 030325a in PyGMTSAR2 for performance evaluation on your hardware using the attached notebook: phasediff.ipynb.zip Please share your processed notebook in PDF (printable) format for reviewing the progress bars.

It should be significantly faster when using a processing chunk size of 2048 and with NetCDF compression disabled. The old configuration was tuned for a chunk size of 512 or 1024, and NetCDF compression was always enabled. Potentially, using Dask chunks that are four times larger allows Dask to work effectively on hosts with 4x4 times more cores compared to my 8-core development host, making it effective on 128+ core hosts.

P.S. The DEM is intentionally incomplete (check the upper left corner on the plot below) to verify that PyGMTSAR2 has no issues in this case, unlike GMTSAR where a single NODATA pixel on DEM crashes the processing completely.

Is there a specific list of scenes you would like me to test on?

[AlexeyPechnikov]

You can choose from any scenes.

[AlexeyPechnikov]

@SteffanDavies Maybe do you have some results? I’ve shared my timings for GMTSAR example: https://www.linkedin.com/feed/update/urn:li:activity:7117094805393858560/

[AlexeyPechnikov]

New PyGMTSAR provides processing that is 10 times faster than GMTSAR on the same case, S1A_Stack_CPGF_T173. This notebook is already available on Google Colab: https://www.linkedin.com/feed/update/urn:li:activity:7117369530749730816/

[SteffanDavies]

I was having some issues setting up but I noticed you no longer import SBAS from pygmtsar

[SteffanDavies]

Help command for S1.scan_slc() function listing "basedir" as argument despite being unrecognized.

[SteffanDavies]

sbas.align() ValueError: no path specified:

`--------------------------------------------------------------------------- _RemoteTraceback Traceback (most recent call last) _RemoteTraceback: """ Traceback (most recent call last): File "/home/steffan/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/joblib/externals/loky/process_executor.py", line 463, in _process_worker r = call_item() File "/home/steffan/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/joblib/externals/loky/process_executor.py", line 291, in call return self.fn(*self.args, self.kwargs) File "/home/steffan/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/joblib/parallel.py", line 589, in call return [func(*args, *kwargs) File "/home/steffan/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/joblib/parallel.py", line 589, in return [func(args, kwargs) File "/home/steffan/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/pygmtsar/Stack_align.py", line 121, in align_ref_subswath self.make_s1a_tops(subswath, debug=debug) File "/home/steffan/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/pygmtsar/Stack_reframe_gmtsar.py", line 127, in make_s1a_tops self.ext_orb_s1a(subswath, stem, date, debug=debug) File "/home/steffan/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/pygmtsar/Stack_reframe_gmtsar.py", line 47, in ext_orb_s1a orbit = os.path.relpath(df['orbitpath'].iloc[0], self.basedir) File "/home/steffan/miniconda3/envs/pygmtsar2/lib/python3.9/posixpath.py", line 454, in relpath raise ValueError("no path specified") ValueError: no path specified """

The above exception was the direct cause of the following exception:

ValueError Traceback (most recent call last) Cell In[64], line 1 ----> 1 sbas.align()

File ~/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/pygmtsar/Stack_align.py:832, in Stack_align.align(self, dates, n_jobs, debug) 829 # prepare reference scene 830 #self.stack_ref() 831 with self.tqdm_joblib(tqdm(desc='Aligning Reference', total=len(subswaths))) as progress_bar: --> 832 joblib.Parallel(n_jobs=n_jobs)(joblib.delayed(self.align_ref_subswath)(subswath, debug=debug) for subswath in subswaths) 834 # prepare secondary images 835 with self.tqdm_joblib(tqdm(desc='Aligning Repeat', total=len(dates_rep)*len(subswaths))) as progress_bar:

File ~/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/joblib/parallel.py:1952, in Parallel.call(self, iterable) 1946 # The first item from the output is blank, but it makes the interpreter 1947 # progress until it enters the Try/Except block of the generator and 1948 # reach the first yield statement. This starts the aynchronous 1949 # dispatch of the tasks to the workers. 1950 next(output) -> 1952 return output if self.return_generator else list(output)

File ~/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/joblib/parallel.py:1595, in Parallel._get_outputs(self, iterator, pre_dispatch) 1592 yield 1594 with self._backend.retrieval_context(): -> 1595 yield from self._retrieve() 1597 except GeneratorExit: 1598 # The generator has been garbage collected before being fully 1599 # consumed. This aborts the remaining tasks if possible and warn 1600 # the user if necessary. 1601 self._exception = True

File ~/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/joblib/parallel.py:1699, in Parallel._retrieve(self) 1692 while self._wait_retrieval(): 1693 1694 # If the callback thread of a worker has signaled that its task 1695 # triggered an exception, or if the retrieval loop has raised an 1696 # exception (e.g. GeneratorExit), exit the loop and surface the 1697 # worker traceback. 1698 if self._aborting: -> 1699 self._raise_error_fast() 1700 break 1702 # If the next job is not ready for retrieval yet, we just wait for 1703 # async callbacks to progress.

File ~/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/joblib/parallel.py:1734, in Parallel._raise_error_fast(self) 1730 # If this error job exists, immediatly raise the error by 1731 # calling get_result. This job might not exists if abort has been 1732 # called directly or if the generator is gc'ed. 1733 if error_job is not None: -> 1734 error_job.get_result(self.timeout)

File ~/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/joblib/parallel.py:736, in BatchCompletionCallBack.get_result(self, timeout) 730 backend = self.parallel._backend 732 if backend.supports_retrieve_callback: 733 # We assume that the result has already been retrieved by the 734 # callback thread, and is stored internally. It's just waiting to 735 # be returned. --> 736 return self._return_or_raise() 738 # For other backends, the main thread needs to run the retrieval step. 739 try:

File ~/miniconda3/envs/pygmtsar2/lib/python3.9/site-packages/joblib/parallel.py:754, in BatchCompletionCallBack._return_or_raise(self) 752 try: 753 if self.status == TASK_ERROR: --> 754 raise self._result 755 return self._result 756 finally:

ValueError: no path specified`

[AlexeyPechnikov]

Help command for S1.scan_slc() function listing "basedir" as argument despite being unrecognized.

Thanks, the docstring is fixed.

sbas.align() ValueError: no path specified:

It appears that orbit files are missing. Nevertheless, this helps in checking the stack dataframe representation for any missing paths.

[SteffanDavies]

sbas.align() ValueError: no path specified:

It appears that orbit files are missing. Nevertheless, this helps in checking the stack dataframe representation for any missing paths.

I don't think this is the case, It works using only 1 Subswath. Could you please test 3 subswath on your system?

[SteffanDavies]

Stack.restore() function looks for a "Stack.pickle" file, however the file generated is called "stack.pickle"

[SteffanDavies]

@mobigroup I am having some issues advancing with PyGMTSAR2. 1 year stack (28 SLC) with 1 Subswath is taking far too much RAM in interferogram processing:

What do you suggest, considering this is a relatively SMALL test dataset?

PS: I have tried to reduce chunksize and increase swap, it makes no difference:

[AlexeyPechnikov]

function looks for a "Stack.pickle" file

Thanks, fixed.

ValueError: no path specified`

I think I found and fixed the issue, see 98b578a7e55b33cddb04340e50075e5df8bf1768 Some Sentinel-1 library versions have issues with dates (resulting in a 1-day offset), which can cause downloaded orbits not to match the scenes.

I am having some issues advancing with PyGMTSAR2

Please be patient. As I mentioned before, I'm not ready to share large processing examples this year. For now, I'm replacing the previous examples with new ones that are clearer and easier to read to ensure that the new PyGMTSAR works everywhere.

Regardless, processing a single subswath for 200 interferograms can be completed on an Apple iMac M1 (16GB RAM, 4 fast cores + 4 slow cores) in about two hours. However, it's worth noting that iterative processing for 30-50 interferograms is significantly faster. We can employ iterative processing in a loop. I hope that recent changes in the Python language interpreter will make it possible to process hundreds of interferograms without any workarounds within a year. However, at the moment, we still need some straightforward wrapper code for stacks larger than the usual 100-200 interferograms.

[AlexeyPechnikov]

@SteffanDavies

You can try this function for interferograms stack processing:

sbas.intf_multilook_save_stack(baseline_pairs, 'intf90m', resolution=90.)
intf90m = sbas.open_stack('intf90m')
intf90m

It takes 3 hours to process 1059 interferograms when the size of S1_...._ALL_F....grd files is 531MB (about a half of a single subswath) on Apple Air.

[SteffanDavies]

@mobigroup Thanks I will try this

[AlexeyPechnikov]

@SteffanDavies Please use the latest commit in the 'pygmtsar2' branch and ensure all related libraries are upgraded. This is important. The upcoming commits will introduce incompatible changes.

[AlexeyPechnikov]

@SteffanDavies Some updates, maybe? By the way, I have shared the Dask scheduler configuring recommendations and the new interferogram processing explanation: https://www.patreon.com/posts/new-pygmtsar-91285612

[SteffanDavies]

@SteffanDavies Some updates, maybe? By the way, I have shared the Dask scheduler configuring recommendations and the new interferogram processing explanation: https://www.patreon.com/posts/new-pygmtsar-91285612

compute_interferogram_multilook works well, I tested this on 32GB laptop running 2 workers on a 2 year stack.

So I have been reading the code and your latest posts etc. Am I correct in saying that for combined PS / SBAS processing, PS is used as input weight for Gaussian Filter? If this is the case, the latest code, including compute_interferogram_multilook does not allow for using weights in interferogram generation.

[AlexeyPechnikov]

... including compute_interferogram_multilook does not allow for using weights in interferogram generation.

Yes, you are right. You can add this argument on your side for now. There are several methods to use the pixel stability measure, and weighted interferogram processing is the primary one.

By the way, I've added Alaska Satellite Facility (ASF) downloading code to extract the necessary polarizations and subswaths without the need for downloading the full scenes locally. When the old approach takes 8 minutes, the new method takes less than 2 minutes and requires less disk space. For more details, please see this notebook: https://colab.research.google.com/drive/1d9RcqBmWIKQDEwJYo8Dh6M4tMjJtvseC?usp=sharing

[SteffanDavies]

Converting from ra2ll and then ll2ra returning different size grid from original (off by 1):

[AlexeyPechnikov]

At the end of the day, I have conducted a comparison between high-resolution SBAS and PSI results, and they align with expectations. For more details, see my LinkedIn post: https://www.linkedin.com/feed/update/urn:li:activity:7132407659806248960/

[AlexeyPechnikov]

@SteffanDavies Have you completed your roads project? It should be nice using the PS-SBAS + PSI technique.

[SteffanDavies]

@SteffanDavies Have you completed your roads project? It should be nice using the PS-SBAS + PSI technique.

Did you manage to reproduce this problem?

Converting from ra2ll and then ll2ra returning different size grid from original (off by 1):

[AlexeyPechnikov]

@SteffanDavies Yes, but I found I had not committed the changes. Thanks for the reminder, you can try the update now.

[SteffanDavies]

@SteffanDavies Yes, but I found I had not committed the changes. Thanks for the reminder, you can try the update now.

Please update intf_multilook_save_stack() function to receive weight parameter. In the meantime I will add this to the local code but it is missing as an argument to pass to the multilooking function.

Otherwise there will be memory issues.

[SteffanDavies]

@SteffanDavies Yes, but I found I had not committed the changes. Thanks for the reminder, you can try the update now.

This problem persists in the most recent Pygmtsar2 docker container. It makes using OSM impossible due to transformations between radar and ll. Was it committed in a more recent version than the Docker container?

[AlexeyPechnikov]

Please use xr.reindex_like() or xr.interp_like() to unify the grids.