support mesh3d data

[SymbolixAU]

TODO

• work on any 'mesh', not just 4 points (e.g., triangles)
• quadmesh uses vb and ib, triangle uses vb and it (?) - use the geometry_columns to specify these names
• triangle mesh isn't showing in Rstudio Browser?
• close polygons don't need this
• colourvalues lists
• mesh JS to colour vertices (Colour path segments / polygon vertices  #177) (if possible )
• calculate_bounds error (in js console), due to bbox
• mesh3d objects
• quadmesh objects
• maintain material$col attribute on mesh3d objects
• (if no holes) use SolidPolygon layer (and therefore binary data format)

[mdsumner]

I had a go at this, converting the quad-type mesh3d into the form shape$data uses in add_polygons.

I'm confused about whether we need a feature for each geometry or not, and I'm confused about conversion to actual json class, so I'll simply share the R-way I'd try to go about this ... perhaps this helps unpacking the mesh from my side.

## convert a 3-column matrix to text embedded triples in json arrays
coord_to_json <- function(x) paste0(unlist(lapply(split(t(x), rep(seq_len(nrow(x)), each = ncol(x))), function(cds) sprintf("[%f,%f,%f]", cds[1], cds[2], cds[3]))), collapse = ",")
## add that first coordinate to the end
coord1 <- function(x) rbind(x, x[1, ])
format_mesh <- function(x, ...) {
	UseMethod("format_mesh")
}
#' format mesh3d into JSON suitable for mapdeck ... WIP 
format_mesh.mesh3d <- function(x, ...) {
	## switch on x$primitivetype == "quad"/"triangle"
	stopifnot(x$primitivetype == "quad")
	feature_template <- '[{"type":"Feature","properties":{"elevation":0,"fill_colour":"#440154FF","stroke_colour":"#440154FF"},%s}]'
	polygon_template <- '"geometry":{"geometry":{"type":"Polygon","coordinates":[[%s]]}}'
	## assume the expanded coords for each primitive
	coords <- sprintf(polygon_template, unlist(lapply(seq_len(ncol(qm$ib)), function(iquad)  coord_to_json(coord1(t(qm$vb[1:3, qm$ib[, iquad] ]))))))
	paste0(sprintf(feature_template, coords, collapse = ","), collapse = ",")
}


library(quadmesh)
library(raster)
## dummy raster
rr <- setExtent(raster::raster(matrix(sample(1:12), 3)), raster::extent(0, 4, 0, 3))
qm <- quadmesh(rr * 10000)  ## something suitably exaggerated

## confused about conversion to 'json' class ...
format_mesh(qm)

Also, it will definitely save the effort of conversion to sf polygons, but it's essentially in the same form here - 5 coordinates for each quad, rather than a vertex array and primitive-index array - but perhaps this will help us find other ways.

[SymbolixAU]

There are a couple of ways I could implement this, and either here or in spatialwidget. Either way, somewhere the quadmesh object needs to be reshaped or parsed into the pseudo-GeoJSON.

1. I could leave spatialwidget as it is and make mapdeck reshape the quadmesh data into the right shape to just be accepted into spatialwidget.

2. Or I could write a custom spatialwidget::api::create_geojson() c++ function to handle quadmesh (as I've done in the example.

I think 1) is more maintainable and easiest to implement, but possibly could incur a slight penalty because of the reshaping. However I think this is the approach I'm going to take.

@mdsumner for what it's worth, I've got a skeleton quadmesh -> sf cpp converter here on branch issue94. It's by no means complete, and I don't need all the sf construction steps, but I made it anyway to make a start. If you think it's useful for a qm_to_sf() function inside quadmesh feel free to take it.

p <- mapdeck:::rcpp_mesh_geojson( qm )
> p
Geometry set for 12 features 
geometry type:  POLYGON
dimension:      XYZM
bbox:           xmin: 0 ymin: 0 xmax: 0 ymax: 0
epsg (SRID):    NA
proj4string:    
First 5 geometries:
POLYGON ZM ((0 3 -24970 1, 1 3 15029.99 1, 1 2 ...
POLYGON ZM ((1 3 15029.99 1, 2 3 17527.49 1, 2 ...
POLYGON ZM ((2 3 17527.49 1, 3 3 35002.5 1, 3 2...
POLYGON ZM ((3 3 35002.5 1, 4 3 104965 1, 4 2 1...
POLYGON ZM ((0 2 -4966.67 1, 1 2 55021.66 1, 1 ...

[mdsumner]

Oh cool thanks! (We definitely need your sf-constructor API, that's a great example for me to learn from - I've found following you elsewhere here has been quite tough. )

Definitely this should be for mesh3d and not quadmesh in particular, the only thing to think about there is that mesh3d can hold either quad or triangle or both types of primitive, but very many other 3d types will extend mesh3d or be trivially converted.

(quadmesh only adds the projection crs and the original raster spec for limited round-tripping)

[SymbolixAU]

prototype example (I may use add_polygon() in stead, rather than add_mesh()). The avg_z fill colour is also a place-holder until I get colourvalues to support lists

library(quadmesh)
library(raster)
## dummy raster
rr <- setExtent(raster::raster(matrix(sample(1:12), 3)), raster::extent(0, 4, 0, 3))
qm <- quadmesh(rr * 10000)  ## something suitably exaggerated


library(mapdeck)
set_token( read.dcf("~/.googleAPI", fields = "MAPBOX"))

mapdeck() %>%
	add_mesh(
		data = qm
		, fill_colour = "avg_z"
                 , vertex = "vb"
                 , index = "ib"
	)

[SymbolixAU]

And using your hawaii exmaple

library(quadmesh)
data("hawaii", package = "marmap")
qm <- hawaii %>% 
	marmap::as.raster() %>% 
	raster::aggregate(fact = 2) %>% 
	quadmesh()

mapdeck() %>%
	add_mesh(
		data = qm
		, fill_colour = "avg_z"
                 , vertex = "vb"
                 , index = "ib"
	)

[SymbolixAU]

Triangle example

tri2 <- qm
tri2$it <- triangulate_quads(tri2$ib)
tri2$ib <- NULL
tri2$primitivetype <- "triangle"

mapdeck() %>%
	add_mesh(
		data = tri2
		, fill_colour = "avg_z"
		, vertex = "vb"
		, index = "it"
	)

[SymbolixAU]

@mdsumner is there a way to know which of the names attributes on a quadmesh object refer to the vertices and indices?

I'm thinking something akin to the way sf objects store an attribute sf_column

sf <- spatialwidget::widget_melbourne

attr(sf, "sf_column")
"geometry"

where this lets us know the column name of the geometries.

Can we infer anything about these attributes to tell us the vb object are the vertices, and the ib object are the indices?

attributes( qm )
$names
[1] "vb"              "ib"              "primitivetype"   "material"        "normals"         "texcoords"      
[7] "raster_metadata" "crs"            

$class
[1] "quadmesh" "mesh3d"   "shape3d" 

[mdsumner]

I could certainly add that to quadmesh, but mesh3d will be more widespread - to me the value of quadmesh()/triangmesh() is that it creates these rgl types from spatial data - though maybe I don't understand?

mesh3d will always have vb, and either it (triangle indices) or ib (quad indices), and sometimes both, so isn't the mesh3d class, primitivetype, and presence of it/ib enough?

[SymbolixAU]

mesh3d will always have vb, and either it (triangle indices) or ib (quad indices), and sometimes both, so isn't the mesh3d class, primitivetype, and presence of it/ib enough?

yeah if this is a consistent rule with these objects then it's enough for me to work with.

I'm just not familiar enough with these objects to know how they're built / stored or what they represent. But I'm getting there...

[mdsumner]

Oh I see, there's a bit in the quadmesh vignette - but I agree, it's extremely hard to understand these types, that's an acquired assumption on my part - I kind of went very quickly from total confusion to near-mastery. I'll try to distill a clear description! Especially material properties, I recently learnt a lot more about those

[mdsumner]

The mesh3d manual page has a pithy description: https://rdrr.io/rforge/rgl/man/mesh3d.html It lacks how the material properties are stored though it does describe the helper arg meshColor that controls interpretation, so I'll focus a bit on that.

[SymbolixAU]

so cool!

( I tried using the whole of etopo but the data didn't exist )

library(quadmesh)
qm <- quadmesh::quadmesh( quadmesh::etopo )
qm$vb[3, ] <- qm$vb[3, ] * 100

mapdeck() %>%
	add_mesh(
		data = qm
	)

mapdeck() %>%
	add_mesh(
		data = qm
		, palette = "magma"
	)

[mdsumner]

haha, I'll get you an entire coverage!

[mdsumner]

Here is GEBCO 2019 reduced 100x, saved as compressed .rds 1.1Mb ( and zipped so I could attach it here)

GEBCO_lowres.zip

Attachment GEBCO_lowres.zip derived from GEBCO 2019 Grid, provided by GEBCO Compilation Group (2019) GEBCO 2019 Grid (doi:10.5285/836f016a-33be-6ddc-e053-6c86abc0788e) https://www.gebco.net/data_and_products/gridded_bathymetry_data/gebco_2019/gebco_2019_info.html

[SymbolixAU]

sweet

[mdsumner]

This is cool! I've had a look at generalizing for any mesh3d object - so you can leave that to me if you like, I'll also write up how the colour properties are applied.

[dcooley]

hold off on the colour properties; I'm going to change it once colourvalues accepts lists. That way I can use every 'z' value, not just the average. Then every vertex will get its own colour and every polygon will have a gradient fill.

I'll happily take a PR on the generalised mesh3d objects though

[mdsumner]

cool cool, coming soon - I got as far as seeing the C++ is specific to quad types, so I'll push what I have and create some example mesh3d data sets

[mdsumner]

Oh, I see it does already deal with vb/it at the C++ - I'm just not seeing anything - still exploring.

[SymbolixAU]

You might be experiencing

triangle mesh isn't showing in Rstudio Browser?

If you open it in a browser do you get the plot?

[mdsumner]

It does work with triangle type, so I think I'm doing something wrong when creating the mesh3d- so I'll retreat for now and make sure my examples are proper!

[mdsumner]

Ah! Indeed, plus I was mixing up Melbourne and Hobart there ....

I find with triangles it doesn't focus the view though does with quad - more than enough for me to go on now, thanks!

[mdsumner]

Here I get some elevation data from Mapbox (no bathymetry is available), after some updates to quadmesh and reproj on github. I'm using mds-mesh3d which just now is up to date with master here.

There's two comment/uncomment options for a Tas or US region, and for mesh as triangles or quads.

## remotes::install_github(c("hypertidy/ceramic", "hypertidy/quadmesh"))
## remotes::install_github("mdsumner/mapdeck@mds-mesh3d")
library(quadmesh)
library(ceramic)
library(raster)
library(mapdeck)
## choose Tas/Vic
im <- cc_elevation(extent(140, 150, -45, -36), zoom = 3); im <- im * 70
## or choose Yosemite
#im <- cc_elevation(extent(-120, -119, 37.5, 38.5), zoom = 6) * 5
tm <- triangmesh(im)
qm <- quadmesh(im)
## project mesh to longlat (choose quad or triangle mesh here)
mesh <- reproj(tm, 4326, source = projection(im))
#mesh <- reproj(qm, 4326, source = projection(im))

tok <- Sys.getenv("MAPBOX_API_KEY")
mapdeck(token = tok) %>%
	add_mesh(
		data = mesh
	)

Context on changes in reproj/quadmesh:

• reproj needs only reproject the vertices of the mesh, and now assumes longlat if there's no crs or class, and the data seems sensible
• quadmesh has a reproj method but it now results in the class being removed (back to mesh3d) and there being no raster_metadata (because the extent now refers to the old coordinates)

(Together this gives a lot more ability to get data into mesh3d form. )

[SymbolixAU]

Thanks also for your write-up on rgl. I now see how the colours are used, and therefore probably should replicate the behaviour if a $colour attribute is supplied.

[mdsumner]

IMPORTANT ANNOUNCEMENT:

the primitivetype property was removed completely and you should not expect it to be a part of mesh3d.

r-forge/rgl@87bcf74#diff-fb1a2ca62de6cfcab71eef95a8103b38

I'll have to clean up quadmesh a fair bit to remove it, you're expected to look only for it or ib arrays for the type/s of primitives available (mesh3d can have both, and there might be other types but I'm not sure).

I'll be back with more about the constructors: tmesh3d(), qmesh3d() which are clearly now the ways you are supposed to build these objects (not from templates like I always have before).

[SymbolixAU]

@mdsumner do mesh objects ever have holes in them?

[mdsumner]

Definitely!

library(quadmesh)
v <- volcano
v[v > 165] <- NA
qm <- quadmesh(v)
rgl::shade3d(qm, col = "grey"); rgl::rglwidget()

"But are they really holes?" I have the mighty boosh in my ear

You can arbitrarily remove triangles or vertices, it doesn't matter - though you end up with a mess of course if not careful.

rgl::clear3d()
qm$vb[1, qm$vb[1, ] > 30] <- NA
rgl::shade3d(qm, col = "grey"); rgl::rglwidget()

[SymbolixAU]

But are they really holes?

cool, this is what I was hoping for; they are 'missing' quads/triangles, rather than 'holes' in the sense of sf POLYGONs

[mdsumner]

That's right, I don't have CRAN code to turn sf polygons into mesh3d but hypertidy/anglr can do it

• one of the magic mesh things is you can have the triangles that define the holes, but simply consider them "invisible" - that's kind of how polygon triangulation occurs, the triangles get culled out or flagged

[mdsumner]

Are you looking for a tidy-ish encoding much like the sf_to_df()? I wonder if you are seeking some kind of entity definition that makes sense. (If so, I agree! mesh3d is cool but is more a rendering format than a data set, anglr::DEL and silicate::TRI are first-cut attempts at something sensible, but I feel the need to do something better). If not, nvm and carry on ;)

[mdsumner]

Oh I see it's about SolidPolygon, sounds promising!

I know you must be a bit busy rn with CRAN (plus general crazy world), but I added a mesh_plot.mesh3d method in quadmesh - so you can do a 2D sanity-check of any mesh you are working with:

library(quadmesh)  ## need latest remotes::install_github("hypertidy/quadmesh")
v <- volcano
v[v > 165] <- NA
qm <- quadmesh(v)
## make it a regular mesh3d (this is the conversion world I need to flesh out)
class(qm) <- c("mesh3d", "shape3d")
qm$raster_metadata <- qm$crs <- NULL

mesh_plot(qm)

## triangles too
tm <- triangmesh(v)
## regular mesh3d
class(tm) <- c("mesh3d", "shape3d")
tm$raster_metadata <- qm$crs <- NULL

mesh_plot(tm)

^{Created on 2020-03-19 by the reprex package (v0.3.0)}

[dcooley]

Are you looking for a tidy-ish encoding much like the sf_to_df()?

Oh I see it's about SolidPolygon, sounds promising!

Yes and yes. If I can get all coordinates in a column then I can use the SolidPolygonLayer, which accepts binary data (deck.gl's definition), so will plot a lot faster.

[mdsumner]

I don't totally get the column thing, but is this what it needs? A SolidPolygon is the explicit path around the coordinates?

Fourth column in xyzi classifies each path:

#x <- quadmesh:::as_mesh3d_matrix(matrix(1:2, 2))
#dput(x)

x <- structure(list(vb = structure(c(0, 1, 1, 1, 1, 1, 1.5, 1, 2, 
                                1, 2, 1, 0, 0, 1, 1, 1, 0, 1.5, 1, 2, 0, 2, 1), 
        .Dim = c(4L, 6L)), 
        ib = structure(c(1L, 2L, 5L, 4L, 2L, 3L, 6L, 5L), .Dim = c(4L,2L)), 
        material = list(color = c("#440154FF", "#218F8DFF", "#FDE725FF",
                                  "#440154FF", "#218F8DFF", "#FDE725FF")), 
        normals = NULL, texcoords = NULL, 
               meshColor = "vertices"), class = c("mesh3d", "shape3d"))


## handle it or ib
id <- x$ib

xyzi <- cbind(t(x$vb[1:3, id]), rep(seq_len(ncol(id)), each = nrow(id)))
xyzi
     [,1] [,2] [,3] [,4]
[1,]    0    1  1.0    1
[2,]    1    1  1.5    1
[3,]    1    0  1.5    1
[4,]    0    0  1.0    1
[5,]    1    1  1.5    2
[6,]    2    1  2.0    2
[7,]    2    0  2.0    2
[8,]    1    0  1.5    2

[dcooley]

In deck.gl terminology, a SolidPolygonLayer is a polygon without holes. A polygon is just one javascript-array of coordinates.

And to make use of their binary data format I need to flatten the coordinates into a single vector

[x,y,x,y,x,y,x,x,y,...]

Which is trivial if I have a data.frame in R with x, y, (z) columns, which I now have working for the POINT plotting layers scatterplot, column and pointcloud.

So if the mesh object can be made into a single column directly then this will save some conversion time - mesh --> df --> vector.

---

But, it may also be useful to have a similar mesh_to_df() function somewhere anyway to do that initial mesh --> df step

[mdsumner]

So iiuc each i-th SolidPolygon is

## i is individually a column-index of x$ib or x$it
as.vector(x$vb[1:2, x$ib[,i]])  # quad

as.vector(x$vb[1:2, x$it[,i]])  # triangle

and an intermediate df with all of them is

idx <- x$ib ## or x$it
data.frame(x = x$vb[1, idx], y = x$vb[2, idx], z = x$vb[3, idx], solid_polygon_id = rep(seq_len(ncol(idx)), each = nrow(idx)))

To close each polygon you would first update the index with idx <- rbind(idx, idx[1, ])

Or am I still misunderstanding?

[dcooley]

yes I believe your understanding is correct.

[mdsumner]

Thoughts while they're fresh, I just used this code successfully:

library(mapdeck)
mapdeck() %>% add_mesh(data = melbourne_mesh)
library(ceramic)
cc <- cc_elevation(raster::extent(147.5, 147.6, -42.1, -42.0), 
                   zoom = 9)

mesh <- reproj::reproj(quadmesh::triangmesh(cc), 4326)
mesh$vb[3, ] <- mesh$vb[3, ] * 35
mesh$material$col <- colourvalues::color_values(mesh$vb[3, mesh$it])
mapdeck() %>% add_mesh(data = mesh)

Caveats: having troubles with raster and ceramic, so I'm dev versions of those - the code above may not work. The plan is for hypertidy/anglr to have a comprehensive as.mesh3d() suite, the question of reprojection is a bit open since mesh3d does not store crs - but quadmesh::triangmesh does - for now I think raw coordinate reproj+replacement is safest. I'd like to move away from quadmesh for this.

• the material property on mesh3d is color (not col, we probably should be flexible to this though)
• primitivetype is no longer part of the spec (mesh3d can have it and ib, so presence is the trigger)
• there's a meshColor property that dictates whether colours apply to vertices (and are continuous) or to faces (discrete)

I'll build a set of mesh3d examples so we can use that as a standard, something we can fix if they are specified wrongly.

The reason for this note is that mapdeck's expectations for a mesh3d, and the code I have to create them (in quadmesh) are not aligned, and not consistent with the modern mesh3d spec. So it always takes me ages to get things sorted and it's time to share an easy-example pool.

[mdsumner]

Here's a few examples of quad and triangle mesh3d:

https://github.com/hypertidy/anglr/tree/issue89/data-raw/mapdeck

Don't worry about the code so much, that's in flux and lots of hypertidy-cleanout - there's a saved list of meshes, in strictly modern and proper rgl form (I hope).