Attempt at plotting a raster from Rasters.jl in a different projection.

projected-mesh-image.jl

using Rasters
using GLMakie
using Proj
using Extents
using ImageCore: RGB, N0f8
using GeometryBasics: uv_mesh, Tesselation, FaceView
import GeometryBasics as GB
Band = Rasters.Band

function transform_mesh!(trans, m)
    # mutate the nodes
    for i in eachindex(m.position)
        p = m.position[i]
        m.position[i] = Point2f(trans(p))
    end
    return m
end

# TODO upstream
function trans_bounds(
    trans::Proj.Transformation,
    bbox::Extent,
    densify_pts::Integer = 21,
)::Extent
    xlims, ylims = Proj.bounds(trans, bbox.X, bbox.Y; densify_pts)
    return Extent(X = xlims, Y = ylims)
end

# from https://www.naturalearthdata.com/downloads/50m-raster-data/
uri = "https://github.com/visr/ribasim-artifacts/releases/download/v0.1.0/50-natural-earth-1-downsampled.tif"
# raster = Raster(string("/vsicurl/", uri); lazy=false)
# raster = Raster(string("/vsicurl/", uri); lazy=true)
raster = Raster("50-natural-earth-1-downsampled.tif"; lazy = false)

# TODO like NDVI @bandmath(red=2/1, green=1, blue=3)
rgb = (red = 2, green = 1, blue = 3)

# reinterpret the UInt8 (0-255) raster as N0f8 (0.0-1.0), used in Images
A = reinterpret.(N0f8, raster)
# This also works on lazy Rasters, but viewing them as an image will be very slow
# since it is done one pixel at a time.
img = RGB.(view(A, Band(1)), view(A, Band(2)), view(A, Band(3)))

# TODO add Proj.Transformation constructors for WellKnownText{CRS} etc.?
source_crs = convert(String, crs(img))
trans = Proj.Transformation(source_crs, "+proj=moll", always_xy = true)
makietrans(point) = Point2f(trans(point))

# we need to transform the bounds using Proj to set good limits on the axis
bbox = extent(img)
dest_bbox = trans_bounds(trans, bbox)
limits = (dest_bbox.X..., dest_bbox.Y...)

# make a Rect2f of the world
rect = Rect2f(bbox.X[1], bbox.Y[1], bbox.X[2] - bbox.X[1], bbox.Y[2] - bbox.Y[1])

coarsen = 10
x = lookup(raster, X)
y = lookup(raster, Y)

Rasters.order(raster)
Rasters.order(x)
Rasters.order(y)

gridmesh = uv_mesh(Tesselation(rect, (length(x) รท coarsen, length(y) รท coarsen)))

wireframe(gridmesh)

# with bottom=south, etc:
# the mesh (x,y) starts at south-west going east
# the raster (x,y) starts at north-east going west
# somehow reversing Y and permuting aligns it, but that
# makes the raster start at south-west going north
# can we construct the mesh in the same direction as the raster?
color = reverse(permutedims(img), dims = Y)

# rect = Rect2f(180, 90, -360, -180) # rotr90
# rect = Rect2f(180, -90, -360, 180) # rotr90 flipud
# rect = Rect2f(-180, 90, 360, -180) # rotl90 flipud
# rect = Rect2f(-180, -90, 360, 180)  # rotl90
# wireframe(rect)
ts = Tesselation(rect, (length(x) รท coarsen, length(y) รท coarsen))
gm = uv_mesh(ts)
typeof(gm)
gm.position
gm.uv
getfield(gm, :simplices)
getfield(gm, :simplices)[30]
propertynames(gm)
gm[1]
gm[30]

# Create the mesh ourselves such that it is ordered like the grid.
# This is not supported by uv_mesh(Tesselation(rect))
# Try FaceView constructor. Maybe this can also be lazy?
x

# view with steprange
DimPoints(img)
DimPoints(img)
collect(DimPoints(img))

length(gm)
size(gm)
begin
    fig, ax, p =
        mesh(gm; color = color, shading = false, axis = (; aspect = DataAspect(), limits))
    p.transformation.transform_func[] = Makie.PointTrans{2}(makietrans)
    fig
    display(fig)
end

gm[1]
gm[30]


img
gridmesh[1]
gridmesh[30]
x[1]
y[1]
img[:, 1]
img[:, 100]
img[1, :]
size(parent(img))
img[2]
size(img)
lookup(img)