diff --git a/Project.toml b/Project.toml index a1a35b08..134198cc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SimpleSDMLayers" uuid = "2c645270-77db-11e9-22c3-0f302a89c64c" authors = ["Timothée Poisot "] -version = "0.2.1" +version = "0.2.2" [deps] GDAL = "add2ef01-049f-52c4-9ee2-e494f65e021a" diff --git a/src/SimpleSDMLayers.jl b/src/SimpleSDMLayers.jl index 81d826ba..8cfa826c 100644 --- a/src/SimpleSDMLayers.jl +++ b/src/SimpleSDMLayers.jl @@ -28,7 +28,7 @@ export coarsen include(joinpath("recipes", "recipes.jl")) # HACK but this fixes the export of clip when GBIF or others are loaded -clip(p::T) where {T <: SimpleSDMLayer} = nothing +clip(::T) where {T <: SimpleSDMLayer} = nothing export clip function __init__() diff --git a/src/bioclimaticdata/worldclim.jl b/src/bioclimaticdata/worldclim.jl index 7d3cf8a4..0bf09978 100644 --- a/src/bioclimaticdata/worldclim.jl +++ b/src/bioclimaticdata/worldclim.jl @@ -1,5 +1,5 @@ """ - worldclim(layers::Vector{Int64}; resolution::AbstractString="10", path::AbstractString="assets") + worldclim(layers::Vector{T}; resolution::AbstractString="10", path::AbstractString="assets") where {T <: Integer} Download and prepare WorldClim 2.0 bioclimatic variables, and returns them as an array of `SimpleSDMPredictor`s. Layers are called by their number, from 1 to 19. @@ -41,10 +41,11 @@ will be much faster. | 19 | Precipitation of Coldest Quarter | """ -function worldclim(layers::Vector{Int64}; resolution::AbstractString="10", path::AbstractString="assets") - @assert all(1 .≤ layers .≤ 19) +function worldclim(layers::Vector{T}; resolution::AbstractString="10", path::AbstractString="assets") where {T <: Integer} + all(1 .≤ layers .≤ 19) || throw(ArgumentError("The number of the layers must all be between 1 and 19")) isdir(path) || mkdir(path) - @assert resolution ∈ ["2.5", "5", "10"] + resolution ∈ ["2.5", "5", "10"] || throw(ArgumentError("The resolution argument ($(resolution) must be 2.5, 5, or 10")) + codes = [lpad(code, 2, "0") for code in layers] paths = [joinpath(path, "wc2.0_bio_$(resolution)m_$(code).tif") for code in codes] @@ -85,15 +86,15 @@ function worldclim(layers::Vector{Int64}; resolution::AbstractString="10", path: end """ - worldclim(layer::Int64; x...) + worldclim(layer::T; x...) where {T <: Integer} Return a single layer from WorldClim 2.0. """ -worldclim(layer::Int64; x...) = worldclim([layer]; x...)[1] +worldclim(layer::T; x...) where {T <: Integer} = first(worldclim([layer]; x...)) """ - worldclim(layers::UnitRange{Int64}; x...) + worldclim(layers::UnitRange{T}; x...) where {T <: Integer} Return a range of layers from WorldClim 2.0. """ -worldclim(layers::UnitRange{Int64}; x...) = worldclim(collect(layers); x...) +worldclim(layers::UnitRange{T}; x...) where {T <: Integer} = worldclim(collect(layers); x...) diff --git a/src/lib/basics.jl b/src/lib/basics.jl index d0878682..3b3ea52b 100644 --- a/src/lib/basics.jl +++ b/src/lib/basics.jl @@ -1,29 +1,37 @@ """ - latitudes(p::T) where {T <: SimpleSDMLayer} + latitudes(layer::T) where {T <: SimpleSDMLayer} Returns an iterator with the latitudes of the SDM layer passed as its argument. +This returns the latitude at the center of each cell in the grid. """ -function latitudes(p::T) where {T <: SimpleSDMLayer} - grid_size = stride(p; dims=2) - centers = range(p.bottom+grid_size; stop=p.top-grid_size, length=size(p, 1)) +function latitudes(layer::T) where {T <: SimpleSDMLayer} + grid_size = stride(layer; dims=2) + centers = range(layer.bottom+grid_size; stop=layer.top-grid_size, length=size(layer, 1)) return centers end """ - longitudes(p::T) where {T <: SimpleSDMLayer} + longitudes(layer::T) where {T <: SimpleSDMLayer} Returns an iterator with the longitudes of the SDM layer passed as its argument. +This returns the longitudes at the center of each cell in the grid. """ -function longitudes(p::T) where {T <: SimpleSDMLayer} - grid_size = stride(p; dims=1) - centers = range(p.left+grid_size; stop=p.right-grid_size, length=size(p, 2)) +function longitudes(layer::T) where {T <: SimpleSDMLayer} + grid_size = stride(layer; dims=1) + centers = range(layer.left+grid_size; stop=layer.right-grid_size, length=size(layer, 2)) return centers end -function are_compatible(l1::FT, l2::ST) where {FT <: SimpleSDMLayer, ST <: SimpleSDMLayer} - @assert size(l1) == size(l2) - @assert l1.top == l2.top - @assert l1.left == l2.left - @assert l1.bottom == l2.bottom - @assert l1.right == l2.right +""" + _layers_are_compatible(l1::X, l2::Y) where {X <: SimpleSDMLayer, Y <: SimpleSDMLayer} + +This returns + +""" +function _layers_are_compatible(l1::X, l2::Y) where {X <: SimpleSDMLayer, Y <: SimpleSDMLayer} + size(l1) == size(l2) || throw(ArgumentError("The layers have different sizes")) + l1.top == l2.top || throw(ArgumentError("The layers have different top coordinates")) + l1.left == l2.left || throw(ArgumentError("The layers have different left coordinates")) + l1.bottom == l2.bottom || throw(ArgumentError("The layers have different bottom coordinates")) + l1.right == l2.right || throw(ArgumentError("The layers have different right coordinates")) end diff --git a/src/lib/geotiff.jl b/src/lib/geotiff.jl index 938c9cc7..99958878 100644 --- a/src/lib/geotiff.jl +++ b/src/lib/geotiff.jl @@ -1,3 +1,9 @@ +""" + geotiff(tiff_file; T::Type=Float64) + +The geotiff function reads a geotiff file, and returns it as a matrix of type +`T`. +""" function geotiff(tiff_file; T::Type=Float64) # Register GDAL drivers GDAL.gdalallregister() diff --git a/src/lib/overloads.jl b/src/lib/overloads.jl index 6fdaeca0..fda0c049 100644 --- a/src/lib/overloads.jl +++ b/src/lib/overloads.jl @@ -9,92 +9,92 @@ import Base: eltype import Base: convert """ - Base.convert(::Type{SimpleSDMResponse}, p::T) where {T <: SimpleSDMPredictor} + Base.convert(::Type{SimpleSDMResponse}, layer::T) where {T <: SimpleSDMPredictor} Returns a response with the same grid and bounding box as the predictor. """ -function Base.convert(::Type{SimpleSDMResponse}, p::T) where {T <: SimpleSDMPredictor} - return copy(SimpleSDMResponse(p.grid, p.left, p.right, p.bottom, p.top)) +function Base.convert(::Type{SimpleSDMResponse}, layer::T) where {T <: SimpleSDMPredictor} + return copy(SimpleSDMResponse(layer.grid, layer.left, layer.right, layer.bottom, layer.top)) end """ - Base.convert(::Type{SimpleSDMPredictor}, p::T) where {T <: SimpleSDMResponse} + Base.convert(::Type{SimpleSDMPredictor}, layer::T) where {T <: SimpleSDMResponse} Returns a predictor with the same grid and bounding box as the response. """ -function Base.convert(::Type{SimpleSDMPredictor}, p::T) where {T <: SimpleSDMResponse} - return copy(SimpleSDMPredictor(p.grid, p.left, p.right, p.bottom, p.top)) +function Base.convert(::Type{SimpleSDMPredictor}, layer::T) where {T <: SimpleSDMResponse} + return copy(SimpleSDMPredictor(layer.grid, layer.left, layer.right, layer.bottom, layer.top)) end """ - Base.convert(::Type{Matrix}, p::T) where {T <: SimpleSDMLayer} + Base.convert(::Type{Matrix}, layer::T) where {T <: SimpleSDMLayer} Returns the grid as an array. """ -function Base.convert(::Type{Matrix}, p::T) where {T <: SimpleSDMLayer} - return copy(p.grid) +function Base.convert(::Type{Matrix}, layer::T) where {T <: SimpleSDMLayer} + return copy(layer.grid) end """ - Base.eltype(p::T) where {T <: SimpleSDMLayer} + Base.eltype(layer::T) where {T <: SimpleSDMLayer} Returns the type of the values stored in the grid. """ -function Base.eltype(p::T) where {T <: SimpleSDMLayer} - return eltype(p.grid) +function Base.eltype(layer::T) where {T <: SimpleSDMLayer} + return eltype(layer.grid) end """ - Base.size(p::T) where {T <: SimpleSDMLayer} + Base.size(layer::T) where {T <: SimpleSDMLayer} Returns the size of the grid. """ -function Base.size(p::T) where {T <: SimpleSDMLayer} - return size(p.grid) +function Base.size(layer::T) where {T <: SimpleSDMLayer} + return size(layer.grid) end """ - Base.size(p::T, i...) where {T <: SimpleSDMLayer} + Base.size(layer::T, i...) where {T <: SimpleSDMLayer} Returns the size of the grid alongside a dimension. """ -function Base.size(p::T, i...) where {T <: SimpleSDMLayer} - return size(p.grid, i...) +function Base.size(layer::T, i...) where {T <: SimpleSDMLayer} + return size(layer.grid, i...) end """ - Base.stride(p::T; dims::Union{Nothing,Integer}=nothing) where {T <: SimpleSDMLayer} + Base.stride(layer::T; dims::Union{Nothing,Integer}=nothing) where {T <: SimpleSDMLayer} Returns the stride, *i.e.* the length, of cell dimensions, possibly alongside a side of the grid. """ -function Base.stride(p::T; dims::Union{Nothing,Integer}=nothing) where {T <: SimpleSDMLayer} - lon_stride = (p.right-p.left)/size(p, 2)/2.0 - lat_stride = (p.top-p.bottom)/size(p, 1)/2.0 +function Base.stride(layer::T; dims::Union{Nothing,Integer}=nothing) where {T <: SimpleSDMLayer} + lon_stride = (layer.right-layer.left)/size(layer, 2)/2.0 + lat_stride = (layer.top-layer.bottom)/size(layer, 1)/2.0 dims == nothing && return (lon_stride, lat_stride) dims == 1 && return lon_stride dims == 2 && return lat_stride end """ - Base.eachindex(p::T) where {T <: SimpleSDMLayer} + Base.eachindex(layer::T) where {T <: SimpleSDMLayer} Returns the index of the grid. """ -function Base.eachindex(p::T) where {T <: SimpleSDMLayer} - return eachindex(p.grid) +function Base.eachindex(layer::T) where {T <: SimpleSDMLayer} + return eachindex(layer.grid) end """ Extracts a value from a layer by its grid position. """ -function Base.getindex(p::T, i::Int64) where {T <: SimpleSDMLayer} - return p.grid[i] +function Base.getindex(layer::T, i::Int64) where {T <: SimpleSDMLayer} + return layer.grid[i] end """ - Base.getindex(p::T, i::R, j::R) where {T <: SimpleSDMLayer, R <: UnitRange} + Base.getindex(layer::T, i::R, j::R) where {T <: SimpleSDMLayer, R <: UnitRange} Extracts a series of positions in a layer, and returns a layer corresponding to the result. This is essentially a way to rapidly crop a layer to a given subset @@ -105,19 +105,19 @@ as its argument, but this can be changed using `convert`. Note that this functio performs additional checks to ensure that the range is not empty, and to also ensure that it does not overflows from the size of the layer. """ -function Base.getindex(p::T, i::R, j::R) where {T <: SimpleSDMLayer, R <: UnitRange} +function Base.getindex(layer::T, i::R, j::R) where {T <: SimpleSDMLayer, R <: UnitRange} i_min = isempty(i) ? max(i.start-1, 1) : i.start - i_max = isempty(i) ? max(i.stop+2, size(p, 1)) : i.stop + i_max = isempty(i) ? max(i.stop+2, size(layer, 1)) : i.stop j_min = isempty(j) ? max(j.start-1, 1) : j.start - j_max = isempty(j) ? max(j.stop+2, size(p, 2)) : j.stop + j_max = isempty(j) ? max(j.stop+2, size(layer, 2)) : j.stop i_fix = i_min:i_max j_fix = j_min:j_max return T( - p.grid[i_fix,j_fix], - minimum(longitudes(p)[j_fix])-stride(p)[1], - maximum(longitudes(p)[j_fix])+stride(p)[1], - minimum(latitudes(p)[i_fix])-stride(p)[2], - maximum(latitudes(p)[i_fix])+stride(p)[2] + layer.grid[i_fix,j_fix], + minimum(longitudes(layer)[j_fix])-stride(layer)[1], + maximum(longitudes(layer)[j_fix])+stride(layer)[1], + minimum(latitudes(layer)[i_fix])-stride(layer)[2], + maximum(latitudes(layer)[i_fix])+stride(layer)[2] ) end @@ -125,105 +125,105 @@ end Given a layer and a latitude, returns NaN if the latitude is outside the range, or the grid index containing this latitude if it is within range """ -function _match_latitude(p::T, l::K) where {T <: SimpleSDMLayer, K <: AbstractFloat} - l > p.top && return NaN - l < p.bottom && return NaN - return findmin(abs.(l .- latitudes(p)))[2] +function _match_latitude(layer::T, lat::K) where {T <: SimpleSDMLayer, K <: AbstractFloat} + lat > layer.top && return NaN + lat < layer.bottom && return NaN + return findmin(abs.(lat .- latitudes(layer)))[2] end """ Given a layer and a longitude, returns NaN if the longitude is outside the range, or the grid index containing this longitude if it is within range """ -function _match_longitude(p::T, l::K) where {T <: SimpleSDMLayer, K <: AbstractFloat} - l > p.right && return NaN - l < p.left && return NaN - return findmin(abs.(l .- longitudes(p)))[2] +function _match_longitude(layer::T, lon::K) where {T <: SimpleSDMLayer, K <: AbstractFloat} + lon > layer.right && return NaN + lon < layer.left && return NaN + return findmin(abs.(lon .- longitudes(layer)))[2] end """ - Base.getindex(p::T, longitude::K, latitude::K) where {T <: SimpleSDMLayer, K <: AbstractFloat} + Base.getindex(layer::T, longitude::K, latitude::K) where {T <: SimpleSDMLayer, K <: AbstractFloat} Extracts the value of a layer at a given latitude and longitude. If values outside the range are requested, will return `NaN`. """ -function Base.getindex(p::T, longitude::K, latitude::K) where {T <: SimpleSDMLayer, K <: AbstractFloat} - i = _match_longitude(p, longitude) - j = _match_latitude(p, latitude) +function Base.getindex(layer::T, longitude::K, latitude::K) where {T <: SimpleSDMLayer, K <: AbstractFloat} + i = _match_longitude(layer, longitude) + j = _match_latitude(layer, latitude) isnan(i) && return NaN isnan(j) && return NaN - return p.grid[j, i] + return layer.grid[j, i] end """ - Base.getindex(p::T; left=nothing, right=nothing, top=nothing, bottom=nothing) where {T <: SimpleSDMLayer, K <: Union{Nothing,AbstractFloat}} + Base.getindex(layer::T; left=nothing, right=nothing, top=nothing, bottom=nothing) where {T <: SimpleSDMLayer, K <: Union{Nothing,AbstractFloat}} Returns a subset of the argument layer, where the new limits are given by `left`, `right`, `top`, and `bottom`. Up to three of these can be omitted, and if so these limits will not be affected. """ -function Base.getindex(p::T; left=nothing, right=nothing, top=nothing, bottom=nothing) where {T <: SimpleSDMLayer} +function Base.getindex(layer::T; left=nothing, right=nothing, top=nothing, bottom=nothing) where {T <: SimpleSDMLayer} for limit in [left, right, top, bottom] if !isnothing(limit) @assert typeof(limit) <: AbstractFloat end end - imax = _match_longitude(p, isnothing(right) ? p.right : right) - imin = _match_longitude(p, isnothing(left) ? p.left : left) - jmax = _match_latitude(p, isnothing(top) ? p.top : top) - jmin = _match_latitude(p, isnothing(bottom) ? p.bottom : bottom) + imax = _match_longitude(layer, isnothing(right) ? layer.right : right) + imin = _match_longitude(layer, isnothing(left) ? layer.left : left) + jmax = _match_latitude(layer, isnothing(top) ? layer.top : top) + jmin = _match_latitude(layer, isnothing(bottom) ? layer.bottom : bottom) any(isnan.([imin, imax, jmin, jmax])) && throw(ArgumentError("Unable to extract, coordinates outside of range")) - return p[jmin:jmax, imin:imax] + return layer[jmin:jmax, imin:imax] end """ - Base.getindex(p::T, n::NT) where {T <: SimpleSDMLayer, NT <: NamedTuple} + Base.getindex(layer::T, n::NT) where {T <: SimpleSDMLayer, NT <: NamedTuple} Returns a subset of the argument layer, where the new limits are given in a NamedTuple by `left`, `right`, `top`, and `bottom`, in any order. Up to three of these can be omitted, and if so these limits will not be affected. """ -function Base.getindex(p::T, n::NT) where {T <: SimpleSDMLayer, NT <: NamedTuple} +function Base.getindex(layer::T, n::NT) where {T <: SimpleSDMLayer, NT <: NamedTuple} l = isdefined(n, :left) ? n.left : nothing r = isdefined(n, :right) ? n.right : nothing t = isdefined(n, :top) ? n.top : nothing b = isdefined(n, :bottom) ? n.bottom : nothing - Base.getindex(p; left = l, right = r, top = t, bottom = b) + Base.getindex(layer; left=l, right=r, top=t, bottom=b) end """ - Base.getindex(p1::T1, p2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer} + Base.getindex(layer1::T1, layer2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer} Extract a layer based on a second layer. Note that the two layers must be *compatible*, which is to say they must have the same bounding box and grid size. """ -function Base.getindex(p1::T1, p2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer} - SimpleSDMLayers.are_compatible(p1, p2) - return p1[left=p2.left, right=p2.right, bottom=p2.bottom, top=p2.top] +function Base.getindex(layer1::T1, layer2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer} + SimpleSDMLayers._layers_are_compatible(layer1, layer2) + return layer1[left=layer2.left, right=layer2.right, bottom=layer2.bottom, top=layer2.top] end """ - Base.setindex!(p::SimpleSDMResponse{T}, v::T, i...) where {T} + Base.setindex!(layer::SimpleSDMResponse{T}, v::T, i...) where {T} Changes the value of a cell, or a range of cells, as indicated by their grid positions. """ -function Base.setindex!(p::SimpleSDMResponse{T}, v::T, i...) where {T} - @assert typeof(v) <: eltype(p.grid) - p.grid[i...] = v +function Base.setindex!(layer::SimpleSDMResponse{T}, v::T, i...) where {T} + typeof(v) <: eltype(layer.grid) || throw(ArgumentError("Impossible to set a value to a non-matching type")) + layer.grid[i...] = v end """ - Base.setindex!(p::T, v, lon::Float64, lat::Float64) where {T <: SimpleSDMResponse} + Base.setindex!(layer::T, v, lon::Float64, lat::Float64) where {T <: SimpleSDMResponse} Changes the values of the cell including the point at the requested latitude and longitude. """ -function Base.setindex!(p::SimpleSDMResponse{T}, v::T, lon::Float64, lat::Float64) where {T} - i = _match_longitude(p, lon) - j = _match_latitude(p, lat) - p[j,i] = v +function Base.setindex!(layer::SimpleSDMResponse{T}, v::T, lon::Float64, lat::Float64) where {T} + i = _match_longitude(layer, lon) + j = _match_latitude(layer, lat) + layer[j,i] = v end """ @@ -231,16 +231,16 @@ end Returns a `SimpleSDMResponse` of the same dimensions as the original layer, with `NaN` in the same positions. The rest of the values are replaced by the output -of `zero(eltype(p.grid))`, which implies that there must be a way to get a zero +of `zero(eltype(layer.grid))`, which implies that there must be a way to get a zero for the type. If not, the same result can always be achieved through the use of `copy`, manual update, and `convert`. """ -function Base.similar(l::T) where {T <: SimpleSDMLayer} - emptygrid = similar(l.grid) +function Base.similar(layer::T) where {T <: SimpleSDMLayer} + emptygrid = similar(layer.grid) for i in eachindex(emptygrid) - emptygrid[i] = isnan(l.grid[i]) ? NaN : zero(eltype(l.grid)) + emptygrid[i] = isnan(layer.grid[i]) ? NaN : zero(eltype(layer.grid)) end - return SimpleSDMResponse(emptygrid, l.left, l.right, l.bottom, l.top) + return SimpleSDMResponse(emptygrid, layer.left, layer.right, layer.bottom, layer.top) end """ @@ -248,7 +248,7 @@ end Returns a new copy of the layer, which has the same type. """ -function Base.copy(l::T) where {T <: SimpleSDMLayer} - copygrid = copy(l.grid) - return T(copygrid, copy(l.left), copy(l.right), copy(l.bottom), copy(l.top)) +function Base.copy(layer::T) where {T <: SimpleSDMLayer} + copygrid = copy(layer.grid) + return T(copygrid, copy(layer.left), copy(layer.right), copy(layer.bottom), copy(layer.top)) end diff --git a/src/lib/types.jl b/src/lib/types.jl index 319cfb16..aedfdb4a 100644 --- a/src/lib/types.jl +++ b/src/lib/types.jl @@ -26,6 +26,12 @@ struct SimpleSDMPredictor{T} <: SimpleSDMLayer top::AbstractFloat end +""" + SimpleSDMPredictor(grid::Matrix{T}) where {T} + +If only a matrix is given to `SimpleSDMPredictor`, by default we assume that it +covers the entire range of latitudes and longitudes. +""" function SimpleSDMPredictor(grid::Matrix{T}) where {T} return SimpleSDMPredictor(grid, -180., 180., -90., 90.) end @@ -42,6 +48,12 @@ mutable struct SimpleSDMResponse{T} <: SimpleSDMLayer top::AbstractFloat end +""" + SimpleSDMResponse(grid::Matrix{T}) where {T} + +If only a matrix is given to `SimpleSDMPredictor`, by default we assume that it +covers the entire range of latitudes and longitudes. +""" function SimpleSDMResponse(grid::Matrix{T}) where {T} return SimpleSDMResponse(grid, -180., 180., -90., 90.) end diff --git a/src/recipes/recipes.jl b/src/recipes/recipes.jl index e0c82f16..558d0084 100644 --- a/src/recipes/recipes.jl +++ b/src/recipes/recipes.jl @@ -22,7 +22,7 @@ test 2 if get(plotattributes, :seriestype, :scatter) in [:scatter, :histogram2d] @assert eltype(l1) <: Number @assert eltype(l2) <: Number - SimpleSDMLayers.are_compatible(l1, l2) + SimpleSDMLayers._layers_are_compatible(l1, l2) valid_i = filter(i -> !(isnan(l1[i])|isnan(l2[i])), eachindex(l1.grid)) l1.grid[valid_i], l2.grid[valid_i] end