Merge pull request #21 from EcoJulia/feature/landcover

EarthEnv landcover + WorldClim 2.1 + CHELSA V1

.gitignore

@@ -2,6 +2,7 @@
 *.jl.cov
 *.jl.*.cov
 
+assets
 test/assets
 test/gallery
 

Project.toml

@@ -1,18 +1,18 @@
 name = "SimpleSDMLayers"
 uuid = "2c645270-77db-11e9-22c3-0f302a89c64c"
 authors = ["Timothée Poisot <[email protected]>"]
-version = "0.2.2"
+version = "0.3.0"
 
 [deps]
-GDAL = "add2ef01-049f-52c4-9ee2-e494f65e021a"
+ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
 
 [compat]
-GDAL = "1.0"
+ArchGDAL = "0.4"
 HTTP = "0.8"
 RecipesBase = "0.7, 0.8, 1.0"
 Requires = "1.0"

docs/Project.toml

@@ -3,3 +3,8 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 GBIF = "ee291a33-5a6c-5552-a3c8-0f29a1181037"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+
+[compat]
+GBIF = "0.2"

docs/make.jl

@@ -2,6 +2,7 @@ push!(LOAD_PATH, joinpath("..", "src"))
 
 using Documenter, SimpleSDMLayers
 using GBIF
+using Statistics
 
 makedocs(
     sitename = "Simple SDM Layers",
@@ -17,7 +18,11 @@ makedocs(
         ],
         "Examples" => [
             "Temperature data" => "examples/temperature.md",
-            "GBIF integration" => "examples/gbif.md"
+            "GBIF integration" => "examples/gbif.md",
+            "Importing raster data" => "examples/import.md",
+            "Sliding window analysis" => "examples/slidingwindow.md",
+            "Landcover data" => "examples/landcover.md",
+            "Landcover consensus" => "examples/consensus.md"
         ]
     ]
 )

docs/src/examples/consensus.md

@@ -0,0 +1,83 @@
+# Landcover consensus
+
+In this example, we will create a consensus map of landcover for Corsica based
+on the EarthEnv data, and measure the variation within each pixel using the
+variance. The first step is to load the packages we need, and create a bounding
+box:
+
+```@example cons
+using SimpleSDMLayers
+using Plots
+
+bbox = (left=8.25, right=10.0, bottom=41.2, top=43.2)
+```
+
+We will then do two things. First, get the first layer of landcover (see the
+help of `landcover` for a list of the layers), and then create a datacube,
+organized around dimensions of latitude, longitude, and layer value - we will
+only focus on the 11 first variables, since we do not want the information on
+open water (layer 12):
+
+```@example cons
+lc = landcover(1; full=false, bbox...)
+use = fill(NaN32, size(lc)..., 11)
+```
+
+At this point, we will simply fill in the first "slice" of our datacube with
+values from the layer:
+
+```@example cons
+for (i,e) in enumerate(lc.grid)
+    coord = (CartesianIndices(size(lc.grid))[i].I..., 1)
+    if !isnothing(e)
+        use[coord...] = e
+    end
+end
+```
+
+The next step is to repeat this process for all other layers, filling the
+appropriate data cube slice:
+
+```@example cons
+for layer in 2:11
+    lc = landcover(layer; full=false, bbox...)
+    for (i,e) in enumerate(lc.grid)
+        coord = (CartesianIndices(size(lc.grid))[i].I..., layer)
+        if !isnothing(e)
+            use[coord...] = e
+        end
+    end
+end
+```
+
+To perform the actual analysis, we will define a `get_most_common_landuse` function, which will return the index of the layer with the highest score:
+
+```@example cons
+function get_most_common_landuse(f)
+    f[isnan.(f)] .= 0.0
+    sum(f) == 0 && return NaN
+    return last(findmax(f))
+end
+
+function shannon(x)
+    v = filter(!isnan, x)
+    length(v) == 0 && return NaN
+    v = v ./ sum(v)
+    return -sum(v.*log2.(v))
+end
+```
+
+```@example cons
+consensus = mapslices(get_most_common_landuse, use; dims=3)[:,:,1]
+entropy = mapslices(shannon, use; dims=3)[:,:,1]
+
+consensus = SimpleSDMResponse(consensus, lc)
+entropy = SimpleSDMResponse(entropy, lc)
+```
+
+```@example cons
+p1 = plot(consensus, c=:Paired_11, frame=:grid)
+p2 = plot(entropy, c=:Greys, frame=:grid)
+
+plot(p1, p2, size=(900, 400))
+```

docs/src/examples/import.md

@@ -0,0 +1,64 @@
+# Importing your own data
+
+It is possible to import your own rasters into a `SimpleSDMLayer` object. This
+requires defining a new type and two "helper" functions, which might seem a
+little bit convoluted, but helps *immensely* underneath in case you want to also
+*download* rasters from the web with different arguments. In this example, we
+will look at a data file produced by the *OmniScape* package, and which
+represents landscape connectivity in the Laurentians region of Québec. This
+example will also show how we can use the `broadcast` operation to modify the
+values of a raster.
+
+```@example temp
+using SimpleSDMLayers
+using Plots
+using StatsBase
+```
+
+The file comes with the package itself, so we can read it directly - this is a
+geotiff file, where values are floating point numbers representing connectivity.
+
+```@example temp
+file = joinpath(dirname(pathof(SimpleSDMLayers)), "..", "data", "connectivity.tiff")
+```
+
+To import this file as a `SimpleSDMLayer`, we need to create a type
+(`MyConnectivityMap`), and declare a method for `latitudes` and `longitudes` for
+this type, where the output is the range of latitudes and longitudes. This might
+seem cumbersome, but remember: it can be automated, and if you do not declare a
+`latitude` and `longitude` method, it will be assumed that the raster covers the
+entire globe. From a end-user perspective, it also removes the need to pass the
+bounding box of your layer as an argument, and to focus instead of the region of
+interest.
+
+```@example temp
+struct MyConnectivityMap <: SimpleSDMLayers.SimpleSDMSource end
+SimpleSDMLayers.latitudes(::Type{MyConnectivityMap}) = (45.34523, 47.38457)
+SimpleSDMLayers.longitudes(::Type{MyConnectivityMap}) = (-75.17734,-72.36486)
+```
+
+Now that this is done, we can read this file as a `SimpleSDMResponse` using the
+`raster` function:
+
+```@example temp
+mp = SimpleSDMLayers.raster(SimpleSDMResponse, MyConnectivityMap(), file)
+```
+
+Because this file has raw values, which are not necessarily great for plotting,
+we will transform it to quantiles, using the `StatsBase.ecdf` function.
+
+```@example temp
+qfunc = ecdf(convert(Vector{Float64}, filter(!isnothing, mp.grid)))
+```
+
+And we can now broadcast this function to the layer:
+
+```@example temp
+qmap = broadcast(qfunc, mp)
+```
+
+Finally, we are ready for plotting:
+
+```@example temp
+plot(qmap, frame=:grid, c=:YlGnBu)
+```

docs/src/examples/landcover.md

@@ -0,0 +1,43 @@
+# Getting landcover data
+
+In this example, we will look at landcover data, specifically the proportion of
+urban/built area in Europe; the entire dataset is very large to fit in memory,
+as it has a resolution of about 1 kilometre squared. Therefore, we will take
+advantage of the ability to only load the part that matters by passing the
+limits of a bounding box.
+
+```@example urban
+using SimpleSDMLayers
+urban = landcover(9; left=-11.0, right=31.1, bottom=29.0, top=71.1)
+```
+
+This dataset is returning data as `UInt8` (as it represents a proportion of the
+pixel occupied by the type), but this is not something that can be plotted efficiently. So in the
+next step, we will manipulate this object a little bit to have something more
+workable.
+
+Let's start by preparing a new grid, with the same dimensions, but a friendlier
+type, and then we can then fill these values using a simple rule of using either
+`NaN` or the converted value:
+
+```@example urban
+n_urban_grid = zeros(Float32, size(urban));
+for (i,e) in enumerate(urban.grid)
+  n_urban_grid[i] = isnothing(e) ? NaN : Float32(e)
+end
+```
+
+We can now overwrite our `urban` object as a layer:
+
+```@example urban
+urban = SimpleSDMPredictor(n_urban_grid, urban)
+```
+
+Note that the previous instruction uses a shortcut where the bounding box from a
+new `SimpleSDMLayer` is drawn from the bounding box for an existing layer. With
+this done, we can show the results:
+
+```@example urban
+using Plots
+heatmap(urban, c=:terrain)
+```

