support mesh3d data

[SymbolixAU]

TODO

• [x] work on any 'mesh', not just 4 points (e.g., triangles)
• [x] 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?
• [x] ~~close polygons~~ don't need this
• [x] colourvalues lists
• [ ] mesh JS to colour vertices (#177) (if possible )
• [x] calculate_bounds error (in js console), due to bbox
• [x] 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.

https://github.com/rforge/rgl/commit/87bcf7445cc7d7a28df7076cbf314128c6b29f07#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 ;)