Skip to content

Commit

Permalink
fix merge
Browse files Browse the repository at this point in the history
  • Loading branch information
mdsumner committed Apr 19, 2023
2 parents 8cd4eb6 + dfd9433 commit 0220581
Show file tree
Hide file tree
Showing 51 changed files with 226 additions and 4,357 deletions.
5 changes: 1 addition & 4 deletions .Rbuildignore
Original file line number Diff line number Diff line change
@@ -1,12 +1,9 @@
^codecov\.yml$
^appveyor\.yml$
^\.travis\.yml$
^data-raw$
^CODE_OF_CONDUCT\.md$
^README\.Rmd$
^ceramic\.Rproj$
^\.Rproj\.user$
^docs$
^gdalwmscache$
^\.Rhistory$
^cran-comments\.md$
^CRAN-RELEASE$
Expand Down
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
.Rproj.user
.Rhistory
gdalwmscache
3 changes: 2 additions & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,8 @@ Imports:
reproj (>= 0.4.2),
utils,
crsmeta,
terra
terra,
vapour
Suggests:
covr,
spelling,
Expand Down
3 changes: 3 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,8 @@

# ceramic 0.7.0

* Begin move to use GDAL for the read, separate tile downloading from raster input.

* Fixed `cc_kingston` location.

* Removed all references to non-supported mapbox styles, we now have 'mapbox.satellite' and and 'mapbox.terrain-rgb'. Others can be used with a custom styles URL from your own mapbox account.
Expand Down
180 changes: 180 additions & 0 deletions R/gdal.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,180 @@
loc_extent <- function(x, buffer, dimension = NULL) {
if (missing(x) && missing(buffer)) {
## get the whole world at zoom provided, as a neat default
x <- cbind(0, 0)
buffer <- c(20037508, 20037508)

}
bbox_pair <- spatial_bbox(x, buffer)
ext <- as.vector(bbox_pair$extent)

## widest
dx_dy <- diff(ext)[c(1, 3)]
if (is.null(dimension)) {
dimension <- rep(min(dev.size("px")), 2L)
was_one <- TRUE
} else {
was_one <- length(dimension) == 1
}
# print(dimension)
dimension <- rep(dimension, length.out = 2L)
if (was_one) dimension[which.min(dx_dy)] <- 0L
c(ext, dimension)
}


format_out <- function (x)
{
dimension <- x$dimension
projection <- x$projection
x <- x$extent
if (inherits(x, "SpatRaster")) {
if (!requireNamespace("terra"))
stop("terra package required but not available, please install it")
x <- terra_out(x)
if (is.na(x$projection) && x$lonlat)
x$projection <- "OGC:CRS84"
x$type <- "terra"
}
if (inherits(x, "BasicRaster")) {
if (!requireNamespace("raster"))
stop("raster package required but not available, please install it")
x <- raster_out(x)
if (is.na(x$projection) && x$lonlat)
x$projection <- "OGC:CRS84"
x$type <- "raster"
}
if (inherits(x, "numeric")) {
if (is.null(projection))
projection <- "OGC:CRS84"
if (is.null(dimension))
dimension <- as.integer(256 * sort(c(1, diff(x[1:2])/diff(x[3:4])),
decreasing = FALSE))
lonlat <- grepl("lonlat", projection) || grepl("4326",
projection) || grepl("4269", projection) || grepl("OGC:CRS84",
projection) || grepl("GEOGCS", projection)
x <- list(extent = x, dimension = dimension, projection = projection,
lonlat = lonlat, type = "matrix")
}
x$lonlat <- vapour::vapour_crs_is_lonlat(x$projection)
x
}

