hypertidy / textures

simple texture plot
GNU General Public License v3.0
12 stars 0 forks source link

texture mapping as a raster warper #3

Open mdsumner opened 4 years ago

mdsumner commented 4 years ago

Raster warping is an expensive and drastic remodelling of a grid to a new grid - and it's lossy. Texture mapping is a cheap and effective and potentially lossless way of projecting an image arbitrarily onto graphical primitives (quads or triangles). If those primitives are fine enough we can do pretty effective image reprojection - dynamically, there's no remodelling or copy it's just the graphics engine figuring out what to colour each pixel within a triangle based on barycentric interpolation.

Usually we don't have a geospatial context, the image is in [0,1,0,1] and we map it to our surface, but this can be done to convert between map projections. (This is a feature of the quadmesh and anglr packages, which allow textures draped onto rasters or polygon surfaces).

Here in textures package I want to explain this weird trick from a less specialist perspective.

This works, but what's with the weird zoom scale (I just guessed until it worked, and it works with different crs)

library(rgl)
library(raster)
library(anglr)
library(textures)

## create a quadmesh and map an image texture to it 
## (no tricks, we are in [0,1,01] like a standard bitmap)
tfile <- tempfile(fileext = ".png")
quad0 <- quad_texture(c(64, 64), texture = tfile)
png::writePNG(ga_topo$img/255, tfile)

## a geospatial coastline , in 2D spatial and 3D form
aus <- subset(anglr::simpleworld, sovereignt == "Australia")
aus_merc <- sp::spTransform(aus, ga_topo$crs)
aus_wire <- copy_down(DEL0(aus_merc), 1000)

## convert the geometry of the mesh to geospatial, these are Mercator coordinates of the image
## as a spatial raster (the texture is already mapped, so this doesn't affect that)
quad0$vb[1L,] <- scales::rescale(quad0$vb[1,], to = ga_topo$extent[c("xmin", "xmax")])
quad0$vb[2L,] <- scales::rescale(quad0$vb[2L,], to = ga_topo$extent[c("ymin", "ymax")])

## a different projection (Lambert Conformal Conic)
prj2 <- "+proj=lcc +lon_0=110 +lat_0=-30 +lat_1=-10 +lat_2=-40"

## reproject Mercator to another projection (this makes the raster "curvy" and would require
## destructive resampling to create a new raster in this projection
quad0$vb[1:2, ] <- t(reproj::reproj(t(quad0$vb[1:2, ]), target = prj2,
                                    source = ga_topo$crs)[,1:2])

## plot the mesh
rgl::plot3d(quad0, specular = "black")
## add the coastline (it's in LCC now, not Mercator)
wire3d(reproj::reproj(aus_wire, prj2), add = T)
set_scene(zoom = 0.485)  ## hack the view-scale (why does this almost work)
aspect3d(1, 1, 0)

## render the scene to a raster brick (3-layer RGB in 0:255)
r <- setExtent(flip(t(brick(rgl.pixels()*255)), "y"), par3d()$bbox[1:4])

## plot the spatial raster brick (it's Lambert Conformal Conic - LCC)
plotRGB(r)
## add the coastline by projecting to LCC (thereby proving that we achieved our goal)
plot(spTransform(aus_merc, prj2), add = T)

This is a geospatial brick in Lambert Conformal Conic (see prj2) rendered from an rgl window, with a proper spatial transform (mercator -> lcc) applied to overplot:

image