docs/src/examples/slidingwindow.md

@@ -0,0 +1,30 @@
+# Sliding window analysis
+
+In this example, we will get precipitation data from Québec, and use a sliding
+window analysis to smooth them out. The beginning of the code should now be
+familiar:
+
+```@example slide
+using SimpleSDMLayers
+using Plots
+using Statistics
+
+precipitation = worldclim(12; left=-80.0, right=-56.0, bottom=44.0, top=62.0)
+```
+
+The sliding window works by taking all pixels *within a given radius* (expressed
+in kilometres) around the pixel of interest, and then applying the function
+given as the second argument to their values. Empty pixels are removed. In this
+case, we will do a summary across a 100 km radius around each pixel:
+
+```@example slide
+averaged = slidingwindow(precipitation, Statistics.mean, 100.0)
+```
+
+We can finally overlap the two layers -- the result of sliding window is a
+little bit smoother than the raw data.
+
+```@example slide
+plot(precipitation, c=:alpine)
+contour!(averaged, c=:white, lw=2.0)
+```

docs/src/examples/temperature.md

@@ -15,7 +15,7 @@ visualize these data:
 
 ```@example temp
 using Plots, StatsPlots
-heatmap(temperature, c=:magma, frame=:box)
+heatmap(temperature, c=:cividis, frame=:box)
 xaxis!("Longitude")
 yaxis!("Latitude")
 ```
@@ -33,7 +33,7 @@ can also make a quick heatmap to see what the region looks like:
 
 ```@example temp
 temperature_europe = temperature[left=-11.0, right=31.1, bottom=29.0, top=71.1]
-heatmap(temperature_europe, c=:magma, aspectratio=1, frame=:box)
+heatmap(temperature_europe, c=:cividis, aspectratio=1, frame=:box)
 ```
 
 The next step will be to coarsen these data, which requires to give the number
@@ -60,7 +60,7 @@ temperature_europe_coarse = coarsen(temperature_europe, Statistics.mean, (2, 2))
 One again, we can plot these data:
 
 ```@example temp
-heatmap(temperature_europe_coarse, aspectratio=1, c=:magma, frame=:box)
+heatmap(temperature_europe_coarse, aspectratio=1, c=:cividis, frame=:box)
 ```
 
 Finally, we can compare our different clipping and approximations to the overall

docs/src/man/data.md

@@ -1,8 +1,23 @@
-# Bioclimatic data
+# Datasets
 
-The package offers access to bioclimatic datasets - they are downloaded, saved
-to the disk, and then read locally.
+The package offers access to bioclimatic and other datasets - they are
+downloaded, saved to the disk, and then read locally. Please note that some of
+them require a lot of memory, so make sure your machine can handle them.
+
+## Worldclim 2.1
 
 ```@docs
 worldclim
 ```
+
+## CHELSA V1
+
+```@docs
+bioclim
+```
+
+## EarthEnv landcover
+
+```@docs
+landcover
+```

docs/src/man/operations.md

@@ -2,6 +2,7 @@
 
 ```@docs
 coarsen
+slidingwindow
 latitudes
 longitudes
 clip

docs/src/man/overloads.md

@@ -22,6 +22,12 @@ Base.maximum
 Base.minimum
 ```
 
+## From `Broadcast`
+
+```@docs
+broadcast
+```
+
 ## From `Statistics`
 
 ```@docs