terra_out <- function (x)
{
if (isS4(x) && inherits(x, "SpatRaster")) {
x <- try(list(extent = c(terra::xmin(x), terra::xmax(x),
terra::ymin(x), terra::ymax(x)), dimension = dim(x)[2:1],
projection = terra::crs(x), lonlat = terra::is.lonlat(x),
terra = TRUE), silent = TRUE)
if (inherits(x, "try-error"))
stop("cannot use terra grid")
}
x
}
gdal_mapbox <- function (extent = c(-180, 180, -90, 90), ..., dimension = NULL,
projection = "OGC:CRS84", resample = "near", source = NULL)
{
xraster <- extent
x <- format_out(list(extent = extent, dimension = dimension,
projection = projection))
src1 <- "<GDAL_WMS><Service name=\"TMS\"><ServerUrl>https://api.mapbox.com/v4/mapbox.satellite/${z}/${x}/${y}.jpg?access_token=%s</ServerUrl></Service><DataWindow><UpperLeftX>-20037508.34</UpperLeftX><UpperLeftY>20037508.34</UpperLeftY><LowerRightX>20037508.34</LowerRightX><LowerRightY>-20037508.34</LowerRightY><TileLevel>22</TileLevel><TileCountX>1</TileCountX><TileCountY>1</TileCountY><YOrigin>top</YOrigin></DataWindow><Projection>EPSG:3857</Projection><BlockSizeX>256</BlockSizeX><BlockSizeY>256</BlockSizeY><BandsCount>3</BandsCount>"
src2 <- sprintf("<UserAgent>%s</UserAgent><Cache /></GDAL_WMS>",
getOption("HTTPUserAgent"))
src <- paste0(src1, src2)


if (is.null(source)) {
rso <- sprintf(src, get_api_key())
}
else {
rso <- source
}
if (is.na(x$projection)) {
message("no projection specified, calling warper without a target projection: results not guaranteed")
x$projection <- ""
}
info <- vapour::vapour_raster_info(rso[1])
bands <- 1:3
if (info$bands < 3) {
bands <- 1L
}
vals <- vapour::gdal_raster_image(rso, target_ext = x$extent,
target_dim = x$dimension, target_crs = x$projection,
resample = resample, ..., bands = bands)

x$type <- "terra"
xraster <- terra::rast(terra::ext(attr(vals, "extent")),
ncol = attr(vals, "dimension")[1L],
nrow = attr(vals, "dimension")[2L], crs = "EPSG:3857", vals = vals[[1]])
xraster
}




gdal_aws <- function (extent = c(-180, 180, -90, 90), ..., dimension = NULL,
projection = "OGC:CRS84", resample = "near", source = NULL)
{
xraster <- extent
x <- format_out(list(extent = extent, dimension = dimension,
projection = projection))

src <- sprintf("<GDAL_WMS>\n <Service name=\"TMS\">\n <ServerUrl>https://s3.amazonaws.com/elevation-tiles-prod/geotiff/${z}/${x}/${y}.tif</ServerUrl>\n </Service>\n <DataWindow>\n <UpperLeftX>-20037508.340000</UpperLeftX>\n <UpperLeftY>20037508.340000</UpperLeftY>\n <LowerRightX>20037508.340000</LowerRightX>\n <LowerRightY>-20037508.340000</LowerRightY>\n <TileLevel>15</TileLevel>\n <TileCountX>1</TileCountX>\n <TileCountY>1</TileCountY>\n <YOrigin>top</YOrigin>\n </DataWindow>\n <Projection>EPSG:3857</Projection>\n <BlockSizeX>512</BlockSizeX>\n <BlockSizeY>512</BlockSizeY>\n <BandsCount>1</BandsCount>\n <UserAgent>%s</UserAgent>\n</GDAL_WMS>",
getOption("HTTPUserAgent"))



if (is.null(source)) {
rso <- src
} else {
rso <- source
}
if (is.na(x$projection)) {
message("no projection specified, calling warper without a target projection: results not guaranteed")
x$projection <- ""
}
vals <- vapour::gdal_raster_data(rso, target_ext = x$extent,
target_dim = x$dimension, target_crs = x$projection,
resample = resample, ..., bands = 1L)


xraster <- terra::rast(terra::ext(attr(vals, "extent")),
ncol = attr(vals, "dimension")[1L],
nrow = attr(vals, "dimension")[2L], crs = "EPSG:3857", vals = vals[[1]])
xraster
}


gdal_terrainrgb <- function (extent = c(-180, 180, -90, 90), ..., dimension = NULL,
projection = "OGC:CRS84", resample = "near", source = NULL)
{
xraster <- extent
x <- format_out(list(extent = extent, dimension = dimension,
projection = projection))

src1 <- "<GDAL_WMS><Service name=\"TMS\"><ServerUrl>https://api.mapbox.com/v4/mapbox.terrain-rgb/${z}/${x}/${y}.png?access_token=%s</ServerUrl></Service><DataWindow><UpperLeftX>-20037508.34</UpperLeftX><UpperLeftY>20037508.34</UpperLeftY><LowerRightX>20037508.34</LowerRightX><LowerRightY>-20037508.34</LowerRightY><TileLevel>22</TileLevel><TileCountX>1</TileCountX><TileCountY>1</TileCountY><YOrigin>top</YOrigin></DataWindow><Projection>EPSG:3857</Projection><BlockSizeX>256</BlockSizeX><BlockSizeY>256</BlockSizeY><BandsCount>3</BandsCount>"
src2 <- sprintf("<UserAgent>%s</UserAgent><Cache /><ZeroBlockHttpCodes>404</ZeroBlockHttpCodes></GDAL_WMS>",
getOption("HTTPUserAgent"))

src <- paste0(src1, src2)
if (is.null(source)) {
rso <- sprintf(src, get_api_key())
} else {
rso <- source
}
if (is.na(x$projection)) {
message("no projection specified, calling warper without a target projection: results not guaranteed")
x$projection <- ""
}
vals <- vapour::gdal_raster_data(rso, target_ext = x$extent,
target_dim = x$dimension, target_crs = x$projection,
resample = resample, ..., bands = 1:3)

d <- -10000 + ((vals[[1]] * 256 * 256 + vals[[2]] * 256 + vals[[3]]) * 0.1)

xraster <- terra::rast(terra::ext(attr(vals, "extent")),
ncol = attr(vals, "dimension")[1L],
nrow = attr(vals, "dimension")[2L], crs = "EPSG:3857", vals = d)
xraster
}
79 changes: 36 additions & 43 deletions R/locale.R
Original file line number Diff line number Diff line change
Expand Up @@ -22,13 +22,15 @@
#' `"https://api.mapbox.com/styles/v1/mdsumner/cjs6yn9hu0coo1fqhdqgw3o18/tiles/512/{zoom}/{x}/{y}"`
#'
#' Currently must be considered in-development.
#'
#' @param loc a longitude, latitude pair of coordinates, or a spatial object
#' @param buffer with in metres to extend around the location, ignored if 'loc' is a spatial object with extent
#' @param type character string of provider imagery type (see Details)
#' @param ... arguments passed to internal function, specifically `base_url` (see Details)
#' @param zoom desired zoom for tiles, use with caution - if `NULL` is chosen automatically
#' @param max_tiles maximum number of tiles to be read into memory - if `NULL` is set by zoom constraints
#' @param debug optionally print out files that will be used
#' @param zoom deprecated (use `dimension`)
#' @param max_tiles deprecated
#' @param debug deprecated
#' @param dimension one or two numbers, used to determine the number of pixels width, height - set one to zero to let GDAL figure it out, or leave as NULL to get something suitable
#'
#' @return A [raster::brick()] object, either 'RasterBrick' with three layers (Red, Green, Blue) or with
#' a single layer in the case of [cc_elevation()].
Expand All @@ -51,82 +53,73 @@
#' plotRGB(im)
#' }
cc_location <- function(loc = NULL, buffer = 5000,
type = "mapbox.satellite", ..., zoom = NULL, max_tiles = NULL, debug = FALSE) {
if (!is.null(zoom) && !is.null(max_tiles)) stop("'zoom' and 'max_tiles' cannot be both set, one must be NULL")
if (is.null(zoom) && is.null(max_tiles)) max_tiles <- 16L
locdata <- get_tiles(x = loc, buffer = buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug)
if (debug) {
return(invisible(NULL))
}
suppressWarnings(make_raster(locdata))
type = "mapbox.satellite", ..., zoom = NULL, max_tiles = NULL, debug = FALSE, dimension = NULL) {
if (!is.null(zoom) || !is.null(max_tiles)) message("'zoom' and 'max_tiles' are ignored")
#if (is.null(zoom) && is.null(max_tiles)) max_tiles <- 16L
#locdata <- get_tiles(x = loc, buffer = buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug)

locdata <- loc_extent(loc, buffer, dimension)

#if (debug) {
# return(invisible(NULL))
#}
d <- gdal_mapbox(extent = locdata[1:4], dimension = as.integer(locdata[5:6]), projection = "EPSG:3857")
d
}
#' @name cc_location
#' @export
cc_macquarie <- function(loc = c(158.93835,-54.49871), buffer = 5000,
type = "mapbox.satellite", ..., zoom = NULL, max_tiles = NULL, debug = FALSE) {
cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug)
type = "mapbox.satellite", ..., zoom = NULL, max_tiles = NULL, debug = FALSE, dimension = NULL) {
cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug, dimension = dimension)
}

#' @name cc_location
#' @export
cc_davis <- function(loc = c(77 + 58/60 + 3/3600,
-(68 + 34/60 + 36/3600)),
buffer = 5000, type = "mapbox.satellite", ..., zoom = NULL, max_tiles = NULL, debug = FALSE) {
buffer = 5000, type = "mapbox.satellite", ..., zoom = NULL, max_tiles = NULL, debug = FALSE, dimension = NULL) {
# 68 34 36 S 77 58 03 E
cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug)

cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug, dimension = dimension)
}
#' @name cc_location
#' @export
cc_mawson <- function(loc = c(62 + 52/60 + 27/3600,
-(67 + 36/60 + 12/3600)), buffer = 5000, type = "mapbox.satellite", ..., zoom = NULL, max_tiles = NULL, debug = FALSE) {
-(67 + 36/60 + 12/3600)), buffer = 5000, type = "mapbox.satellite", ..., zoom = NULL, max_tiles = NULL, debug = FALSE, dimension = NULL) {
# 67 36 12 S 62 52 27 E

cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug)

cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug, dimension = dimension)
}

#' @name cc_location
#' @export
cc_casey <- function( loc = cbind(110 + 31/60 + 36/3600,
-(66 + 16/60 + 57/3600)), buffer = 5000, type = "mapbox.satellite", ...,zoom = NULL, max_tiles = NULL,debug = FALSE) {
-(66 + 16/60 + 57/3600)), buffer = 5000, type = "mapbox.satellite", ...,zoom = NULL, max_tiles = NULL,debug = FALSE, dimension = NULL) {
#66 16 57 S 110 31 36 E

cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug)

cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug, dimension = dimension)
}
#' @name cc_location
#' @export
cc_heard <- function(loc = c(73 + 30/60 + 30/3600,
-(53 + 0 + 0/3600)), buffer = 5000, type = "mapbox.satellite",...,zoom = NULL, max_tiles = NULL, debug = FALSE) {
-(53 + 0 + 0/3600)), buffer = 5000, type = "mapbox.satellite",...,zoom = NULL, max_tiles = NULL, debug = FALSE, dimension = NULL) {
# 53 S 73 30 E.

cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug)

cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug, dimension = dimension)
}
#' @name cc_location
#' @export
cc_kingston <- function(loc = c(147.2901,
-42.98682), buffer = 5000, type = "mapbox.satellite", ...,zoom = NULL, max_tiles = NULL, debug = FALSE) {
cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug)

-42.98682), buffer = 5000, type = "mapbox.satellite", ...,zoom = NULL, max_tiles = NULL, debug = FALSE, dimension = NULL) {
cc_location(loc, buffer, type = type, ..., zoom = zoom, max_tiles = max_tiles, debug = debug, dimension = dimension)
}

#' @name cc_location
#' @export
cc_elevation <- function(loc = NULL, buffer = 5000, type = NULL, ...,zoom = NULL, max_tiles = NULL, debug = FALSE) {
rgbunpack <- FALSE
cc_elevation <- function(loc = NULL, buffer = 5000, type = NULL, ...,zoom = NULL, max_tiles = NULL, debug = FALSE, dimension = NULL) {
locdata <- loc_extent(loc, buffer, dimension)

if (is.null(type)) {
rgbunpack <- TRUE
type <- "mapbox.terrain-rgb"
}
dat <- cc_location(loc, buffer = buffer, type = type, zoom = zoom, max_tiles = max_tiles, debug = debug, ...)
if (rgbunpack) {
projection(dat) <- NA ## prevent yelling about lost datums that we aren't specifying
dat <- -10000 + ((dat[[1]] * 256 * 256 + dat[[2]] * 256 + dat[[3]]) * 0.1)

raster::crs(dat) <- sp::CRS(.merc(), doCheckCRSArgs = FALSE)
dat <- gdal_terrainrgb(extent = locdata[1:4], dimension = as.integer(locdata[5:6]), projection = "EPSG:3857")

} else {
dat <- gdal_aws(extent = locdata[1:4], dimension = as.integer(locdata[5:6]), projection = "EPSG:3857")
}
dat

}
2 changes: 1 addition & 1 deletion R/spatial.R
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ spatial_bbox <- function(loc, buffer = NULL) {
xmax = bb_points_lonlat[2,1], ymax = bb_points_lonlat[2,2])
user_points <- bb_points

list(tile_bbox = tile_bbox, user_points = user_points)
list(tile_bbox = tile_bbox, user_points = user_points, extent = as.vector(bb_points))
}
spex_to_pt <- function(x) {
pt <- cbind(mean(c(raster::xmax(x), raster::xmin(x))),
Expand Down
Loading

0 comments on commit 0220581

Please sign in to comment.