Merge pull request #56 from EcoJulia/tp/rescale

Add a rescale function

Project.toml

@@ -1,7 +1,7 @@
 name = "SimpleSDMLayers"
 uuid = "2c645270-77db-11e9-22c3-0f302a89c64c"
 authors = ["Timothée Poisot <[email protected]>", "Gabriel Dansereau <[email protected]>"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"

docs/src/man/operations.md

@@ -7,4 +7,6 @@ latitudes
 longitudes
 clip
 mask
+rescale!
+rescale
 ```

docs/src/man/overloads.md

@@ -42,4 +42,5 @@ broadcast
 Statistics.mean
 Statistics.median
 Statistics.std
+Statistics.quantile
 ```

src/SimpleSDMLayers.jl

@@ -5,6 +5,7 @@ using HTTP
 using RecipesBase
 using ZipFile
 using Requires
+using Statistics
 
 include(joinpath("lib", "types.jl"))
 export SimpleSDMLayer, SimpleSDMResponse, SimpleSDMPredictor
@@ -33,7 +34,8 @@ export worldclim, bioclim, landcover
 include(joinpath("operations", "coarsen.jl"))
 include(joinpath("operations", "sliding.jl"))
 include(joinpath("operations", "mask.jl"))
-export coarsen, slidingwindow, mask
+include(joinpath("operations", "rescale.jl"))
+export coarsen, slidingwindow, mask, rescale!, rescale
 
 include(joinpath("recipes", "recipes.jl"))
 

src/datasets/ascii.jl

@@ -25,7 +25,6 @@ function ascii(file::AbstractString, datatype::Type{T}=Float64) where {T <: Numb
     for line_id in data_start:data_end
         M[:,nrows-(line_id-(data_start))] = parse.(datatype, split(lines[line_id]))
     end
-    println(permutedims(M))
     # Put data in the grid
     grid = convert(Matrix{Union{datatype,Nothing}}, permutedims(M))
     for i in eachindex(M)
@@ -67,6 +66,11 @@ function ascii(layer::SimpleSDMPredictor{T}, file::AbstractString; nodata::T=con
     return file
 end
 
+"""
+    ascii(layer::SimpleSDMResponse{T}, file::AbstractString; nodata::T=convert(T, -9999)) where {T <: Number}
+
+Writes a `layer` to a grid file, with a given `nodata` value. The layer must store numbers.
+"""
 function ascii(layer::SimpleSDMResponse{T}, file::AbstractString; nodata::T=convert(T, -9999)) where {T <: Number}
     return ascii(convert(SimpleSDMPredictor, layer), file; nodata=nodata)
 end

src/lib/overloads.jl

@@ -80,6 +80,7 @@ function Base.stride(layer::T; dims::Union{Nothing,Integer}=nothing) where {T <:
    dims == 1 && return lon_stride
    dims == 2 && return lat_stride
 end
+Base.stride(layer::T, i::Int) where {T<:SimpleSDMLayer} = stride(layer; dims=i)
 
 """
     Base.eachindex(layer::T) where {T <: SimpleSDMLayer}
@@ -299,7 +300,7 @@ Returns the non-`nothing` values of a layer.
 """
 function Base.collect(l::T) where {T <: SimpleSDMLayer}
     v = filter(!isnothing, l.grid)
-    return convert(Vector{eltype(l)}, v)    
+    return convert(Vector{typeof(v[1])}, v)    
 end
 
 """
@@ -334,4 +335,13 @@ function Base.hcat(l1::T, l2::T) where {T <: SimpleSDMLayer}
     new_grid = hcat(l1.grid, l2.grid)
     RT = T <: SimpleSDMPredictor ? SimpleSDMPredictor : SimpleSDMResponse
     return RT(new_grid, l1.left, l2.right, l1.top, l1.bottom)
+end
+
+"""
+    quantile(layer::T, p) where {T <: SimpleSDMLayer}
+
+Returns the quantiles of `layer` at `p`, using `Statistics.quantile`.
+"""
+function Statistics.quantile(layer::T, p) where {T <: SimpleSDMLayer}
+    return quantile(collect(layer), p)
 end

src/operations/rescale.jl

@@ -0,0 +1,71 @@
+_rescale(x, m, M) = (x .- minimum(x))./(maximum(x)-minimum(x)).*(M-m).+m
+
+"""
+    rescale!(layer::TI, template::TJ) where {TI <: SimpleSDMLayer, TJ <: SimpleSDMLayer}
+
+Changes the values of the layer given as its first argument, so that it has the
+same *range* as the values of the layer given as its second argument.
+Modification is done in-place.
+"""
+function rescale!(layer::TI, template::TJ) where {TI <: SimpleSDMLayer, TJ <: SimpleSDMLayer}
+    return rescale!(layer, extrema(template))
+end
+
+"""
+    rescale!(layer::TI, t::Tuple{T,T}) where {TI <: SimpleSDMLayer, T <: Number}
+
+Changes the values of the layer given as its first argument, so that it has the
+same *range* as the values given as a tuple of values. Modification is done
+in-place.
+"""
+function rescale!(layer::TI, t::Tuple{T,T}) where {TI <: SimpleSDMLayer, T <: Number}
+    occ = findall(!isnothing, layer.grid)
+    layer.grid[occ] .= _rescale(layer.grid[occ], t...)
+end
+
+"""
+    rescale(layer::TI, template::TJ) where {TI <: SimpleSDMLayer, TJ <: SimpleSDMLayer}
+
+Copying version of `rescale!`.
+"""
+function rescale(layer::TI, template::TJ) where {TI <: SimpleSDMLayer, TJ <: SimpleSDMLayer}
+    l = copy(layer)
+    return rescale!(l, extrema(template))
+end
+
+"""
+    rescale(layer::TI, t::Tuple{T,T}) where {TI <: SimpleSDMLayer, T <: Number}
+
+Copying version of `rescale!`.
+"""
+function rescale(layer::TI, t::Tuple{T,T}) where {TI <: SimpleSDMLayer, T <: Number}
+    l = copy(layer)
+    return rescale!(l, t)
+end
+
+
+"""
+    rescale!(layer::T, p::Vector{Real}) where {T <: SimpleSDMLayer}
+
+Rescale the values of a `layer` so that they match with the quantiles given in
+`p`. Internally, this uses the `Statistics.quantile` function.
+"""
+function rescale!(layer::T, p::Vector{TI}) where {T <: SimpleSDMLayer, TI <: AbstractFloat}
+    q = reverse!(quantile(layer, p))
+    occupied = findall(!isnothing, layer.grid)
+    v = collect(layer)
+    for i in 1:length(q)
+        layer.grid[occupied[findall(x -> x <= q[i], v)]] .= reverse(p)[i]
+    end
+    return layer
+end 
+
+"""
+    rescale(layer::T, p::Vector{Real}) where {T <: SimpleSDMLayer}
+
+Copying version of `rescale!`.
+"""
+function rescale(layer::T, p::Vector{TI}) where {T <: SimpleSDMLayer, TI <: AbstractFloat}
+    l = copy(layer)
+    return rescale!(l, p)
+end

test/overloads.jl

@@ -46,7 +46,7 @@ Y[0.1,0.1] = 4.0
 
 Z = convert(SimpleSDMPredictor, Y)
 V = collect(Z)
-@test typeof(V) == Vector{eltype(Z)}
+@test typeof(V) == Vector{Float64}
 
 # hcat / vcat
 l1 = worldclim(1, left=0.0, right=10.0, bottom=0.0, top=10.0)

test/rescale.jl

@@ -0,0 +1,17 @@
+module SSLTestRescale
+using SimpleSDMLayers
+using Test
+
+M1 = rand(Float64, (5, 10)).*4
+S1 = SimpleSDMPredictor(M1, 0.0, 1.0, 0.0, 1.0)
+
+M2 = (rand(Float64, (5, 10)).+1).*5
+S2 = SimpleSDMPredictor(M2, 0.0, 1.0, 0.0, 1.0)
+
+@test extrema(rescale(S1, S2)) == extrema(S2)
+@test extrema(rescale(S2, S1)) == extrema(S1)
+
+rescale!(S1, [0.0, 0.5, 1.0])
+@test sort(unique(collect(S1))) == [0.0, 0.5, 1.0]
+
+end

test/runtests.jl

@@ -9,6 +9,7 @@ tests = [
    "basics" => "basics.jl",
    "overloads" => "overloads.jl",
    "iteration" => "iteration.jl",
+   "rescale" => "rescale.jl",
    "generated" => "generated.jl",
    "import" => "dataread.jl",
    "worldclim" => "worldclim.jl",