poisotlab / simplesdmlayers.jl Goto Github PK

According to @gabrieldansereau in #21 it looks like the size is inverted:

julia> latitudes(temperature) |> length
2160

julia> longitudes(temperature) |> length
1080

julia> stride(temperature)
(0.16666666666666666, 0.04158950617283951)

In particular, the stride should have the same values. This is because ArchGDAL can read things in a buffer with different dimensions, and this is in the process of being fixed.

It would be nice to have some built-in support for DataFrames, a bit like the clip() function for GBIF occurrences and possibly a convert() method to easily switch between layers and DataFrames. I'm thinking of some kind of optional dependency, like for GBIF (if possible, not sure if it can easily be done). I have some code I'll try to port here as soon as I can.

Get the values by removing the nothing and return as an array

Currently, we can extract a layer given a range of integers, which are coordinates on the grid. For most applications, we want to extract given a range of latitudes / longitudes. This would be a getindex method with a declaration like

getindex(p::T, longitudes::Tuple{K,K}, latitudes::Tuple{K,K})

This would return a SimpleSDMLayer (of type T), where the endpoints are contained in the grid. Note that because the grid is discrete, the endpoints given in arguments will not be the endpoints of the actual layer.

Question for @mkborregaard

I've added a recipe for a scatterplot (scatter(l1, l2) plots the values of the cells in two layers, and for large maps it would be better to have this as 2d histogram -- is this doable through recipes?

(should be simple to add a recipe)

expanding on #59

one thing i found myself needing is wanting to create a SDMLayer that contained 1 at all locations with given occurrence data, and 0 elsewhere. the not optimized example i have is

getOccupancyLayer(envLayer::SimpleSDMPredictor, occupancy) = begin
    occLayer = similar(envLayer)
    for o in occupancy
        long,lat = longitude(o), latitude(o) 
        occLayer[long,lat] = convert(eltype(occLayer), 1)
    end
    occLayer
end

_pixel_score(x) = x >= 0.5 ? 2.0(1.0-x) : 2.0x

and then add

prediction = BIOCLIM(predictors, models)

rescale!(prediction, collect(0.0:0.05:1.0))

cutoff = broadcast(x -> x > 0.05, prediction)

plot(prediction, c=:lightgrey)
mask(cutoff, prediction) |> plot!
scatter!(longitudes(records), latitudes(records), lab="")

When reading from geotiff, the nodata value can be ignored.

Something to create a similar layer but with a different type (e.g. same as a worldclime layer but with Bool values)

similar(layer::ST, T) where {ST <: SimpleSDMLayer}

or something to that effect.

@tpoisot I ran into this when trying to update some older scripts. Subsetting a layer based on coordinates does not return the same dimensions as subsetting when reading from a tif file. I would expect the two options to return the same thing.

coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)

t1 = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...)
t2 = SimpleSDMPredictor(WorldClim, BioClim, 1)[coords]

size(t1) # 331 x 571
size(t2) # 332 x 571

The dimensions of t1 match what I had in my earlier scripts.

I found that the breaking change occurred a while ago between v0.4.2 and v0.4.3, more specifically in #58 where you did many edits to _match_latitude() and those functions. t2 uses _match_latitude() and getindex() for the subsetting, while t1 uses geotiff().

I'm not sure how to properly fix this. I tried playing with the stride in _match_latitude() (something like last(findmin(abs.(lat .- latitudes(layer) **.- stride(layer)**)), or even using range instead of LinRange in latitudes(). Both fix my problem but cause the tests to fail, and surely break other things.

You tried many things in #58, do you see how to fix this?

Part of the plotting recipes involves adding NaN, which means that we need to have something that is floating-point. This implies that there is a need to convert the layer to something more amenable to plotting in the recipe. This can be done with a more general convert method to change the inner type.

I don't see any reference to community guidelines for people wishing to contribute. This is the text in the reviewer checklist for reference:

Community guidelines: Are there clear guidelines for third parties wishing to 1) Contribute to the software 2) Report issues or problems with the software 3) Seek support

Hi folks, I had a few feature ideas that have come up when using the package, happy to help implement any/all of them if they seem good.

