nara_map.md

## imagery source
dsn <- "<GDAL_WMS><Service name=\"TMS\"><ServerUrl>http://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/${z}/${y}/${x}</ServerUrl></Service><DataWindow><UpperLeftX>-20037508.34</UpperLeftX><UpperLeftY>20037508.34</UpperLeftY><LowerRightX>20037508.34</LowerRightX><LowerRightY>-20037508.34</LowerRightY><TileLevel>17</TileLevel><TileCountX>1</TileCountX><TileCountY>1</TileCountY><YOrigin>top</YOrigin></DataWindow><Projection>EPSG:900913</Projection><BlockSizeX>256</BlockSizeX><BlockSizeY>256</BlockSizeY><BandsCount>3</BandsCount><MaxConnections>10</MaxConnections><Cache /><ZeroBlockHttpCodes>204,404,403</ZeroBlockHttpCodes></GDAL_WMS>"

## return a function to transform coordinates to nara space (for a given extent)
nr_trans <- function(nara, extent = c(0, ncol(nara), 0, nrow(nara))) {
  function(xy) {
    cbind(scales::rescale(xy[,1], from = extent[1:2], to = c(0, ncol(nara))),
          scales::rescale(xy[,2], from = extent[3:4], to = c(nrow(nara), 0)))
    
  }
}

library(vapour)
im <- gdal_raster_nara(dsn, target_dim = dev.size("px"), 
                       bands = c(1, 2, 3, 3), target_crs = "+proj=laea +lon_0=148", 
                       target_ext = c(-1, 1, -1, 1) * .65 * pi * 6378137)
m <- do.call(cbind, maps::map(plot = F)[1:2])
m[m[,1] > 180, ] <- NA
m <- reproj::reproj_xy(m, attr(im, "projection"), 
                       source = "EPSG:4326")

tr <- nr_trans(im[[1]], attr(im, "extent"))
## make segments, explicit in 4 columns
p2s <- function(x) {
  out <- cbind(utils::head(x, -1L),  utils::tail(x, -1))
 colnames(out) <- c("x0", "x1", "y0", "y1") 
 out
}
segs <- p2s(tr(m))
segs <- na.omit(segs) # [seq(1, nrow(segs), by = 5), ])
v <- im[[1]]

library(nara)  ## coolbutuseless/nara
v <- nr_line(v, segs[,1], segs[,2], segs[,3], segs[,4])
#v <- nr_line(v, 1, nrow(v)/2, ncol(v), nrow(v)/2)
#v <- nr_line(v, ncol(v)/2, 1, ncol(v)/2, nrow(v))

ximage::ximage(v, attr(im, "extent"), asp = 1)

for rough graticule something like this will do it

lonlat <- reproj::reproj_xy(vaster::xy_from_cell(dim(v), attr(im, "extent"), 1:prod(dim(v))), 
                            "EPSG:4326", source = attr(im, "projection"))
lon <- matrix(lonlat[,1], nrow(v), byrow = TRUE)
lat <- matrix(lonlat[,2], nrow(v), byrow = TRUE)
cl <- contourLines(seq(1, ncol(lat)), seq(1, nrow(lat)), t(lat))

for (i in seq_along(cl)) {
  segs <- p2s(tr(do.call(cbind, cl[[i]])))
  v <- nr_line(v, segs[,1], segs[,2], segs[,3], segs[,4])
}

is there a bug in vaster?

cl <- contourLines(vaster::x_from_cell(dim(v), attr(im, "extent"), 1:ncol(v)),
                   vaster::y_from_cell(dim(v), attr(im, "extent"), nrow(v):1),
                    t(lat))