introduce rgbplot

Author: Max Freudenberg

Project.toml

@@ -18,6 +18,7 @@ Flatten = "4c728ea3-d9ee-5c9a-9642-b6f7d7dc04fa"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 Missings = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
@@ -28,6 +29,7 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 Adapt = "2, 3.0"

src/Rasters.jl

@@ -20,14 +20,16 @@ import Adapt,
        Flatten,
        GeoInterface,
        HDF5,
+       ImageCore,
        PolygonInbounds,
        ProgressMeter,
        Missings,
        Mmap,
        NCDatasets,
        RecipesBase,
        Reexport,
-       Setfield
+       Setfield,
+       Statistics
 
 Reexport.@reexport using DimensionalData, GeoFormatTypes, RasterDataSources
 
@@ -40,7 +42,8 @@ using DimensionalData: Name, NoName
 using .Dimensions: StandardIndices, DimTuple
 using .LookupArrays: LookupArrayTuple 
 
-using RecipesBase: @recipe, @series
+using RecipesBase: @recipe, @series, @userplot
+using Statistics: quantile
 using Base: tail, @propagate_inbounds
 
 using Setfield: @set, @set!

src/plotrecipes.jl

@@ -193,6 +193,58 @@ end
     end
 end
 
+# plot multispectral images as rgb
+@userplot RGBPlot
+
+@recipe function f(p::RGBPlot; bands=1:3, low=0.02, high=0.98)
+    r = first(p.args)
+    if !(typeof(r) <: AbstractRaster)
+        error("Plotting RGB images is only implemented for AbstractRasters.")
+    end
+    if !hasdim(r, Band)
+        error("Raster must have a band dimension.")
+    end
+    if !(0 <= low <= 1) || !(0 <= high <= 1)
+        error("'low' and 'high' have to be between 0 and 1 (inclusive).")
+    end
+    img = Float32.(copy(r[Band([bands...])]))
+    normalize!(img, low, high)
+    img = permutedims(img, (Band, X, Y))
+    img = DD.reorder(img, X=>DD.ForwardOrdered, Y=>DD.ForwardOrdered)
+    plot_image = ImageCore.colorview(RGB, img)
+
+    yguide, xguide = DD.label(dims(img, (Y,X)))
+
+    y, x = map(Rasters._prepare, dims(img, (Y,X)))
+
+    rdt = DD.refdims_title(img; issingle=true)
+    :title --> (rdt === "" ? Rasters._maybename(img) : Rasters._maybename(img) * "\n" * rdt)
+    :xguide --> xguide
+    :xrotation --> 45
+    :yguide --> yguide
+    :tick_direction --> :out
+    :framestyle --> :box
+
+    if all(d -> lookup(d) isa Mapped, (x, y))
+        :xlims --> mappedbounds(x)
+        :ylims --> mappedbounds(y)
+        :aspect_ratio --> :equal
+    else
+        :xlims --> bounds(img, x)
+        :ylims --> bounds(img, y)
+        bnds = bounds(img, (X, Y))
+        # # TODO: Rethink this....
+        s1, s2 = map(((l, u),) -> (u - l), bnds) ./ (size(img, X), size(img, Y))
+        square_pixels = s2 / s1
+        :aspect_ratio --> square_pixels
+    end
+
+    :seriestype --> :image
+    :yflip --> false
+    DD.index(x), DD.index(y), plot_image'
+end
+
+
 # Plots.jl heatmaps pixels are centered.
 # So we should center the index, and use the projected value.
 _prepare(d::Dimension) = d |> _maybe_shift |> _maybe_mapped
@@ -248,3 +300,26 @@ function DD.refdims_title(refdim::Band; issingle=false)
     end
 end
 
+function eachband_view(r::Raster)
+    nbands = size(r, Band)
+    return (view(r, Band(n)) for n in 1:nbands)
+end
+
+function normalize!(raster, low=0.1, high=0.9)
+    if !hasdim(raster, Band)
+        l = quantile(skipmissing(raster), low)
+        h = quantile(skipmissing(raster), high)
+        raster .-= l
+        raster ./= h - l + eps(float(eltype(raster)))
+        raster .= clamp.(raster, zero(eltype(raster)), one(eltype(raster)))
+    else
+        for band in eachband_view(raster)
+            l = quantile(skipmissing(band), low)
+            h = quantile(skipmissing(band), high)
+            band .-= l
+            band ./= h - l + eps(float(eltype(raster)))
+            band .= clamp.(band, zero(eltype(raster)), one(eltype(raster)))
+        end
+    end
+    return raster
+end

src/skipmissing.jl

@@ -35,10 +35,13 @@ end
 
 _missing(x, itr) = ismissing(x) || x == missingval(itr)  # mind the order, as comparison with missing returns missing
 function _missing(x::AbstractFloat, itr)
-    if isnan(missingval(itr))
-        return ismissing(x) || isnan(x)
+    # x is an AbstractFloat here and hence cannot be nothing or missing
+    if isnothing(missingval(itr))
+        return false
+    elseif isnan(missingval(itr))
+        return isnan(x)
     else
-        return ismissing(x) || x == missingval(itr)
+        return x == missingval(itr)
     end
 end