Make RasterIndex CRS-aware#31
Conversation
Co-authored-by: Deepak Cherian <dcherian@users.noreply.github.com>
|
I integrated @scottyhq this PR doesn't fully supersede your PR #29, which adds some valuable instructions in README as well as a nice example of integration with |
|
Seems nice to me. Can we add some magic to |
I probably won't get to it next week, so if you'd like feel free to take it over! Otherwise leave it open and I'll adjust it later |
| ) | ||
|
|
||
| def join(self, other: RasterIndex, how: JoinOptions = "inner") -> RasterIndex: | ||
| if not self._proj_crs_equals(cast(CRSAwareIndex, other), allow_none=True): |
There was a problem hiding this comment.
confused here, shouldn't alignment check the spatial_ref index and error there too. Why does RasterIndex need the check? Is it because the order in which alignment is checked is unpredictable?
There was a problem hiding this comment.
It is needed for the case where RasterIndex has a CRS defined but the dataset / dataarray has no spatial_ref coordinate with xproj.CRSIndex.
There was a problem hiding this comment.
do we need to support that case. I was assuming we could just defer to xproj for all that.
There was a problem hiding this comment.
Unfortunately we cannot easily detect within the raster index here whether or not there is a related CRSIndex that will be used in the alignment. There's also the possibility that the internal CRS of the raster index is not in-sync with the CRSIndex. By design, xproj.CRSIndex and CRS-aware indexes (like RasterIndex) are not tightly coupled (I haven't found a way to do it nicely).
I guess the big question is: do we need to add a CRS property to RasterIndex at all? We add it here because at some point it will be useful for .sel, etc.?
There was a problem hiding this comment.
I guess the big question is: do we need to add a CRS property to RasterIndex at all? We add it here because at some point it will be useful for .sel, etc.?
Yeah I agree. My inclination is that we are going to bundle in CRS of some sort, then we may as wrap the Xproj index. So perhaps RasterIndex should stay ignorant of CRS and handle only the coordinate transform pieces.
We can handle the CRS based sel, either as an accessor function that knows about both XProjIndex & RasterIndex; or a new Index that composes the two. SInce they are orthogonal in concept, I'm not sure there's much value to a composed index.
What does XVec use the CRS for? Remapping the geometries passed in .xvec.query?
There was a problem hiding this comment.
We don't use CRS much in Xvec. I have implemented as we needed to store it but it kind of hangs around for now. There's certainly scope to make more use of it though.
There was a problem hiding this comment.
There's certainly scope to make more use of it though.
Just curious, how would you use it?
There was a problem hiding this comment.
I have the feeling that we could lost ourselves again in long discussions without making real progress, though.
Yeah I agree. I responded here . But I think we should move forward with this (#22 and #26 are good-enough reasons for me). We should document these choices, and release something that can be experimented with.
There was a problem hiding this comment.
Just curious, how would you use it?
Depends how much work do we want to do for a user under the hood. If we follow geopandas model, we should use it to check if the CRS of two objects that are being used together (e.g. zonal stats, possibly indexing) match and raise otherwise. If we want to do one more step, we can automatically re-project when they do not align, do the zonal stats using the transformed geom and use the original geom as the resulting index.
|
I just copied the notebook example from #29 here. Now we hit xarray-contrib/xproj#33 but it can be easily fixed. At least trying to align matching transforms but non-matching CRSs still raises an error.
Although this PR adds a CRS property to RasterIndex, it is opt-in so by default RasterIndex is ignorant of CRS. I suggest to leave it as-is (so that RasterIndex is nicely integrated with the latest version of Xproj) but not advertise too much that RasterIndex has its own CRS property (i.e., in the example notebook we just set an This should be ready then. |
|
Let's merge and iterate. |
4 participants
As a preliminary to #29.
I mostly copied the CRS-handling logic from
xvec.GeometryIndex(although instead of doing this, I'm thinking of factoring it out intoxproj.ProjIndexMixinand suggest that bothxvec.GeometryIndexandrasterix.RasterIndexexplicitly inherit from that mixin class so we can maintain that CRS-handling logic in one common place).Also added a few tests (RasterIndex still needs more tests, in a follow-up PR).