library(rvest)
library(xml2)
library(sf)
library(geojsonio)
url = "http://www.landxml.org/webapps/landxmlsamples.aspx"

url_file = paste0("http://landxml.org/schema/LandXML-2.0/",
                  "samples/Carlson%20Software/nunosurf.xml")
xml_file = file.path("~/Downloads", basename(url_file))
if (!file.exists(xml_file)) {
  download.file(url_file, xml_file)
}

res = read_xml(xml_file) %>%
  xml2::xml_ns_strip()
all_names = res %>% xml_children() %>% xml_name

surfs = res %>%
  xml_nodes(xpath = "//Surfaces")
mat_table = res %>%
  xml_nodes(xpath = "//MaterialTable") %>%
  xml_children() %>%
  xml_attrs()
mat_table = do.call(rbind, mat_table)
materials = res %>%
  xml_nodes(xpath = "//MaterialTable") %>%
  xml_children() %>%
  xml_text()
txt_table = res %>%
  xml_nodes(xpath = "//TextureImageTable") %>%
  xml_children() %>%
  xml_attrs()
txt_table = do.call(rbind, txt_table)
textures = res %>%
  xml_nodes(xpath = "//TextureImageTable") %>%
  xml_children() %>%
  xml_text()

units = res %>%
  xml_nodes(xpath = "//Units")
coord = res %>%
  xml_nodes(xpath = "//CoordinateSystem")
xml_attrs(coord)[[1]]
app = res %>%
  xml_nodes(xpath = "//Application")
app = xml_attrs(app)[[1]]
if (length(surfs) == 1) {
  surfs = surfs[[1]]
}
surfs = as_list(surfs)
surf_name = attributes(surfs$Surface)$name
surf_type = attributes(surfs$Surface$Definition)$surfType

stopifnot(all(names(surfs$Surface$SourceData) == "Boundaries"))
bounds = surfs$Surface$SourceData$Boundaries
bounds = lapply(bounds, function(x) x$PntList3D)
bounds = lapply(bounds, function(x) {
  if (length(x) != 1) {
    stop(paste0("Boundaries are more complex then ",
                "one long vector, failing"))
  }
  strsplit(trimws(x), " ")[[1]]
})


#################################
# Point and face data
#################################
stopifnot(all(names(surfs$Surface$Definition) %in% c("Pnts", "Faces")))
stopifnot(length(names(surfs$Surface$Definition)) == 2)
get_values = function(pts, check_name = "P") {
  stopifnot(all(names(pts) == check_name))
  pt_ids = sapply(pts, function(x) {
    a = attributes(x)
    stopifnot(all(names(a) == "id"))
    a$id
  })
  pt_ids = unname(pt_ids)
  pt_ids = unlist(pt_ids)

  pts = sapply(pts, function(x) {
    if (length(x) != 1) {
      stop(paste0("Points are more complex then ",
                  "one long vector, failing"))
    }
    x = unlist(x)
    if (length(x) != 1) {
      stop(paste0("Points are more complex then ",
                  "one long vector, failing"))
    }
    strsplit(trimws(x), " ")[[1]]
  })
  pts = unname(pts)
  stopifnot(length(pts) == ncol(pt_ids))
  L = list(values = pts)
  L$ids = pt_ids
  L
}
points = get_values(surfs$Surface$Definition$Pnts, "P")
faces = get_values(surfs$Surface$Definition$Faces, "F")

class(faces$values) = "numeric"
class(points$values) = "numeric"
faces = faces$values
points = points$values

faces = cbind(t(faces), 1)
points = t(points)
out = tmesh3d(vertices = points, indices = faces)