georeferenced output

[potrykus]

tmap produces PDFs wih clean vector data. Shouldn't tmap automatically georeference that output - i.e. make its output geo-enabled? Here is a python overview: https://gis.stackexchange.com/questions/49646/is-it-possible-to-georeference-an-existing-un-georeferenced-pdf. This is what I found for the R raster package: https://rpubs.com/Rubio-Polania/1123497.

[Nowosad]

Related discussion: https://github.com/r-tmap/tmap/issues/383

[potrykus]

Thanks. Avenza is an app for mobile use of geo-referenced PDFs (including in a disconnected back-country hiking use case.) They are an end-user platform for US Bureau of Land Management final products https://www.blm.gov/maps/georeferenced-PDFs (expand "How to download and use"). Right now those two major players might be georeferencing maps manually, and may be locked into an old ArcGIS workflow: Avenza-ready Google Earth imagery georeferenced PDF in QGIS. (QGIS is a GIS system that cannot compete with ArcGIS for feature-completeness.)

Given this technological state, tmap georeferenced output would be a powerful and potentially widely-consumed application of the R sf + tmap architecture.

[mtennekes]

Thx for this input @potrykus

Just tested, and it is already possible in sf:

library(sf)
library(tmap)

data(World)

st_write(World, "World.pdf", driver = "PDF")

So now the question is how replace the sf plotting mechanism with tmap. Any ideas @edzer or @tim-salabim ?

If someone could help me with the gdal part, I am happy to implement the other part (passing on the map coordinates of the reference objects).

[edzer]

So now the question is how replace the sf plotting mechanism with tmap.

st_write() only writes data, it doesn't plot: there are no visual elements passed on to it (colors, symbols, line properties etc). If you want to create a GeoPDF from a tmap plot, you'd have to add the geo elements to R's pdf driver - I'd guess these are offset & scale for x and y axes, and a coordinate reference system.

[mtennekes]

Can I use st_write for that purpose? E.g. like this:

tm = tm_shape(World) + tm_poygons()
tmap_save(tm, "World.pdf")

#within tmap_save:
pdf(file)
print(tm)
dev.off()

sf::st_write(?, "World.pdf", offset = , scale = , crs = )

Is so, do you know via which of its arguments these can be passed on (..., dataset_options, layer_options, ... ?

[edzer]

I think st_write() does offset, scale and CRS by itself, as it is an OGR (GDAL) driver; it however takes an sf object, not a tmap object.

[tim-salabim]

I cannot add much here, as mapshot does only simple pdf rendering via png screenshots of the html. The only thing I found is https://gdal.org/en/latest/drivers/vector/pdf.html#feature-style-support

[edzer]

Oh, yes, and there's the raster driver too: https://gdal.org/en/latest/drivers/raster/pdf.html#raster-pdf

[tim-salabim]

Oh, yes, and there's the raster driver too: https://gdal.org/en/latest/drivers/raster/pdf.html#raster-pdf

Yes, but if I understand that correctly, outpu can either be RGB(A) or not. In case it's not I guess it'll be greyscale?

[mdsumner]

would have to get the "plotted" extents of the map for the geotransform, i.e, if it was rendered and written out to file with base graphics the a_ullr for GDAL VRT would be par('usr')[c(1, 4, 2, 3)] .

Tools will complain about coordinates out of range for longlat and potentially other crs (the margin of a plot in Mollweide won't have a valid inverse for example, but I don't think that's a huge problem. I'll explore where the render happens and if the overall frame extent can be gotten. (I wrote something for this elsewhere but I can't find it).

Just to clarify, we're talking about raster graphics yeah? In GDAL terms you need the geotransform, which is a rejig of extent and dimension (shape) mixed together as offset and scale, i.e. equivalent of

f <- function(x, dimension) {
    px <- c(diff(x[1:2])/dimension[1L], -diff(x[3:4])/dimension[2L])
 c(xmin = x[1], xres = px[1], yskew = 0, ymax = x[4], xskew = 0, yres = px[2])   
}
f(par("usr"), dev.size("px"))

but, extent is enough reordered into 'a_ullr' with VRT or vrt::// or write to .wld file even.

[mdsumner]

oops sorry, not "usr", it needs the full expanded plot range - in base might need to use 'plt' to get that (which I think I did before somewhere)

at any rate we need the grid viewport extent but handy to have some more examples I expect

[mdsumner]

grconvertX/Y is your friend for base graphics:

device_extent <- function() {
    c(grconvertX(c(0, 1) , "ndc", "user"), 
      grconvertY(c(0, 1), "ndc", "user"))
}

pdf(tf <- tempfile(fileext = ".pdf"))
maps::map(col = "grey")
## get this information relative to the current plot
dext <- device_extent()
dev.off()
print(dext)

##  (equivalent to a call to gdal_translate {tf} out.vrt -a_ullr ...) since GDAL 3.7
dsn <- sprintf("vrt://%s?a_srs=EPSG:4326&a_ullr=%s", tf, paste0(dext[c(1, 4, 2, 3)], collapse = ","))

terra::plot(terra::rast(dsn))
maps::map(add = TRUE, col = "red")

[mdsumner]

is it possible to get this information? Surely in the pipeline from a tmap object to a rendered graphic we could get information about the frame within the device coordinates? We need this for getting decent images from web servers in #975 and #502 also. With base graphics you can do

ideal_dim <- function(x, ...) {
  if (!missing(x)) plot(x, ...)
  dcx <- grconvertX(c(0, 1), "npc", "device")
  dcy <- grconvertY(c(0, 1), "npc", "device")
  ## "ideal" dimension then is
  ceiling(c(diff(dcx), diff(dcy[2:1])))
}

to get the right dimensions to fit within an existing plot, the xsize and ysize to provide to GDAL (using non-nearest resampling is also important, for maps with text or other features on them).

I see that tm_rgb() and the tm_lines() and other vector plotting families set up different framing to each other, and to base graphics. I don't know how to get these equivalent values, par("usr") would be enough. I'll explore grid package and see if it's somehow possible to interrogate.

[mtennekes]

Thx @mdsumner for following up on this.

I use a grid layout of about 20ish rows and columns, and more if facets are used. Some of these grid cell viewports contain the actual maps with the native coordinates. Which ones exactly is stored internally.

What I can do is to look up this information and pass this on to tmap_grob. Perhaps as attribute? The tmap object itself is not enough, because the rendering depends on the aspect ratio and scale of the map (https://r-tmap.github.io/tmap/articles/adv_margins).

[mdsumner]

I realized I can do this now with the exploring I did, it's the same question - what is the right resolution and extent of data for the warper in a facet, is the same information as what is the extent of the entire rendered canvas. I'll get back to it soonish

[mdsumner]

here's a very rough example, I don't know how to export the graphic via code in a way that is faithful to what I see on the current device, so I do that manually - we could simply write a .wld file with the transform or write a relevant VRT - I also use VRT to expand the RGB because I don't know otherwise how to have terra plot the thing again fathfully

library(tmap)
data(land, World)
tm_shape(land) +
     tm_raster("cover")

## here I don't know how to export the current graphic as-is relative to the current device, 
## so I use the UI to "export->save->PNG" my dev.size("px") is 500x500

tree <- grid::grid.ls(viewports = TRUE, flatten = TRUE, grobs = FALSE)
tree$name
vpname <- grep("vp_map", tree$name, value = TRUE)[1]  ## how to get the right vp ... but this is the one we want
grid::seekViewport(vpname)
vp <- grid::current.viewport()
## not sure how to reset our seek viewport ...


## pix is where we are in the current device
pix <- lapply(grid::deviceLoc(grid::unit(c(0, 1), "npc"), grid::unit(c(0, 1), "npc"), device = T), as.numeric)
dm <- unlist(lapply(pix, \(.x) abs(ceiling(diff(.x)))))
ex <- c(vp$xscale, vp$yscale)

## the size of the pixels
res <- diff(ex)[c(1, 3)]/dm

px <- dev.size("px")

## so overall extent of current device is
ovex <- c(ex[1] - pix$x[1] * res[1], ex[2] + (px[1] - pix$x[2]) * res[1], 
  ex[3] - (px[2] - pix$y[1]) * res[2], ex[4] + (pix$y[2]) * res[2])

on my setup this 'ovex' is c(-187.5,  187.5, -270.2,  104.8 )


r <- rast("vrt://afile.png?expand=rgb")
r <- flip(r)
ext(r) <- ext(ovex)
plotRGB(r)
maps::map(add = T)

see attached afile.png

[mdsumner]

afile.png

[mdsumner]

bundling this up a little more, we get color table with raster but actual RGB with the polygons so there's still a lot to worry about, also I have no idea how to relate the rendered PDF to this (also because of my viewport traversing I don't leave the plot in a correct state that tmap did)

library(tmap)
data(land, World)
tm_shape(land) +
     tm_raster("cover")

render_tmap <- function(filename = tempfile(fileext = ".png")) {
  px <- dev.size("px")
  
  ## here I don't know how to export the current graphic as-is relative to the current device, 
  ## so I use the UI to "export->save->PNG" my dev.size("px") is 500x500
  rstudioapi::savePlotAsImage(filename, width = px[1], height = px[2])

  tree <- grid::grid.ls(viewports = TRUE, flatten = TRUE, grobs = FALSE)
  vpname <- grep("vp_map", tree$name, value = TRUE)[1]  ## how to get the right vp ... but this is the one we want
  grid::seekViewport(vpname)
  vp <- grid::current.viewport()
  ## not sure how to reset our seek viewport ...


  ## pix is where we are in the current device
  pix <- lapply(grid::deviceLoc(grid::unit(c(0, 1), "npc"), grid::unit(c(0, 1), "npc"), device = T), as.numeric)
  dm <- unlist(lapply(pix, \(.x) abs(ceiling(diff(.x)))))
  ex <- c(vp$xscale, vp$yscale)

  ## the size of the pixels 
 res <- diff(ex)[c(1, 3)]/dm


  ## so overall extent of current device is
  ovex <- c(ex[1] - pix$x[1] * res[1], ex[2] + (px[1] - pix$x[2]) * res[1], 
   ex[3] - (px[2] - pix$y[1]) * res[2], ex[4] + (pix$y[2]) * res[2])

  ## in world file coordinates
  res[2] <- -res[2]
  world <- c(xres = res[1L], yskew = 0, xskew = 0, yres = res[2L], 
    xmin = ovex[1L] + res[1L]/2, ymax = ovex[4L] + res[2L]/2)
  worldfile <- gsub("png$", "wld", filename)
  writeLines(as.character(world), worldfile)
  filename
}


## we need this for the Africa example
#data(World)
#Africa = World[World$continent == "Africa", ]
#tm_shape(Africa) + tm_polygons()
#v <- render_tmap()
#plot(rast(sprintf("vrt://%s?expand=rgb", v)))


## this for the raster example
v <- render_tmap()
plotRGB(rast(v))
maps::map(add = T)

So, there's a rough prototype but not that promising

[mtennekes]

Getting back to this. There are two parts of this issue:

1. How to retrieve the bounding box and corresponding viewport of the map?
2. How to save this information as georeference?

I will take care of 1 (should be pretty straightforward).

However, I lack knowledge about 2: e.g. what georeference file types are there? And how is the reference data encoded/written? As additional file (like wld file)? How about georeferenced pdf? Can you share some of your insights @mdsumner ?

[mdsumner]

Just bbox and crs is enough, any GDAL interface will instantiate a dataset or write a file with that info. I only mention sidecar files as a hack. (If you can fold in 1 I'll do 2) 🙏 Just need to be clear on whether it's the the device or an inner window of that, and obviously some tools will complain if you georeferenced an image that had a margin outside a geographic range (but I don't think that's a problem).

[mtennekes]

Thx 🙏

Part 1 done. Now tmap_save silently returns those metadata:

library(tmap)
tm = tm_shape(World) + tm_polygons("HPI")

geo_ref1 = tmap_save(tm, "map1.png", width = 2000, height = 1400)
#> [tip] Consider a suitable map projection, e.g. by adding `+ tm_crs("auto")`.
#> Map saved to map1.png
#> 
#> Resolution: 2000 by 1400 pixels
#> 
#> Size: 6.666667 by 4.666667 inches (300 dpi)
#> 
#> This message is displayed once per session.
str(geo_ref1)
#> List of 7
#>  $ filename   : chr "map1.png"
#>  $ size_pixels: Named num [1:2] 2000 1400
#>   ..- attr(*, "names")= chr [1:2] "width" "height"
#>  $ dpi        : num 300
#>  $ size_inch  : Named num [1:2] 6.67 4.67
#>   ..- attr(*, "names")= chr [1:2] "width" "height"
#>  $ crs        :List of 2
#>   ..$ input: chr "EPSG:4326"
#>   ..$ wkt  : chr "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic Sys"| __truncated__
#>   ..- attr(*, "class")= chr "crs"
#>  $ bbx        : 'bbox' Named num [1:4] -180 -90 180 83.6
#>   ..- attr(*, "names")= chr [1:4] "xmin" "ymin" "xmax" "ymax"
#>   ..- attr(*, "crs")=List of 2
#>   .. ..$ input: chr "EPSG:4326"
#>   .. ..$ wkt  : chr "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic Sys"| __truncated__
#>   .. ..- attr(*, "class")= chr "crs"
#>  $ rel_coords : Named num [1:4] 0.0848 0.3955 0.9152 0.9685
#>   ..- attr(*, "names")= chr [1:4] "xmin" "ymin" "xmin" "xmax"

geo_ref2 = tmap_save(tm, "map2.png", width = 1400, height = 2000)
#> Map saved to map2.png
#> Resolution: 1400 by 2000 pixels
#> Size: 4.666667 by 6.666667 inches (300 dpi)
str(geo_ref2)
#> List of 7
#>  $ filename   : chr "map2.png"
#>  $ size_pixels: Named num [1:2] 1400 2000
#>   ..- attr(*, "names")= chr [1:2] "width" "height"
#>  $ dpi        : num 300
#>  $ size_inch  : Named num [1:2] 4.67 6.67
#>   ..- attr(*, "names")= chr [1:2] "width" "height"
#>  $ crs        :List of 2
#>   ..$ input: chr "EPSG:4326"
#>   ..$ wkt  : chr "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic Sys"| __truncated__
#>   ..- attr(*, "class")= chr "crs"
#>  $ bbx        : 'bbox' Named num [1:4] -180 -90 180 83.6
#>   ..- attr(*, "names")= chr [1:4] "xmin" "ymin" "xmax" "ymax"
#>   ..- attr(*, "crs")=List of 2
#>   .. ..$ input: chr "EPSG:4326"
#>   .. ..$ wkt  : chr "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic Sys"| __truncated__
#>   .. ..- attr(*, "class")= chr "crs"
#>  $ rel_coords : Named num [1:4] 0.0385 0.4713 0.9615 0.7835
#>   ..- attr(*, "names")= chr [1:4] "xmin" "ymin" "xmin" "xmax"

tm2 = tm_shape(NLD_muni) + tm_polygons("income_high")

geo_ref3 = tmap_save(tm2, "map3.png", width = 2000, height = 1400)
#> Map saved to map3.png
#> Resolution: 2000 by 1400 pixels
#> Size: 6.666667 by 4.666667 inches (300 dpi)
str(geo_ref4)
#> Error: object 'geo_ref4' not found

geo_ref4 = tmap_save(tm2, "map4.png", width = 1400, height = 2000)
#> Map saved to map4.png
#> Resolution: 1400 by 2000 pixels
#> Size: 4.666667 by 6.666667 inches (300 dpi)
str(geo_ref4)
#> List of 7
#>  $ filename   : chr "map4.png"
#>  $ size_pixels: Named num [1:2] 1400 2000
#>   ..- attr(*, "names")= chr [1:2] "width" "height"
#>  $ dpi        : num 300
#>  $ size_inch  : Named num [1:2] 4.67 6.67
#>   ..- attr(*, "names")= chr [1:2] "width" "height"
#>  $ crs        :List of 2
#>   ..$ input: chr "Amersfoort / RD New"
#>   ..$ wkt  : chr "PROJCRS[\"Amersfoort / RD New\",\n    BASEGEOGCRS[\"Amersfoort\",\n        DATUM[\"Amersfoort\",\n            E"| __truncated__
#>   ..- attr(*, "class")= chr "crs"
#>  $ bbx        : 'bbox' Named num [1:4] 13565 306846 278026 619232
#>   ..- attr(*, "names")= chr [1:4] "xmin" "ymin" "xmax" "ymax"
#>   ..- attr(*, "crs")=List of 2
#>   .. ..$ input: chr "Amersfoort / RD New"
#>   .. ..$ wkt  : chr "PROJCRS[\"Amersfoort / RD New\",\n    BASEGEOGCRS[\"Amersfoort\",\n        DATUM[\"Amersfoort\",\n            E"| __truncated__
#>   .. ..- attr(*, "class")= chr "crs"
#>  $ rel_coords : Named num [1:4] 0.124 0.346 0.876 0.968
#>   ..- attr(*, "names")= chr [1:4] "xmin" "ymin" "xmin" "xmax"

^{Created on 2025-04-16 with reprex v2.1.1}

[mdsumner]

I have to warm up to this again, maybe in a few days, btw in $rel_coords the names are not meaningful, I think this is xmin, ymin, xmax, ymax in there?

[mtennekes]

I have to warm up to this again, maybe in a few days, btw in $rel_coords the names are not meaningful, I think this is xmin, ymin, xmax, ymax in there?

Yes, exactly. Typos 🙈

The origin (0,0) is bottom left and (1, 1) is top right. This is where the bounding box is placed in the image. I've already subtracted the inner margins, so the bounding box should be valid. In case facets are drawn, it will refer to the first facet.