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: #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")