diff.py btw. time-series and velocity + other changes

[yuankailiu]

Description of proposed changes

Several modifications.

1. diff.py: does (timeseries) - (velocity) = (timeseries_new)

   • implemented by function diff_timeseries_velocity()
   • this commit is motivated by doing correction based on inputs/ITRF.h5 velocity model in the time-series domain. Thus, the commit only allows a time function model m to identify linear, annual, and semiannual functions. Later, we can write a new ts2velo.hdf5_dataset2model() to produce this m based on the dataset_list names in general velocity files.

2. spatial_filter.py: add a median filter from skimage.filters.median

3. image_stitch.py: allow for vmin/vmax, cmap, and display unit

4. plot.py: numTriNonzeroIntAmbiguity.h5 mask out the non-zero values (non-zero is bad pixels) for mask creating

5. gps.py object: GPS time function fit with a user-defined model

Reminders

• Pass Pre-commit check (green)
• Pass Codacy code review (green)
• Pass Circle CI test (green)
• Make sure that your code follows our style. Use the other functions/files as a basis.
• If modifying functionality, describe changes to function behavior and arguments in a comment below the function declaration.
• If adding new functionality, add a detailed description to the documentation and/or an example.

[yunjunz]

Thank you @yuankailiu for the useful PR. For the diff.py time-series and velocity functionality, I modified the code on top of yours. A test based on the SanFranSenDT42 example dataset between timeseries.h5 and velocity.h5 shows a reasonable result, but I don't have a large dataset at hand, could you re-confirm the effectiveness of plate motion correction?

For the other changes, I left a few minor comments below. Take a look and let me know what you think please.

[yuankailiu]

Here is the demonstraion of the correction done in the time-series domain. It is effective.

# original time series (with SET, ERA5, Ion, and topo residual corrections) fitted with time functions 
timeseries2velocity.py timeseries_SET_ERA5_Ion_demErr.h5 --poly-order 1 --periodic 1.0 0.5 -o velocity_test1.h5 --ref-lalo 29.5463 36.0810 --ref-date 20190118

# Difference that with the ITRF14 velocity file of the Arabian plate (i.e., apply plate motion correction)
diff.py timeseries_SET_ERA5_Ion_demErr.h5 inputs/ITRF14.h5

# Fit the new time-series file `timeseries_SET_ERA5_Ion_demErr_diff_ITRF14.h5` with time functions
timeseries2velocity.py timeseries_SET_ERA5_Ion_demErr_diff_ITRF14.h5 --poly-order 1 --periodic 1.0 0.5 -o velocity_test1_ITRF.h5 --ref-lalo 29.5463 36.0810 --ref-date 20190118

# plot the original veloctiy
view.py velocity_test1.h5 velocity -c RdYlBu_r --dem ./inputs/srtm.dem --alpha 0.8 --dem-nocontour --shade-exag 0.02 --shade-min -6000 --shade-max 4000 --mask ./maskTempCoh_0.9.h5 --unit mm --ref-lalo 29.5463 36.0810 --vlim -5 5 --nodisplay --dpi 300 --figtitle no_ITRF_cor

# plot the plate motion corrected velocity
view.py velocity_test1_ITRF.h5 velocity -c RdYlBu_r --dem ./inputs/srtm.dem --alpha 0.8 --dem-nocontour --shade-exag 0.02 --shade-min -6000 --shade-max 4000 --mask ./maskTempCoh_0.9.h5 --unit mm --ref-lalo 29.5463 36.0810 --vlim -5 5 --nodisplay --dpi 300 --figtitle with_ITRF_cor

# plot the ITRF14 plate motion velocity itself
view.py inputs/ITRF14.h5 velocity -c RdYlBu_r --dem ./inputs/srtm.dem --alpha 0.8 --dem-nocontour --shade-exag 0.02 --shade-min -6000 --shade-max 4000 --mask ./maskTempCoh_0.9.h5 --unit mm --ref-lalo 29.5463 36.0810 --vlim -5 5 --nodisplay --dpi 300 --figtitle ITRF_cor

[yunjunz]

Looks all good to me.

[yuankailiu]

Looks all good to me.

Thanks! Appreciate for reviewing and refactoring them! I will merge it.

[scottstanie]

I saw that you added the nice way to ignore nans here, but that it's only for lowpass gaussian right now. Sorry I didn't see this in time, but what if we changed this

MintPy/src/mintpy/spatial_filter.py

Lines 68 to 70 in 6fffa3b

	elif filter_type == "highpass_gaussian":
	lp_data = filters.gaussian(data, sigma=filter_par)
	data_filt = data - lp_data

to be

    elif filter_type == "highpass_gaussian":
        return data - filter_data(data, "lowpass_gaussian", filter_par=filter_par):

so we don't have to redo the logic twice?

[yunjunz]

Good catch @scottstanie. I missed it as well. Could you issue a PR for it?