✅ resizing one layer to the size of another, and ambiguity in `rescale` function meaning

[gottacatchenall]

Description of the to-do item

As far as I can tell, there is not a function to take a raster and change its resolution to match the dimensions of a target raster. Doing this would be straightforward by iterating over the longitudes(a), latitudes(a) in the smaller raster and then access the values of the other raster at that coord, and have the getindex(b::SimpleSDMLayer, lat, long) handle that conversion.

When trying to do this I found rescale(layer1::T, layer2::T) where {T <: SimpleSDMLayer},which scales the values of a given raster to match a target, but not dimensions. To avoid this semantic ambiguity, perhaps the method rescale could rescale both values and dimensions by default with kw arguments to turn either off

[gottacatchenall]

my hacky workaround atm is

function make_same_dims(input::IT, target::TT) where {IT,TT <: SimpleSDMLayer}
	lat, long = collect(longitudes(target)), collect(latitudes(target))
	newgrid = similar(target)
	for lt in lat, lg in long
		newgrid[lt,lg] = input[lt,lg]
	end
	newgrid
end

[tpoisot]

That's a good point, but ideally something we might want to do with some sort of spatial interpolation. There's a way to get the haversine distance, do we can use this. I wouldn't necessarily call it rescale, maybe reproject? There's a little more thinking to do on this one before moving to a PR.

[mkborregaard]

FWIW, GeoData calls this "resample" and forwards the functionality to ArchGDAL.gdalwarp

[tpoisot]

That sounds like something we might do -- my current rough draft is

using SimpleSDMLayers
using StatsPlots
using Statistics

ref = SimpleSDMPredictor(WorldClim, BioClim, 3; left=-20.0, right=40.0, top=60.0, bottom=30.0)

target = SimpleSDMResponse(zeros(Float32, 200, 200), -12.056, 37.065, 33.1234, 55.06)

for p in keys(target)
    if isnothing(ref[p])
        target[p] = nothing
    else
        vals = SimpleSDMLayers._sliding_values(ref, p, 20.0)
        target[p] = isempty(vals) ? nothing : mean(vals)
    end
end

@gottacatchenall , given the info @mkborregaard provided, do you want to see if a version based on ArchGDAL.gdalwarp is possible? It's already a dependency, and we can write/read temp geotiff files...

[gottacatchenall]

@tpoisot Yeah GeoData.resample wraps it with arguments for interpolation rules, it isn't currently a dependency though

i've also written a few functions to take a set of SimpleSDMLayers and return PCA layers which i could include in the same PR

[tpoisot]

Two PRs! Mostly because I'd like to do some solid work on the PCA, to make sure we can drop layers into multivariate stats pretty much directly.

[tpoisot]

Also I wouldn't bring too much stuff from Geo until they figure out whether they use the Geometry types or not...

[gottacatchenall]

@tpoisot i'm not sure if GeoData.jl is affiliated with JuliaGeo, here's its deps, which don't seem too heavy

also in progress on opening prs for both PCA and this, from what i can tell interfacing with Mvstats is pretty easy

[tpoisot]

Re. GeoData that's still a lot of deps. I'm curious to see how far we can go using ArchGDAL only, since we depend on that no matter what.

[gottacatchenall]

As far as I can tell resample only uses ArchGDAL functions.I could port it, but it would effectively be copy/pasting which isn't ideal.

I suppose long-term the julian-esque would be to move all functions that rely on ArchGDAL to a shared lib and have any other heavy dependencies loaded from higher level packages

Thoughts? @tpoisot @mkborregaard @rafaqz

[mkborregaard]

So I don't think it would make sense to take a GeoData dep just to get resample. GeoData is a general raster package, very similar to R's terra (and their old raster package), and as such has a lot of overlapping functionality with SimpleSDMLayers (such as masking, cropping, extracting data, downloading data, aggregating, opening raster files etc etc). GeoData is more focused on being a general raster package, though, whereas in my understanding SimpleSDMLayers is focused on creating a smooth user experience for fitting SDM models.

One could consider having SimpleSDMLayers depend on GeoData, which would give all the raster operations here for free, and I guess make it easier to maintain, but that would change the focus of the package to be explicitly on creating the smooth SDM interface. That depends completely on @tpoisot 's vision for this package.

I think most of the added deps for geodata (e.g. DiskArrays, Flatten, Mmap, ) relate to the functionality for working with on-disk rather than in-memory rasters.

When I originally posted here my suggestion was that you could be inspired (or essentially copy) GeoData's resample function, which should also be doable.

[rafaqz]

I agree with that assessment @mkborregaard . Depending on GeoData.jl would be a larger change in philosophy and more like a decision to reduce the codebase here by half. GeoData generalizes many of the methods here, including resample.

Proper broadcasting, base array methods and the Tables.jl interface in GeoData.jl would also save work here. These things are hard to do properly.

So I wouldn't add a GeoData.jl dep just for resample. You can just copy code from GeoData.jl, and delete all of the generalization code. resample is just a wrapper to warp which is n-dimensionally generalised gdalwarp. Although @gottacatchenall I tend to agree that copy-pasting doesn't seem like a sustainable long-term strategy.

In terms of JuliaGeo, GeoData is no more affiliated than ArchGDAL is - its not a JuliaGeo package. Its my package, partly dog-fed by cesaraustralia but to do generic things better, not specifically to work with cesar packages.

It's also less tightly integrated with GeoInterface.jl than ArchGDAL is. I agree about the GeoInterface/GeometryBasics problem and will likely be involved in the solution in the coming months. But I don't think that's an insurmountable problem currently.

Many of the dependencies like DiskArrays/GeoInterface/ColorTypes/GeoFormatTypes are already ArchGDAL deps. HDF5 and NCDatasets are included because Requires.jl is not a great solution to handle dependency bounds, and we dont have proper glue packages yet. Please chip in there too if you would like the dependency problem resolved, I've used up my quota pushing that...