Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
wants to merge 2 commits into from

Conversation

asinghvi17
Copy link
Collaborator

@asinghvi17 asinghvi17 commented Oct 2, 2024

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

Copy link
Owner

@rafaqz rafaqz left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we not do this already with mapslices when skipmissing=false?

Or is it a performance issue

function _mapzonal(reducer, operator, x, ext::Extents.Extent; skipmissing = true)
cropped = crop(x; to = ext, touches = true)
prod(size(cropped)) > 0 || return missing
# We can't use skipmissing here, since it doesn't work on rasterstacks
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should just make this work, but it's just not totally clear what the rule is - any missing or all missing

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In here I've implemented the "all missing" variant but it's up for debate

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess the question is, does skipmissing mean "skip all values that are missing in general" or "skip all values that are missing because of my mask"? The first is the most useful because there can be functions that don't tolerate missings

@asinghvi17
Copy link
Collaborator Author

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
Copy link
Owner

rafaqz commented Oct 2, 2024

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
Copy link
Collaborator Author

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
Copy link
Owner

rafaqz commented Oct 3, 2024

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
Copy link
Collaborator Author

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

@rafaqz
Copy link
Owner

rafaqz commented Oct 3, 2024

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
Copy link
Collaborator Author

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
Copy link
Collaborator Author

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

@asinghvi17
Copy link
Collaborator Author

asinghvi17 commented Oct 3, 2024

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
Copy link
Owner

rafaqz commented Oct 3, 2024

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
Copy link
Owner

rafaqz commented Oct 3, 2024

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
Copy link
Collaborator Author

asinghvi17 commented Nov 13, 2024

Succeeded by #820

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants