From 051df92ef56d26e54899e7f4f6f169bf64b6e1d5 Mon Sep 17 00:00:00 2001 From: Max Freudenberg Date: Tue, 29 Nov 2022 15:02:02 +0100 Subject: [PATCH 1/2] introduce rgbplot --- Project.toml | 2 ++ src/Rasters.jl | 7 +++-- src/plotrecipes.jl | 75 ++++++++++++++++++++++++++++++++++++++++++++++ src/skipmissing.jl | 9 ++++-- 4 files changed, 88 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 8833c40df..9475a591b 100644 --- a/Project.toml +++ b/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" diff --git a/src/Rasters.jl b/src/Rasters.jl index a4d0eb512..ee957ed52 100644 --- a/src/Rasters.jl +++ b/src/Rasters.jl @@ -20,6 +20,7 @@ import Adapt, Flatten, GeoInterface, HDF5, + ImageCore, PolygonInbounds, ProgressMeter, Missings, @@ -27,7 +28,8 @@ import Adapt, 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! diff --git a/src/plotrecipes.jl b/src/plotrecipes.jl index 179f8b9c2..dcc001e9b 100644 --- a/src/plotrecipes.jl +++ b/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 \ No newline at end of file diff --git a/src/skipmissing.jl b/src/skipmissing.jl index 35988bb33..c0700534e 100644 --- a/src/skipmissing.jl +++ b/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 From 53177eb8cfd182f998ef4d51b481c8801748366e Mon Sep 17 00:00:00 2001 From: Max Freudenberg Date: Tue, 29 Nov 2022 15:27:04 +0100 Subject: [PATCH 2/2] specify ImageCore compat --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 9475a591b..29c76cfd7 100644 --- a/Project.toml +++ b/Project.toml @@ -45,6 +45,7 @@ Flatten = "0.4" GeoFormatTypes = "0.4" GeoInterface = "1" HDF5 = "0.14, 0.15, 0.16" +ImageCore = "<1.0" Missings = "0.4, 1" NCDatasets = "0.10, 0.11, 0.12" PolygonInbounds = "0.2.0"