Calculate intersection weights of line and grid

[maxmargraf]

This feature should implement the grid intersection code from pycomlink that calculates the weights from CML paths over a radar grid. It also should return the radar-along-CML time series calculated from the intersection weights.

There should be a function to get the grid intersection weights, one to calculate the radar_along_CML time series as well as simple options/functions to save or load intersection weights.

It could be discussed if similar code doing the same with point and grid e.g. from wradlib should be implemented along with this feature.

Todos:

• Check if not ascending grid coordinates lead to an error in calc_sparse_intersect_weights_for_several_cmls(). If so, there is the need for a function that checks and transposes the grid prior to calculating the intersection weights - see comment below
• get_grid_time_series_at_intersections() requires the dimension time in its pycmolink implementation but single time steps without time dimension should be supported

[cchwala]

It could be discussed if similar code doing the same with point and grid e.g. from wradlib should be implemented along with this feature.

I think this should come later. For sensor-merging code we would need on common approach (API) to get grid values at the sensors location. But I would implement that as a higher level function that calls the actual grid intersection code (or the code that gets grid values at or close to points).

[cchwala]

There was (maybe) a problem when I did a quick implementation (using the code from pycomlink) in the ragali_prototpy and used it in this notebook.

See this code snippet with my comment.

# TODO: Needs more thought...
# Either there is an error in my code or it has to do with the fact
# that I have to transpose the grids here when using grid_point_location='lower_left'
# because otherwise I get an error "y values must be ascending along axis 0".
# 
# For many CMLs the grid intersections are empty, even though the grid lon-lat 
# and the CML lon-lat seem to fit when viewed on a map...

intersect_weights = rgl.grid_intersection.calc_sparse_intersect_weights_for_several_cmls(
    ds_cmls=ds_cmls,
    x_grid=ds_radar.lon.values,
    y_grid=ds_radar.lat.values,
    #grid_point_location='lower_left',
    grid_point_location='center',
)

# Use grid intersection weights to calculate radar-CML path-averages
radar_along_cml = rgl.grid_intersection.get_grid_time_series_at_intersections(
    #grid_data=ds_radar.rainfall.transpose('time', 'x', 'y').values, # only needed if we have to transpose in the function above
    grid_data=ds_radar.rainfall.values,
    intersect_weights=intersect_weights,
)

When used with grid_point_location='lower_left' I had to (if I recall correctly) do ds_radar.lon.values.T or something like that. But I was not sure how that affected the results. Hence, I used center for the example notebook.

We should check again the logic in the existing code in pycomlink and make sure that the reason for the issue that I found is that the grid coordinate values are ascending or descending. If ascending vs descending is the problem we should raise an error or automatically transpose or flip. Since there are tests for the grid intersection code, we could also try to make the code more general. But that might be harder than just checking the input grid.

[cchwala]

The error that is raised when using grid_point_location='lower_left' is actually an error the results from checking that values are ascending:

    280 if (np.diff(grid[:, :, 1], axis=0) < 0).any():
--> 281     raise ValueError("y values must be ascending along axis 0")
    283 # Upper right
    284 ur_grid = np.zeros_like(grid)

ValueError: y values must be ascending along axis 0

It is not yet clear how to deal with descending coordinate values.