A feature I needed was to convert a set of layers and GBIF occupancy points to a set of features (a matrix of environmental variables at each point in the layer) with corresponding labels (0 for points in the layer without occurrence, 1 at each point with occurrence).

It seems like this feature could have widespread utility for using this with Turing (as I did) or Flux.

My (crude) example code below --- apologies if this already exists in some form

function features(environment, occurrence) 

    xDim, yDim = size(environment[1])
    numberSpatialPoints = xDim*yDim
    numFeatures = length(environment)
        
    featuresMatrix = zeros(numberSpatialPoints, numFeatures)
    labels = [false for i in 1:numberSpatialPoints]
    occupancyLayer = getOccupancyLayer(environment[1], occurrence)
    
    cursor = 1
    for pt in 1:numberSpatialPoints
        if (!isnothing(environment[1][pt] ))
            for f in 1:numFeatures
                featuresMatrix[cursor, f] = environment[f][pt]  
                labels[cursor] = occupancyLayer[pt]
            end
            cursor += 1
        end
    end

    return (featuresMatrix[1:cursor, :], labels[1:cursor])
end

Unless there's something I don't see, our current overload for getindex(layer1, layer2) (used when calling layer1[layer2]) can only return exactly the same layer as layer1, which makes it a bit useless.

Calling _layers_are_compatible(layer1, layer2) means the overload only works when the two layers have the same size AND bounding coordinates. But then we return layer1 with the bounding coordinates of layer2, which are necessarily the same as layer[1].

I think this overload would be more useful if it allowed to subset a layer based on another layer if both have the same stride and the 2nd one is entirely contained in the first one, as in the example below.

l1, l2 = SimpleSDMPredictor(WorldClim, BioClim, 1:2)
l3 = SimpleSDMPredictor(WorldClim, BioClim, 1; left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
stride(l1) == stride(l3) # same stride

l4 = l1[l2] # works
l4 == l1 # but it's exactly the same

l5 = l1[l3] # doesn't work as layers have different sizes, but should work
l5 == l3 # should return true

This is the template in Documenter.jl -- might be a good idea to have Plots as a dependency for the documentation too...

jobs:
  include:
    - stage: "Documentation"
      julia: 1.0
      os: linux
      script:
        - julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd()));
                                               Pkg.instantiate()'
        - julia --project=docs/ docs/make.jl
      after_success: skip

Perform a simple check on bounding box and stride, and return the result, this will help with tiles

Love the new mask function @tpoisot ! It will be incredibly useful.

I think it would be nice if it worked on DataFrames too, not just on GBIFrecords elements. That would be in line with the previous DataFrames integration I've added, where I mostly mimicked the features from the GBIF integration so that both can be used in a similar way.

I created this function for my project a while ago and I think it works well enough. I'll update it to be similar to mask and add it in a PR, if you agree.

continuing on #59

latitude(r::GBIFRecord) = r.latitude longitude(r::GBIFRecord) = r.longitude
as an interface was nice, allthough latitudes/longitudes are defined already for layers

This is related to openjournals/joss-reviews#2872.

So I've managed to install the package, but running through some of the examples in the docs I get this error. Do you have any idea what might be going? Something on my end or something wrong with the site?

julia> precipitation = worldclim(12)
[ Info: Downloading wc2.0_bio10m.zip
ERROR: HTTP.ExceptionRequest.StatusError(404, "GET", "/data/worldclim/v2.0/tif/base/wc2.0_10m_bio.zip", HTTP.Messages.Response:
"""
HTTP/1.1 404 Not Found
Date: Thu, 17 Dec 2020 17:42:41 GMT
Server: Apache/2.4.18 (Ubuntu)
Last-Modified: Mon, 12 Oct 2020 23:20:08 GMT
ETag: "df-5b1818a0b093c;5b18143b1d0ec"
Accept-Ranges: bytes
Content-Length: 223
Content-Type: text/html

<!DOCTYPE html>
<html>
<head>
  <meta http-equiv="Refresh" content="5; URL=https://worldclim.org/">
</head>

<body>

<h1>Error Loading Page (404)</h1>

The page or file you are looking for is not available.

</body>
</html>""")

This would be super convenient (and easy to write). Something like

write(layer, file; left, right, bottom, top)

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

The reason is that the layers are too big for git - the fix is simple, we need to remove the assets folder at the end of the doc. This should be a feature in SimpleSDMDataSources as well.

@mkborregaard suggested an iterator in #48 , which is a good idea

According to my (limited) understanding of the recipes system, calling heatmap should set the correct values for x and y, through the following lines:

https://github.com/EcoJulia/SimpleSDMLayers.jl/blob/6b2bdf6676f8588e4b5196e13d9c83014bcbaf24/src/recipes/recipes.jl#L5-L7

And yet, when calling it with a layer, as in

https://github.com/EcoJulia/SimpleSDMLayers.jl/blob/6b2bdf6676f8588e4b5196e13d9c83014bcbaf24/test/plots.jl#L15

I get a figure with the correct z, but the wrong x and y (basically the positions):

So clearly, there's a step I'm missing -- @mkborregaard do you know the right incantation to get this to work?

While fixing #100, I realized geotiff() also modifies the source layer when writing, even without re-reading. Probably just because there's a call to replace!() in the function, which modifies the layer in place.

l1 = SimpleSDMPredictor(WorldClim, BioClim, 1)
length(l1) # 808053

geotiff(tempname(), l1)
length(l1) # 2332800

All of this discussion about CHELSA in #31, as well as the amount of data managing code I had to write for #21, makes me wonder if we would not be better off with a SimpleSDMDataSources.jl package. In term, I'd like to offer access to more than bioclim on the CHELSA side, and maybe other earthenv layers as well.

@gabrieldansereau and @mkborregaard what are your opinion on that? I don't think this tentative new package would be used on its own, but would instead be a way to separate the data download from the layer manipulation code.

geotiff does not use the specified nodata value when writing to a tif file, whether it's the default -9999 or a custom one. It sets the nodata value for the file with ArchGDAL.setnodatavalue! correctly, but does not use it in the grid itself and uses NaN instead.

Because of this, it is always necessary to call replace(geotiff(SimpleSDMPredictor, file), NaN => nothing) when re-reading a layer after writing it, which in my sense could be avoided.

l = SimpleSDMPredictor(WorldClim, BioClim, 1)

# With -9999 as nodata
f1 = tempname()
geotiff(f1, l) # use -9999 as nodata by default
mp1 = geotiff(SimpleSDMResponse, f1)

length(mp1) == length(l) # false, but should be true
mp1 == l # false, same thing
replace(mp1, NaN => nothing) == l1 # true, but could be avoided

# With custom nodata
f2 = tempname()
geotiff(f2, l; nodata=-3.4f38) # this is the nodata from the WorldClim files
mp2 = geotiff(SimpleSDMPredictor, f2)

length(mp2) == length(l) # false
mp2 == l # false
replace(mp2, NaN => nothing) == l1 # true

With ArchGDAL it is possible to only load a subset of the grid, not the entire grid. Because we know the total range, it should be possible to use something like

datasource(layer; left=10.0, right=20.0, bottom=0.0, top=30.0)

to only get a subset. This would result in allocating a smaller array, which is definitely helpful when working with either landcover or bioclim.

Technically all of the pieces are already there since the package does this sort of operations when cropping.

This will add Requires as a dependency, and would greatly facilitate all of the GBIF-based pipelines.

Tagging @mkborregaard

I added support for plots in this comment: 41c4b6b

So far, this works with density and histogram (show the values) and heatmap (show the map). I am not sure how to proceed to make it work so that plot defaults to a heatmap - any tips?

hcat and vcat return inversed bottom & top bounding coordinates. Because of this, the results of latitudes (and other functions which depend on it) may be wrong. Somehow, this created a circular error and the tests still passed, but it is clearly not the intended behaviour.

Example from the tests:

julia> l1 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=0.0, top=10.0);
julia> l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=10.0, top=20.0);
julia> l3 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=10.0, right=20.0, bottom=0.0, top=10.0);

