Resampling geocoded data

[rzinke]

I need to resample time-series data that are already geocoded. That is, I need to crop/expand map bounds and change the pixel size. Setting the geographic bounds from earlier in the workflow (e.g., during ISCE or ARIA-tools processing steps) is not an option -- I need to work directly with the HDF5 MintPy outputs. From what I can tell, subset.py is not capable of changing the pixel size. I see that geocode.py calls on the mintpy.objects.resample class. Although this class is not ready to resample geocoded data, there are placeholders for potential future functionality.

If I were to build out the geocoded-to-geocoded resampling methods for the resample class, could I get some guidance and feedback during the process? For example, my initial questions are:

• Do the pixel coordinates define the upper-left corners of the pixels, or the pixel centers?
• What other pitfalls should I avoid that could introduce a sub-pixel shift?
• Is Pyresample the preferred method of resampling the geocoded data?
• For scipy.interpolate RGI, how is the best way to define a new grid (e.g., using np.arange, np.linspace, etc)?

Please let me know if this is something that MintPy is interested in. If so, I would appreciate any help I can get based on community knowledge and lessons learned to help bring me up to speed, as time and funding are limited as always.

Thank you for any thoughts.

[yunjunz]

Hi @rzinke, resampling the geocoded data/matrix (a special case of resampling) is easy to do using the skimage.transform.resize function, which is much faster than pyresample and scipy.interpolate. You could find many example usage of this in mintpy, such as this one:

MintPy/src/mintpy/plate_motion.py

Lines 102 to 108 in 40c0d90

	# resample coarse grid back to the initial fine grid
	print(f'resample plate motion from corase back to original grid: {lats.shape} -> {shape_geo}'
	' via skimage.transform.resize ...')
	kwargs = dict(order=1, mode='edge', anti_aliasing=True, preserve_range=True)
	ve = resize(ve_low, shape_geo, **kwargs)
	vn = resize(vn_low, shape_geo, **kwargs)
	vu = resize(vu_low, shape_geo, **kwargs)

The Y/X_FIRST represents the upper left corner of the upper left pixel, as defined in https://mintpy.readthedocs.io/en/latest/api/attributes/#required_attributes, so they do not need to change during this process, following the philosophy of "a pixel is an area".

If you want to crop/expand the map extent, I would suggest to resize first to the desired pixel size, then run subset.py, which support both crop/expand operations.

[rzinke]

Thanks for the thoughts and suggestion! I'll give it a try.