add los2slope to project LOS displacement to the downslope component
#1171
Conversation
|
💖 Thanks for opening this pull request! Please check out our contributing guidelines. 💖 |
|
@taliboliver @mgovorcin Could you take a look at the PR when you get a chance please? Here are a few things to check after a quick look:
@ashokdahal Thank you for the PR! Here is a guide for adding a new executable script in mintpy. I believe a few things are missing in this PR. |
yunjunz
changed the title
los2slope to project LOS displacement to the downslope component
Hi @yunjunz sounds good, I'll try to take a look at it in the coming days. |
There was a problem hiding this comment.
@ashokdahal thanks of the PR. I would suggest to double check the projection equation (please see my comment).
When I run the sample data, I am getting the downslope value at the center of the landslide [y: 188, x:188] as -20mm/yr derived from the -14mm/yr LOS value. I would expect it to be smaller than dlos, once projected.
| return filtered_dem | ||
|
|
||
|
|
||
| def smooth_dem(dem, sigma): |
There was a problem hiding this comment.
I do not see the need to wrap smoothing in the function, as it is directly called from the module, so suggest to remove this function smooth_dem, and replace LINE 308, with
# Apply Gaussian smoothing
smoothed_dem = gaussian_filter(dem, sigma=sigma)
| cn = np.sin(np.deg2rad(aspect)) * np.cos(np.deg2rad(slope)) | ||
| cu = np.cos(np.deg2rad(slope)) | ||
|
|
||
| G = ve * ce + vu * cu + vn * cn |
There was a problem hiding this comment.
Here, 3 components of the unit_vector are sum up, I think projection to slope should be done differently:
on the example for East_West:
hv_proj. * vt_proj
ve = -sin(azimuth_angle) * sin(inc_angle),
you would use above to project LOS to EW and assume vn == vu =0. Hv_proj projects the vector onto EW direction, and vt_proj projects in on the horizontal plane. So I think you should project hv_proj to the slope aspect direction, and vt_proj down slope: e.g hv_proj = cos(azimuth + aspect) * -1 , vt_proj= cos(slope + incidence), and downslope = dlos * hv_proj * vt_proj. Please double check the math
There was a problem hiding this comment.
@mgovorcin The assumption that Ashok is using is different than projecting LOS first to horizontal and then to the downslope direction. The assumption here is that all of the displacement is parallel to the downslope direction and we are only seeing a part of that in the LOS direction, so the result is larger than the LOS value.
There was a problem hiding this comment.
@EJFielding thanks for the explanation, I would suggest to add the assumptions to the function description
There was a problem hiding this comment.
I think this could be tested using a simplified example that can be added to the tests directory as @yunjunz suggested. I made a rough example of this using asc_desc2horz_vert.py as reference. (example attached)
los2slope_test.py.txt
| return G | ||
|
|
||
|
|
||
| def los2slope(dlos, los_inc_angle, los_az_angle, slope, aspect, inps, step=20): |
There was a problem hiding this comment.
suggest to replace inps with the input parameters:
slp_thresh=5.0, asp_thresh=0.0, scaling_factor_thresh=10, step=20
to make this function modular
There was a problem hiding this comment.
I second this suggestion
| dslope = np.zeros((length, width), dtype=np.float32) * np.nan | ||
| scale_factor = np.zeros((length, width), dtype=np.float32) * np.nan | ||
|
|
||
| # 2D incidence / azimuth angle --> invert [window-by-window to speed up] |
There was a problem hiding this comment.
do not see the need to use for loop, and do projection window-by-window as the final output can be done with matrix multiplication element-wise... this is more or the decomposition approach
| geocode.py velocity_msk.h5 -l inputs/geometryRadar.h5 -x 0.00027778 -y -0.00027778 --bbox 32.0 32.5 130.1 130.5 | ||
|
|
||
| # To run this program do the following | ||
| los2slope.py -file ./geo_velocity.h5 -g ./geo_geometryRadar.h5 -o ./geo_velocity_slp.h5 -m 15 --slp_thresh 1.0 --scaling_factor_thresh 10 |
There was a problem hiding this comment.
Its a bit unclear how to use the thresholds for example here -m 15 --slp_thresh 1.0 --scaling_factor_thresh 10 are suggested but the defaults are "-m 25 --slp_thresh 5.0 --scaling_factor_thresh 10".
| cn = np.sin(np.deg2rad(aspect)) * np.cos(np.deg2rad(slope)) | ||
| cu = np.cos(np.deg2rad(slope)) | ||
|
|
||
| G = ve * ce + vu * cu + vn * cn |
There was a problem hiding this comment.
I think this could be tested using a simplified example that can be added to the tests directory as @yunjunz suggested. I made a rough example of this using asc_desc2horz_vert.py as reference. (example attached)
los2slope_test.py.txt
| return G | ||
|
|
||
|
|
||
| def los2slope(dlos, los_inc_angle, los_az_angle, slope, aspect, inps, step=20): |
There was a problem hiding this comment.
I second this suggestion
Description of proposed changes
Liine-of-Sight (LoS) displacement to the downslope component with reference to issue #1170 :
This code was developed with a particular focus on the landslide/slope deformation community. The downslope projection is performed by defining the satellite azimuth unit vector and the hillslope azimuth unit vector then projecting the LoS onto the slope using the dot product. This then leads to a “scaling factor” where the conversion to Displacement_Downslope = Displacement_LoS * Scaling Factor. The Python scripts will execute, but we request a full review of the equations to ensure the output is correct. One challenge we faced is that due to different conventions used, e.g., sometimes 0 is north, sometimes 0 is east, sometimes positive clockwise and sometimes positive anticlockwise, the literature is full of a mix of cosine and sine functions when performing these projections.
The code uses MintPy data as inputs to project the LoS deformation onto the downslope direction. The code used the existing code asc_desc2horz_vert.py as a starting point. The main new added functions are median_filter, smooth_dem, get_slp_asp, and get_design_matrix4slope.
The median_filter filters the existing digital elevation model into a less-noisy version, and its kernel size is controlled via input parameter “-m” which defaults to 25 pixels. This filter helps reduce large downslope scaling factors that result from noise in the DEM. The smooth_dem can be optionally used instead of median_filter in case of speckle-like noise where few pixels have extremely small or large values, which is not applied by default and can be used by using the “-si” parameter. The get_slp_asp function calculates and saves hillslope angle (slope) and hillslope aspect (aspect) as geotiff files using gdal (we used gdal instead of save_gdal util because the gdal function that calculates slope and aspect automatically saves the data as .tif, so it was easier to save on the fly than using save_gdal). The get_design_matrix4slope is where the downslope projection factor (“scaling factor”) is calculated based on Handwerger et al. (2019; https://www.nature.com/articles/s41598-018-38300-0#Sec4) to match the same convention; the aspect is converted to the anti-clockwise direction and heading angle is measured from the north. A full review of the equations and the codes would be greatly appreciated.
To test and run the function, please use the command; the test dataset provided is attached here:
SampleData.zip
To test the code, please run:
python los2slope.py -file ./T122_ASC/geo_velocity_msk.h5 -g ./T122_ASC/geo_geometryRadar.h5 -o ./T122_ASC/geo_velocityDownslope.h5 -m 15 --slp_thresh 5.0 --scaling_factor_thresh 5
For further review/discussion on this as well as review of the formula could @yunjunz, @taliboliver, @mgovorcin, @EJFielding and @alhandwerger provide some feedback?