julia> vl1 = vcat(l1, l2) # latitudes span makes no sense
SDM predictor → 122×61 grid with 5573 Float32-valued cells
  Latitudes     (9.834016393442623, 9.99931693989071)
  Longitudes    (-0.08333333333333333, 9.916666666666666)

julia> ml1 = hcat(l1, l3) # things look fine here, but are they?
SDM predictor → 61×122 grid with 5573 Float32-valued cells
  Latitudes     (-0.08333333333333333, 9.916666666666666)
  Longitudes    (-0.08401639344262293, 19.917349726775956)

julia> latitudes(vl1) # span really makes no sense
122-element LinRange{Float64}:
 9.99932,9.99795,9.99658,9.99522,9.99385,9.99249,…,9.83948,9.83811,9.83675,9.83538,9.83402

julia> vl1.bottom < vl1.top # false
false

julia> latitudes(ml1) # latitudes are reversed??
61-element LinRange{Float64}:
 9.91667,9.75,9.58333,9.41667,9.25,9.08333,8.91667,…,0.583333,0.416667,0.25,0.0833333,-0.0833333

julia> ml1.bottom < ml1.top # false
false

Add a rescale/rescale! function to give a layer the first extrema as another layer

#102 made me realize how easily the bounding coordinates can be inversed by mistake, and how they can pass silently for the same layers.

I don't really see any use in allowing constructing layers with inversed bounds, where top < bottom & right < left, so I'll add an argument error to prevent this.

l1 = SimpleSDMPredictor(WorldClim, BioClim, 1)
l2 = SimpleSDMPredictor(copy(l1.grid), 180.0, -180.0, 90.0, -90.0) # looks the same, but bounds really are inversed
# This should throw an error

l1.grid == l2.grid # same grid
extrema(longitudes(l1)) == extrema(longitudes(l2)) # same longitude extremas (displayed by show method)
extrema(latitudes(l1)) == extrema(latitudes(l2)) # same latitude extremas too

longitudes(l1) == longitudes(l2) # false, inversed
latitudes(l1) == latitudes(l2) # false, inversed

Edit: Fixed in a40267e and 77eca5d

Add union, intersect, setdiff functions for layers

I'm trying to use Revise.jl to make the changes discussed in #74, and I'm running into some issues. I'm new to Revise though.

For now, calling using Revise before using SimpleSDMLayers in a dev project works great when I make changes to any functions from the package, except for the GBIF/DataFrames integrations. I know they are loaded through Requires.jl, but this is supposed to be supported by Revise (see here)

I managed to make it work by removing the joinpath() in https://github.com/EcoJulia/SimpleSDMLayers.jl/blob/db148a57a597dbb7412bae50b2c0168815ad0df1/src/SimpleSDMLayers.jl#L55-L65 to

function __init__()
    @require GBIF="ee291a33-5a6c-5552-a3c8-0f29a1181037" begin
        @info "GBIF integration loaded"
        include("integrations/GBIF.jl")
    end
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        @info "DataFrames integration loaded"
        include("integrations/DataFrames.jl")
    end

end

@tpoisot I know your preference for joinpath, but do you mind leaving that one? Is there another fix?

Will do as part of #90

Probably the same problem as in #99, but I ran into this in a different context so I'll post it in a new issue.

Writing a layer to a .tif file using geotiff() and then re-loading it (again with geotiff) can sometimes return a layer with different dimensions, as in this example:

using SimpleSDMLayers
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)

ref = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...)
geotiff("test.tif", ref)
test = geotiff(SimpleSDMPredictor, "test.tif")

isequal(size(ref), size(test)) # not equal, 331 x 571 vs 331 x 570

isequal(ref.left, test.left) # equal
isequal(ref.bottom, test.bottom) # equal
isequal(ref.top, test.top) # equal
isequal(ref.right, test.right) # not equal

ref.right # -50.0
test.right # -50.16666666665

