Conversation
b57caba to
fb36168
Compare
fb36168 to
d42f267
Compare
|
I think this works with the tests that you added here because the affine transformation is linear (symmetrical) and also because how numpy works with negative indices. I don't think this is will always work without something to handle wraparound, though. |
| indexed = ds.pipe(assign_index) | ||
| # We lose existing attrs when calling ``assign_index``. | ||
| assert_equal(ds.sel(lon=[220, 240]), indexed.sel(lon=[-140, -120])) | ||
| assert_equal(ds.sel(lon=220), indexed.sel(lon=-140)) |
There was a problem hiding this comment.
For example
assert_equal(ds.sel(lon=[0, 0]), indexed.sel(lon=[0, 359.9]))should pass with some proper wraparound but raises
IndexError: index 180 is out of bounds for axis 2 with size 180
for indexed.sel(lon=359.9)
There was a problem hiding this comment.
hmmm.. using method="nearest" by default is quite wrong, but CoordinateTransformIndex requires it
|
I would be nice to also have an option for shifting the range [0, 360] from/to [-180, 180]. |
d42f267 to
f4cd46a
Compare
Initially I thought that is what this PR was addressing. I guess I didn't realize that the affine "automatically" handles out-of-range values (for certain cases if I understand @benbovy's comment above correctly)! I was expecting a function or keyword argument to
I guess CRS is needed to check that coordinates are lon/lat and raise if the CRS is projected (where there is no wrapping)? https://pyproj4.github.io/pyproj/stable/api/crs/crs.html#pyproj.crs.CRS.is_projected In my mind there are two scenarios: 1. Selection using an affine that generates coordinates in [0,360] & user supplies out-of-bounds longitude that needs to be wrapped into the [0,360] range def wrap_coord(lon):
# same as converting from [-180,180] to [0,360] range
return lon % 360da.sel(lon=-120) == ds.isel(lon=240)
da.sel(lon=365) == ds.isel(lon=5)2. Selection on affine that generates coordinates in [-180,180] & user supplies longitude in [0, 360] range that needs to be converted da.sel(lon=240) == da.sel(lon=-120)def convert_to_negative_west(lon):
return ((lon - 180) % 360) - 180 |
|
Naive question: is it common to have unprojected lat/lon raster datasets with a rectilinear / no-rotation affine transform? It is pretty common for model outputs but they don't use affine transforms. |
|
FWIW Apache SIS has a special WraparoundTransform class. I guess one option could be to provide for convenience a similar class that would interoperate with |
|
Yeah, this is a lot more subtle than I initially realized, I won't be able to do it before SciPy. I am going to issue an alpha release now. |

Closes #26
Apparently we inherit this from
affine, so no CRS needed. Does anyone understand this?