[WIP] `mapzonal`, allowing you to apply some operation per raster cell, then `zonal` it

[asinghvi17]

This PR introduces a function mapzonal that looks like this:

mapzonal(reducer, operator, raster; of = geometries)

operator is fed each slice of the raster in all dimensions except X and Y. For a RasterStack that has only X and Y dimensions, for example, operator would be passed the output of rasterstack[X(ind), Y(ind)] for each valid index of X and Y.

If you wanted to calculate NDVI from some super-high-fidelity rasterstack, for example, you would perform the computation on demand in operator.

Then, reducer is passed the vector generated by operator.(each_cell_of_rasterstack). This can be e.g. sum, or even a histogram function.

Motivation

The motivation here is that zonal iterates over the layers of a rasterstack before passing them to the reducer - this can be good, and necessary for the most common workflows, but we needed to operate on the raw values of the individual layers simultaneously.

Some interesting usecases might be:

1. Imagine I have a Raster of soil type and canopy height and I want to know the average canopy height for each soil type within each watershed
2. I have a map of slope-aspect and a map of temperature and I want to know the maximum temperature for each aspect for different mountain ranges
3. I want to calculate the average surface lapse rate (temperature vs elevation) over glacier surfaces given a DEM and raster of surface temperature

TODOs

• [ ] Tests
• [ ] Comprehensive documentation
• [ ] Performance on regular Rasters with >2 dims (it's currently 10x slower, since indexing into a raster with a dimindex of lacking rank causes a new raster to be allocated)
• [ ] Benchmarking

[asinghvi17]

To do this with mapslices you would have to essentially remake zonal anyway (crop + mask and walk down the tree of geometry), but I can check the performance of that

[rafaqz]

I mean inside the zonal function, can you just call mapslices on the RasterStack it gives you?

(The idea is good, I'm just wary of feature creep)

[asinghvi17]

It doesn't give you a rasterstack though, since it iterates over layers...

Unless you mean changing the implementation? That would be breaking I think

[rafaqz]

Hmm yes that's annoying. How do we get the best of both behaviours...

But can we workshop how to get this all into zonal elegantly? I'm happy to break things at this stage, or at least add options rather than more functions

[asinghvi17]

Hmm, it could be a switch that does this internally?

[rafaqz]

bylayer=true ?

False gives you the whole stack, or iterates over NamedTuples when skipmissing=true

But we need to fix skipmissing on RasterStack for that to work

(Guess this might be a useful keyword in many places a function is used on stack values)

[asinghvi17]

There are two scenarios here - "set of rasters" (where you return multiple columns) or "linked rasters" (where you want a slice at those x and y coordinates)

[asinghvi17]

The thing is that methods won't necessarily work the same in both, or work at all...

[asinghvi17]

Hmm, with context bylayer does make sense. How would you treat 3d rasters in zonal then? Or a stack with a 3d timeseries raster and a 2d DEM or something

[rafaqz]

But if we have bylayer=true, skipmissing=false and get passed a complete masked RasterStack in f, will that let you do what you need with mapslices?

[rafaqz]

How  would you treat 3d rasters in zonal then?

With skipmissing=false you already get the whole object, whatever the dimensions are, zonal just crops/masks over X/Y. Otherwise it's a skipmissing iterator over all dimensions so mean will just work.

With bylayer=false stacks can be the same, and you handle slicing however you need to as with any DimStack

[asinghvi17]

Succeeded by #820