Add Rasters.sample

[tiemvanderdeure]

This PR adds a Rasters.sample function in a StatsBase extension, as suggested here: #771

The keywords correspond to extract keywords. Some of the internals are also re-used.

I would suggest not exporting this to avoid confusion with StatsBase.sample, which has a totally different API.

If skipmissing = true, Rasters.sample and extract will always return the same types - otherwise it depends on the element type of the Raster. I tried to avoid creating Union{Missing, ...} as much as possible.

Closes #771

[tiemvanderdeure]

Rasters.sample is about just as fast as extract.

What slows it down in some cases is the conversion that happens if skipmissing = true to get rid of a Union type.

Just something like

A = [missing 7.0f0; 2.0f0 missing]
ga = Raster(A, (X(1.0:1:2.0), Y(1.0:1:2.0)); missingval=missing) 
Rasters.sample(ga, 1_000_000, skipmissing = true)

takes 0.3 seconds because of this:

[tiemvanderdeure]

Tests fail because of [5453] signal 11 (128): Segmentation fault. What a joy.

[rafaqz]

What slows it down in some cases is the conversion that happens if skipmissing = true to get rid of a Union type.

Do we know the answer from the type? What do you think is stopping it propagating?

It could be your g::Type instead of Type{G} I commented on in the other PR. Lets workshop it over there rather than splitting the disscussion. We probably need a benchmark in #790 to show this problem

(The segmentation fault is of course GDAL)

[tiemvanderdeure]

Sorry for all the mess with the rebases. I added an extra keyword to _rowtype to distinguish between whether index/geometry fields may be missing, and wether values extracted from the raster may be missing. The former is never the case in sample. So now we can completely reuse the same internal.

[tiemvanderdeure]

Somehow NamedTuple(x[idx])[K] is not type stable and I can't quite figure out why...

[tiemvanderdeure]

Was going to say that looks dubious. Probably we should subset the stack with st[K] before we get to this method to reduce indexing. Then you can just do:

x[idx]

I tried this but it failed tests, because st[K] does not necessarily have all the dimensions of st...

[tiemvanderdeure]

This should be type stable now in all cases.

I would still like to change the geometry keywords, which now is not very flexible. It could be nice to allow users to specify the type of the geometry column, and which dimensions should be used. So something like geometry = (Y=Y, X=X) would return a Y-X NamedTuple with data from the Y and X dimensions.

[tiemvanderdeure]

Tests are still failing on a segmentation fault in resample 😭

[rafaqz]

Ahh that's interesting about the dimensions. I guess we need to be consistent and return a table of the same length every time.

You can put the resample tests last to get around them for now. But we will need to deal with them before bumping a version.

[rafaqz]

We can totally let users specify the geometry like that too.

[tiemvanderdeure]

Ahh that's interesting about the dimensions. I guess we need to be consistent and return a table of the same length every time.

It's how we designed extract so I think it's good to be consistent. extract(stack, point; names) is not the same as extract(stack[names], points).

[rafaqz]

Maybe we should put that exact line in the extract docs

[tiemvanderdeure]

Explicitly using the eltype of the stack like this was the only way I could find in which this would compile away. But it looks somehow ugly and I don't know if there's a better way?

[rafaqz]

Like what? Is this comment on the wrong line?

[tiemvanderdeure]

Maybe we should put that exact line in the extract docs

Maybe. This feature justifies the having the names keyword in the first place, but it probably has a limited number of use cases.

[tiemvanderdeure]

I added the option to specify a geometry type, even though it now only works for Tuple and NamedTuple. For other types we can't just call convert, so it's not so easy to make it more generic.

[rafaqz]

It's weird this NT trick is needed, maybe DD stack getindex stability is broken on 1.11?

Looks like something needs to be forced somewhere

[tiemvanderdeure]

I know, nothing else worked. Even AbstractRasterStack{K} and then NamedTuple{K} wasn't type stable.

[rafaqz]

No that makes sense, it's the eltype that's not propagating.

We might need some more Base.assume_effects in DD indexing

[rafaqz]

Is there a test for the NamedTuple geometry?

[tiemvanderdeure]

Is there a test for the NamedTuple geometry?

Yes, although it's pretty rudimentary.

    @test typeof(first(
        Rasters.sample(StableRNG(123), test, 2, geometry = (X = X, Y = Y))
    ).geometry) <: NamedTuple{(:X, :Y)}

[rafaqz]

Ok that's fine, I like the syntax of just passing the tuple. Does it work for symbols and constructed dims too?

Probably want to convert those values to constructed basedims for type stability, Types get lost in the tuple

[tiemvanderdeure]

It just calls commondims, so I think constructed dims yes, symbols no. This is parsed and a _True() or _False(), a dimtuple, and a type are passed to _sample, so that should all be type stable.

NamedTuple types still allocate though, idk if/how to avoid that?

[rafaqz]

That's what I mean. For NamedTuple you need to pass through a NamedTuple of constructed dimensions, not types or symbols. That can be type stable making NamedTuple points, the others can't be

DD always tries to convert to constructed dims asap for this reason. If it can inline or constprop to the point of construction it can be totally stable using types. Otherwise it's super bad

[tiemvanderdeure]

I just added a dispatch so Rasters.sample(myraster) also works and returns just a NamedTuple. I don't have anything else to add here. The methods tests are before resample so if it segfaults again then tests have passed.

[rafaqz]

Amazing let's merge

[rafaqz]

Ohh maybe some links in the methods page of the docs would have helped

[tiemvanderdeure]

True, I forgot about that