Looking at the exported tif file, the size and bounding coordinates look fine, so the problem is probably with the reading geotiff call (which returned correct dimensions in #99 on the other hand).

using ArchGDAL
d = ArchGDAL.read("test.tif") # correct raster size, 331 x 571
gt = ArchGDAL.getgeotransform(d)
gt[1] + gt[2]*ArchGDAL.width(d) # correct right bound, ≈ 50.0

In some functions, the layers are called p, and in other they are called l. There isn't any reason to have both, so this is good argument to refactor (preferably to something like layer, just to be extra explicit).

CSV.write refuses to write nothing to a file - I think it would be acceptable to remove all rows with a nothing in values, right?

It needs to be fully done

The following line causes coarsen to crash when the element type of the grid is not floating points:

https://github.com/EcoJulia/SimpleSDMLayers.jl/blob/72db07ab55f596fe1c5ed8e657b37918ed890b21/src/operations/coarsen.jl#L30

A use case is: every cell stores a Set of species name, and we want to do a union of multiple cells to aggregate them over space. This doesn't work. This requires a new test, and to add a condition to only performs this check if the eltype(V) <: AbstractFloat.

#74 made me realize that there is no overload for replace or replace!, which would be good to have

As an alternative to coarsen, it could be really good to implement a sliding window approach, which can be either using blocks, or using all pixels below a certain Haversine distance.

Idea for an example: get all data from Hypomyces lactifluorum, and make the predictions in Québec - this can be done with a random forest, so that's a nice addition to the doc.

@mkborregaard I've solved the issues related to GDAL updates, and to curl doing weird things on some machines (by using HTTP instead) -- I think this is usable for light explorations (@gabrieldansereau is using it at the moment), what is your opinion on a v0.0.1 version?

A couple of ideas that could go in the next patch release. I'll work on them soon-ish

• Add an empty(layer) overload, returning something like SimpleSDMPredictor(fill(nothing, size(layer)), layer)
• Allow some methods to work on AbstractDataFrames instead of DataFrames only, which would be useful when working with GroupedDataFrames grouped on species
• mask
  • Fix the docstring for DataFrames
  • Add some info in the docstring mentioning that, to produce a GBIF/DataFrame range layer, we need to use a full continental background, not a pre-existing range layer of some sort
  • Possibly add an unmask function or something doing the opposite of mask: adding zeros where another has defined values
  • Try to speed up a bit
  • Fix mask(::SimpleSDMLayer, ::DataFrame) when the DataFrame contains observations outside the layer's range (by filtering out the coordinates) (#111)
• Add the union, intersect, setdiff overloads: probably just for Bool layers though (#94)
• Try to speed up mosaic, which is suuuuper slow
• Add geotiff method to call bands to load with a UnitRange or Vector, as for the SimpleSDMPredictor methods
• Fix the convert(::T, layer::SimpleSDMPredictor) method to return a SimpleSDMPredictor, not a SimpleSDMResponse (or update the docstring to match the current behaviour)
• Add a convert(Vector, layer) or vec(layer) overload. Or even a convert(Matrix, layers::Vector{SimpleSDMLayer}) grid, #112
• Allow to write band names with geotiff
• Figure out a way to write Bool layers with geotiff
• xlim & ylim on plot recipe should be based on boundinxbox(layer), not extrema(latitudes(layer))
• Speed up the future data loading calls
• Clarify error message for plot(layer::SimpleSDMLayer{Int64})
• Try to fix geotiff writing with Int8 & Float16
• Add a way to set layer names (a new field for the types?)
• Add instructions to set SDMLAYERS_PATH
• Add info message when downloading data
• Add info for the conversion factor for CHELSA
• Add a more informative error message when trying to load a file which does not exist with geotiff (currently returns an error from GDAL/ArchGDAL GDALError (CE_Failure, code 10): Pointer 'hDS' is NULL in 'GDALGetGeoTransform')

The cellsize is half of what it should be

These three are supposedly the same latitudes/longitudes (after #23) , so this is concerning - this is blocking #21

Interestingly, the longitudes seem to be right, and it's the latitudes that get shifted. This suggests that the problem is with finding the part of the buffer to read in geotiff. Alternatively, worldclim seems to get the latitudes right, so maybe there is some weird GDAL thing going on with bioclim and landcover, which would require to dig into the documentation.

https://www.earthenv.org/landcover