## 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)
Created
June 6, 2024 20:45
-
-
Save mdsumner/8b267893bbd6b2df454baac65129a896 to your computer and use it in GitHub Desktop.
Author
mdsumner
commented
Jun 6, 2024
orthographic, with some more segment clean up for infinite
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>"
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(sds::wms_arcgis_mapserver_ESRI.WorldImagery_tms(), target_dim = dev.size("px"),
bands = c(1, 2, 3, 3), target_crs = "+proj=ortho +lon_0=148",
target_ext = c(-1, 1, -1, 1) * .32 * 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 <- [seq(1, nrow(segs), by = 5), ])
segs <- segs[! rowSums(is.na(segs) | !is.finite(segs)) > 0, ]
v <- im[[1]]
library(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))
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment