Create surface forms from matrix or 'raster' data for flexible plotting and conversion to other mesh types. The functions 'quadmesh' or 'triangmesh' produce a continuous surface as a 'mesh3d' object as used by the 'rgl' package. This is used for plotting raster data in 3D (optionally with texture), and allows the application of a map projection without data loss and many processing applications that are restricted by inflexible regular grid rasters. There are discrete forms of these continuous surfaces available with 'dquadmesh' and 'dtriangmesh' functions.
A quadmesh is a dense mesh describing a topologically continuous surface of 4-corner primitives. This is also known as a cell-based raster but in those contexts the corner coordinates and the continuous nature of the mesh is completely implicit. By making the dense mesh explicit we have access to every corner coordinate (not just the centres) which allows for some extra facilities over raster grids.
This package provides helpers for working with this mesh interpretation of gridded data to enable
You can install:
Quadmesh provides a number of key features that are not readily available to more commonly used geospatial applications.
In raster there is an implicit “area-based” interpretion of the extent
and value of every cell. Coordinates are implicit, and centre-based but
extent reflects a finite and constant width and height.
library(quadmesh)library(raster)#> Loading required package: spr <- raster(matrix(1:12, 3), xmn = 0, xmx = 4, ymn = 0, ymx = 3)qm <- quadmesh(r)op <- par(xpd = NA)plot(extent(r) + 0.5, type = "n", axes = FALSE, xlab = "", ylab = "")plot(r, col = rep(c("palevioletred3", "aliceblue"), 6), add = TRUE)text(coordinates(r), lab = seq_len(ncell(r)))plot(extent(r), add = TRUE)## with quadmesh there are 20 corner coordinatespoints(t(qm$vb[1:2, ]), pch = "x", cex = 2)
Every individual quad is described implicitly by index into the unique
corner coordinates. This format is built upon the
mesh3d class of the
## coordinates, transpose hereqm$vb[1:2, ]#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]#> x 0 1 2 3 4 0 1 2 3 4 0 1 2#> y 3 3 3 3 3 2 2 2 2 2 1 1 1#> [,14] [,15] [,16] [,17] [,18] [,19] [,20]#> x 3 4 0 1 2 3 4#> y 1 1 0 0 0 0 0## indexes, also transposeqm$ib#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]#> [1,] 1 2 3 4 6 7 8 9 11 12 13 14#> [2,] 2 3 4 5 7 8 9 10 12 13 14 15#> [3,] 7 8 9 10 12 13 14 15 17 18 19 20#> [4,] 6 7 8 9 11 12 13 14 16 17 18 19
This separation of geometry (the 20 unique corner coordinates) and topology (12 sets of 4-index) is the key concept of a mesh and is found in many domains that involve computer graphics and modelling.
A quadmesh is centre-based, and each raster pixel occupies a little rectangle, but the centre-based interpretation is better described by a mesh of triangles.
Every centre point of the grid data is a node in this mesh, while the inside corners with their ambiguous value from 4 neighbouring cells are not explicit. Note that a triangulation can be top-left/bottom-right aligned as here, or bottom-left/top-right aligned, or be a mixture.
op <- par(xpd = NA)plot(extent(r) + 0.5, type = "n", axes = FALSE, xlab = "", ylab = "")plot(r, col = rep(c("palevioletred3", "aliceblue"), 6), add = TRUE)tri <- triangmesh(r)points(t(tri$vb[1:2, ]))polygon(t(tri$vb[1:2, head(as.vector(rbind(tri$it, NA)), -1)]))
Quads are trivially converted into triangle form by splitting each in two. Note that this is different again from the centre-based triangle interpretation.
tri2 <- qmtri2$it <- triangulate_quads(tri2$ib)tri2$ib <- NULLtri2$primtivetype <- "triangle"op <- par(xpd = NA)plot(extent(r) + 0.5, type = "n", axes = FALSE, xlab = "", ylab = "")plot(r, col = rep(c("palevioletred3", "aliceblue"), 6), add = TRUE)points(t(tri2$vb[1:2, ]))polygon(t(tri2$vb[1:2, head(as.vector(rbind(tri2$it, NA)), -1)]), lwd = 2, lty = 2)
etopo data set is used to create a plot in a local map
projection. Here each cell is drawn by reprojecting it directly and
individually into this new coordinate system.
library(quadmesh)library(raster)## VicGridprj <- "+proj=lcc +lat_1=-36 +lat_2=-38 +lat_0=-37 +lon_0=145 +x_0=2500000 +y_0=2500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"er <- crop(etopo, extent(110, 160, -50, -20))system.time(mesh_plot(er, crs = prj))#> user system elapsed#> 0.78 0.03 0.82## This is faster to plot and uses much less data that converting explicitly to polygons.library(sf)#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3p <- st_transform(spex::polygonize(er), prj)plot(st_geometry(p), add = TRUE)
system.time(plot(p, border = NA))
#> user system elapsed #> 0.37 0.07 0.48 pryr::object_size(er) #> 37.7 kB pryr::object_size(p) #> 1.68 MB pryr::object_size(quadmesh(er)) #> 169 kB
The data size and timing benefits are more substantial for larger data sets.
We get exactly what we asked for from
mesh_plot, without the complete
remodelling of the data required by grid resampling.
pol <- "+proj=stere +lat_0=-90 +lon_0=147 +datum=WGS84"mesh_plot(etopo, crs = pol)
plot(projectRaster(etopo, crs = pol), col = viridis::viridis(64))
triangmesh types are extensions of the
mesh3d, and so are readily used by that package’s high level functions
class(qm)#>  "quadmesh" "mesh3d" "shape3d"class(tri)#>  "triangmesh" "mesh3d" "shape3d"
spex package has functions
qm_rasterToPolygons_sp() which provide very fast conversion of raster
grids to a polygon layer with 5 explicit coordinates for every cell.
sf now provide an even faster version by using GDAL).
rr <- disaggregate(r, fact = 20)system.time(spex::polygonize(rr))#> user system elapsed#> 0.10 0.00 0.11system.time(raster::rasterToPolygons(rr))#> user system elapsed#> 1.73 0.00 1.74## stars has now improved on spex by calling out to GDAL to do the worksystem.time(sf::st_as_sf(stars::st_as_stars(rr), merge = FALSE, as_points = FALSE))#> user system elapsed#> 0.19 0.02 0.21
Using a triangulation version of a raster grid we can build an index of weightings for a new of of arbitrary coordinates to estimate the implicit value at each point as if there were continuous interpolation across each primitive.
WIP - see
Have feedback/ideas? Please let me know via issues.
Many aspects of this package have developed in conjunction with the angstroms package for dealing with ROMS model output.
Please note that this project is released with a Contributor Code of Conduct. By participating in this project you agree to abide by its terms.
Now using the
reproj package for coordinate transformations.
mesh_plot() method for
TRI objects from silicate/anglr.
dquadmesh to provide the discrete interpretation of a grid.
bary_index to provide a basis for interpolation of a regular grid
dtriangmesh to provide the triangle-interpretation of a raster (bound to centre coordinates).
quadmesh() is now generic with a new method for 'matrix'.
Some fixes and additional features for
mesh_plot() to allow plotting with
qm_as_raster to round-trip to a RasterLayer.
z argument of
quadmesh() may now have a different coordinate system to the
x. The mesh coordinates of
x are silently reprojected if needed and
used to extract values in the native coordinate system of
qsc function to return a very simple six quad rendition of the Quadrilateralized Spherical Cube
mesh_plot function gains true curvilinear grid support by allowing input of
coords, a two-layer
raster of longitude and latitude values for the cell geography. This aligns with the approach used in the
quadmesh function gains a new argument 'texture_filename' to control the output PNG file.
quadmesh functions gains a
texture argument, to map on a
3-layer RGB raster (in 0-255) via a temporary PNG file.
mesh_plot for basic vectorized raster plot with optional near near-lossless
Index arrays are now created with integer mode.
triangulate_quads function to convert quad index to pairs of triangle index.
Added supporting information to package.
llh2xyz is now exported to transform from longitude latitude coordinates to
geocentric xyz coordinates.