Feature parity with R raster

[rafaqz]

Essentially:

• polygon operations like masking, union, intersect etc
• a few helper functions like extend - functions that make the spatial extent larger are missing
• summation helpers like area and distance

[floswald]

was just working with the follow-up of raster in R, stars and thought if only one could do that stuff in julia. Well! Here we are! Great package!
I thought it's useful to add couple more useful features here. I think most will be about documenting things in a certain kind of way such that people with different backgrounds can access the functionality. I could help with that.

• aggregate by a Geom. If I understand correctly the current Rasters.jl method aggregate is kind of low level in that the user is expected to supply the indices where to aggregate over with method. The R solution here is nice because I can have a raster prec and a sf object nc say and do

a = aggregate(prec, by = nc, FUN = max)

as long as they are in the same CRS. Again, I think that might be achievable with the current mask approach etc, but maybe that's worth an extension.

• you are tracking intersection and union etc any ways here Intersect between SpatialLines and AbstractRaster #233

[rafaqz]

Thanks for commenting here! This is the best kind of feedback :)

Yes aggregate could be better. It was written before the other methods, and isn't as slick as say mask.

But if you could help me out a little more with this function, I dont understand how the aggregation amount comes from nc:

a = aggregate(prec, by = nc, FUN = max)

Does nc cover the area that will be aggregated to? So the X/Y aggregation pixels are those that would just cover the extend of nc?

[floswald]

That's exactly it. It can be any polygon or other shape, I guess even a line. It's related to another issue I saw you had I think, where someone wants to intersect with a line. So yes it's just about not having to think about creating a mask and aggregating over that.

[floswald]

For example I have a use case where the Raster represents some measures of soil quality in a country. Then I have 100 circles which are made as buffers around certain points. I want to know soil quality in those 100 circles. nc would just be a dataframe with a sf::geometry column (each row a different polygon representing the relevant circle).

[rafaqz]

I'm still not sure what the return value is though. Do you mean it's a raster with lower resolution? How does that resolution exactly match your 100 circles?

It sounds like you really want zonal stats? If I wanted to get mean values for 100 circles I would crop then mask them individually and then take the mean of the the skipmissing values of each. That's a super fast operation in Raster.jl. I've meant to add a zonal function to do that, it's only 1 line but would make it easier.

(maybe if you can link to the docs - is this a stars function?)

[floswald]

The return value is a vector as long as geometries in the by argument. So 100 numbers in this example.
You are correct with your approach, yes.
Here the stars docs https://r-spatial.github.io/stars/reference/aggregate.stars.html

[rafaqz]

Ok perfect. Thats what the new zonal PR does.

aggregate here returns a raster. Im not sure if it should share a name with zonal like in stars or remain separate? aggregate can do things like take the center value rather than apply a function, which diesnt really make sense for polygons.