Closed alessandrofelder closed 2 years ago
Looks like there's a (numerical?) problem with the domain we are passing to the nearneighbour interpolation in some cases?
nearneighbor [INFORMATION]: Cartesian input grid
nearneighbor [INFORMATION]: Given domain implies x_inc = 0.00100044
nearneighbor [INFORMATION]: Given domain implies y_inc = 0.00100073
nearneighbor [INFORMATION]: Cartesian input grid
nearneighbor [INFORMATION]: Grid dimensions are n_columns = 758, n_rows = **916** # !!!should be 917!!!
nearneighbor [INFORMATION]: Number of sectors = 4, minimum number of filled sectors = 2
nearneighbor [INFORMATION]: Processing input table data
Values of parameters passed to nearneighbor
call:
diff_grid = nearneighbour(
diff,
region=base_grid.region,
spacing=base_grid.spacing,
S=2 * max_spacing,
N=4,
E=NODATA_VAL,
verbose=True,
)
Looks like I can hack it by slightly reducing the spacing passed to nearneighbor
but this is not a viable way forward in itself. I also still don't understand why it happens exactly.
diff_grid = nearneighbour(
diff,
region=base_grid.region,
spacing=0.9995*base_grid.spacing,
S=2 * max_spacing,
N=4,
E=NODATA_VAL,
verbose=True,
)
Full diagnosis suggests underlying mismatch due to how gmt grdsample
by default slightly adapts the increments (passed by -I
) of the resampled grid to fit the domain. We could use "-I<spacing>+e"
as a flag to slightly adapt the domain instead and guarantee the spacing.
Somewhere in L124-148 of remove_restore.py it's possible (I've only managed to reproduce for a larger dataset) for the grid dimensions to get mismatched. With some of @dimitrasal 's data, the base_grid is 917x758, but the after the
nearneighbour
interpolation, the diff_grid becomes 916x758. Haven't found time to go deeper yet, will keep investigating next week, but thoughts/help appreciated.