read_raster_sample: fix antimeridian and implement raster gradients #551
Conversation
| suffix = Path(path).suffix | ||
| if suffix in supported_suffixes: | ||
| path = f"/vsi{supported_suffixes[suffix]}/{path}" | ||
| return str(path) |
There was a problem hiding this comment.
Previously, the prefixing was manually applied for gz-files in three places. But we can easily add support for other archive formats and we can avoid redundant code by introducing this helper function.
Note that the prefix "looks like a path". Since the result with the added prefix is not really a path, I would not convert it to a Path object. And for consistency, I would also convert to a str in the cases where no prefix is added at all.
| @@ -942,14 +942,12 @@ def test_dist_to_coast_nasa(self): | |||
| [1.96475615, 45.23249055], | |||
| ]) | |||
| dists = [-3000, -1393549.5, 48.77] | |||
| dists_lowres = [416.66666667, 1393448.09801077, 1191.38205367] | |||
| # Warning: This will download more than 300 MB of data! | |||
| dists_lowres = [729.1666667, 1393670.6973145, 945.73129294] | |||
There was a problem hiding this comment.
Here is why I had to change the test data:
The definition of the intermediate raster in read_raster_sample changed. Previously, the intermediate raster depended on the interpolation points, which can lead to inconsistencies when the function is called several times with the same intermediate resolution but different interpolation points. Now, the intermediate raster is more consistent since it is always relative to the original raster's origin.
This only affects calls to read_raster_sample where an intermediate_res is specified, i.e. where the user doesn't want to sample from the actual raster data, but from a low resolution version of it. I don't think that we have to be super-consistent here because it's a rather unreliable approach anyway. I assume that applications that require a high level of consistency wouldn't use this feature, but sample directly from the original raster.
climada/util/coordinates.py
Outdated
| gradient : np.array of shape (npoints, 2), optional | ||
| If grad=True, the raster gradient at each of the given coordinate points is returned. | ||
| The first/second value in each row is the derivative in lat/lon direction (lat is first!). |
There was a problem hiding this comment.
I prefer methods returning unambiguous data types, or if need be an Option thereof. So that after y = f(x) y is always of type T (or None). Here y could be an array or a pair of arrays.
In this particular case I would have preferred two methods sharing the job, something like this:
values = read_raster_sample_without_gradient(...)
gradient = read_raster_sample_gradient(values, ...)Since I assume that would be a waste of cpu in some way, I'd suggest to always return a pair
def read_raster_sample(...) -> Tuple[np.array, Option[np.array]]:
...There was a problem hiding this comment.
Good point, let me think about this again, I will come back to you with a solution ...
There was a problem hiding this comment.
I solved this problem by leaving the old functionality of read_raster_sample untouched, and adding a new function read_raster_sample_with_gradients. There is almost no redundant code, because I moved some of the code to a helper function that is used by both functions.
If you like that fix, feel free to merge.
There was a problem hiding this comment.
👍 much better than my suggestion as it keeps backward compatibility!
|
Looks good to me. 😃 Although I went over it only superficially. |
|
🙌 |
…551) * u_coord: _add_gdal_vsi_prefix * read_raster_sample: gradient feature * Fix tests for dist_to_coast_nasa * u_coord._raster_gradient: lat direction first * u_coord.read_raster_sample: support different interp methods for data and its gradient * read_raster_sample_with_gradients * u_coord: pylint unused variable Co-authored-by: emanuel-schmid <[email protected]>
Changes proposed in this PR:
u_coord.read_raster_sampleby making use of the already existing and antimeridian-proof functionu_coord.read_raster_bounds.u_coord.read_raster_sample(calledu_coord.read_raster_sample_with_gradients) that also returns the gradient of the raster at the given points.The raster gradient feature will be used in the upcoming new TC rainfall model, since the rainfall depends on the slope of the terrain.
PR Author Checklist
develop)PR Reviewer Checklist