Plot a georeferenced image as image

[mauro3]

Is it possible to plot a georeferenced image as an image? This is a RGB image.

Running

using GeoData, Plots, ArchGDAL
fl = download("https://data.geo.admin.ch/ch.swisstopo.swissimage-dop10/swissimage-dop10_2020_2671-1161/swissimage-dop10_2020_2671-1161_2_2056.tif", "ortho.tif")
a = geoarray(fl)
plot(a)

gives

but I expect

(the image extents are only approximate)

PyPlot.imshow(a[1:2:end,1:2:end,:]) works (minus the missing coordinates of course).

Here the type of a:

julia> typeof(a)
GDALarray{UInt8, 3, String, Tuple{X{LinRange{Float64}, Projected{Ordered{ForwardIndex, ForwardArray, ForwardRelation}, Regular{Float64}, Intervals{Start}, WellKnownText{GeoFormatTypes.CRS, String}, Nothing}, Metadata{:GDAL, Dict{Any, Any}}}, Y{LinRange{Float64}, Projected{Ordered{ReverseIndex, ReverseArray, ForwardRelation}, Regular{Float64}, Intervals{Start}, WellKnownText{GeoFormatTypes.CRS, String}, Nothing}, Metadata{:GDAL, Dict{Any, Any}}}, Band{UnitRange{Int64}, Categorical{Ordered{ForwardIndex, ForwardArray, ForwardRelation}}, NoMetadata}}, Tuple{}, Symbol, Metadata{:GDAL, Dict{Symbol, Any}}, Nothing, Tuple{Int64, Int64, Int64}}

[rafaqz]

Sure. I think ArchGDAL has image handling now too. The question is how to set/detect that a particular tiff can be shown as an image and how to handle plotting it.

1. a command to load a file as an image?
2. a command/plot recipe to plot a GeoArray as an image? is there seriestype=:image or something in Plots.jl?
3. autodetect and add some trait to the Band dimension? Like the mode is not Categorical but Image?
   Then you could also do A = set(A, Band=RGBimage()) or something, to set it if not detected so normal plot would be an image.

I'm not fully across band coloring in geotiff, but we should probably try to accommodate the standard

[mauro3]

I don't know either what makes a GeoTiff an image. Sounds like 2 is the easiest? 3 seems to be overkill for starters. 1 sounds ok too.

[rafaqz]

Ok we can use some flag like plot(A; series=:image) in the plot recipes, that does seem easy. If you want to look into how bands are usually shown in RGB in GeoTiff and similar we can do it roughly how it's normally done?

[rafaqz]

R raster has plotRGB that does pretty much what you want here:

plotRGB(b, r=1, g=2, b=3)

https://rspatial.org/terra/pkg/6-plotting.html

It's cool that you can set which layer gets which color. But we don't want a plots dep, so better to add a keyword with a NamedTuple

plot(A; rgb_bands=(r=1, g=2, b=3))

Then we just need a conditional check for the rgb_bands (or a different keyword) in our recipe, and combine bands if we find it. @mauro3 does that sound reasonable?

[maxfreu]

Here's what I use to plot. I think it can easily be fleshed out a bit and turned into a recipe. As far as I know there is no image function in Plots, see here. Makie has a bit nicer image plotting, where you can even flip through the images in a series using a Slider, but it is also slower than Plots.

using ImageCore
using Statistics
using Plots
using Rasters
using DimensionalData
const DD = DimensionalData

function eachband(r::Raster)
    bands = dims(r, Band)
    return (view(r, Band(b)) for b in bands)
end

function normalize!(raster, low=0.1, high=0.9)
    for band in eachband(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
    return raster
end

function plot_raster(r::Raster; bands=[1,2,3], low=0.02, high=0.98)
    img = float32.(copy(r[Band([bands...])]))
    normalize!(img, low, high)
    img = permutedims(img, (Band, X, Y))
    img = DimensionalData.reorder(img, DD.ForwardOrdered)
    x = DD.index(reorder(dims(img, X), DD.ForwardOrdered))
    y = DD.index(reorder(dims(img, Y), DD.ForwardOrdered))
    plottable_img = colorview(RGB, parent(img))
    Plots.plot(x,y,plottable_img,
               title = string(name(r)),
               xlabel = label(dims(r, X)),
               ylabel = label(dims(r, Y)),
               )
end

[rafaqz]

If you want to write it up that would be awesome.

We just need a keyword to pass to plot to trigger it - preferably an existing one. Is there still an :image series type in Plots.jl?

And/or imagebands for you bands keyword?

[maxfreu]

I'll write it up when I find the time. Is it possible to add plot recipes for Plots and Makie at the same time without depending on them?

[rafaqz]

Plots yes, just Recipes.jl is enough. See the current plot recipes.jl file.

Would be good to have Makie.jl recipes working too, but I havent had time yet.

[maxfreu]

Would be good to have Makie.jl recipes working too, but I havent had time yet.

Yes, but this will introduce a new dependency on MakieCore, correct?

[rafaqz]

Totally. This use case is a reason that got separated out ;)

MakieOrg/Makie.jl#996