convert_vrt_rct.md

x_from_col <- function(dimension, bbox, col) {
  col[col < 1] <- NA
  col[col > dimension[1L]] <- NA
  xres <- diff(bbox[c(1, 3)]) / dimension[1]
  bbox[1] - xres/2 + col * xres
}
y_from_row <- function(dimension, bbox, row) {
  row[row < 1] <- NA
  row[row > dimension[2]] <- NA
  yres <- diff(bbox[c(2, 4)]) / dimension[2]
  bbox[4] + yres/2 - row * yres
}

##url <- "https://raw.githubusercontent.com/mdsumner/rema-ovr/main/REMA-2m_dem_ovr.vrt"
url <-  "https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt"

library(xml2)
## get requisite information from the VRT about each ComplexSource
xml <- read_xml(url)
dst <- xml |> xml_find_all(".//DstRect")
src <- xml |> xml_find_all(".//SrcRect")
xOff <- as.integer(xml_attr(dst, "xOff"))
yOff <- as.integer(xml_attr(dst, "yOff"))
xSize <- as.integer(dst |> xml_attr("xSize"))
ySize <- as.integer(dst |> xml_attr("ySize"))

## get the overall dataset dim, bbox, crs
wkt <- xml |> xml_find_all("SRS") |> xml_text()

dm <- as.integer(c(xml |> xml_find_all("//VRTDataset") |> xml_attr("rasterXSize"), 
    xml |> xml_find_all("//VRTDataset") |> xml_attr("rasterYSize")))
gt <- as.numeric(unlist(xml |> xml_find_all("//GeoTransform") |> xml_text() |> strsplit(",")))
bbox <- c(gt[1], gt[4] + dm[2] * gt[6], gt[1] + dm[1] * gt[2], gt[4])

 
xmin <- x_from_col(dm, bbox, xOff + 1)
xmax <- x_from_col(dm, bbox, xOff + xSize)
ymax <- y_from_row(dm, bbox, yOff + 1)
ymin <- y_from_row(dm, bbox, yOff + ySize)
  
## all bbox of ComplexSource (or SimpleSource) in VRT:
bb <- cbind(xmin = xmin,  ymin = ymin, xmax = xmax, ymax = ymax)
  
## so, then we can do rct(), which is useful for geos::geos_strtree(rc) an querying that with xy()
rc <- wk::rct(bb[,1], bb[,2], bb[,3], bb[,4])
plot(rc)