## see discussion on R-Sig-Geo
## https://stat.ethz.ch/pipermail/r-sig-geo/2018-February/026328.html
## GDAL details here http://www.gdal.org/frmt_wms.html
library(raster)
## Loading required package: sp
yesterday <- Sys.Date() - 1
TileLevel <- 3
xml_specification <-
paste0('<GDAL_WMS>
<Service name="TMS">
<ServerUrl>
https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/MODIS_Aqua_CorrectedReflectance_Bands721/default/',
format(yesterday, "%Y-%m-%d"),
'/250m/${z}/${y}/${x}.jpg</ServerUrl>
</Service>
<DataWindow>
<UpperLeftX>-4194304</UpperLeftX>
<UpperLeftY>4194304</UpperLeftY>
<LowerRightX>4194304</LowerRightX>
<LowerRightY>-4194304</LowerRightY>
<TileLevel>', TileLevel, '</TileLevel>
<TileCountX>2</TileCountX>
<TileCountY>2</TileCountY>
<YOrigin>top</YOrigin>
</DataWindow>
<Projection>EPSG:3413</Projection>
<BlockSizeX>512</BlockSizeX>
<BlockSizeY>512</BlockSizeY>
<BandsCount>3</BandsCount>
</GDAL_WMS>
')
## raster is *lazy* here, we only get metadata
rastersource <- raster(xml_specification)
localproj <- projection(rastersource)
llproj <- "+init=epsg:4326"
## this is not meaningful, because the truly polar extent of
## the source means this wraps to the entire extent
#local_ex <- projectExtent(raster(extent(c(-135, 170), c(52, 75)),
# crs = llproj), localproj)
## we can crop in index space (note that it's a different orientation
## to the usual xmin/xmax/ymin/ymax
## r1, r2, c1, c2, representing the first and last row and column number)
## and note that any subsetting means the "lazy-ness" contract above
## is broken and data is actually read
rc <- c(1, 150, 1, 200)
## this works!!
r <- crop(rastersource, extent(rastersource, rc[1], rc[2], rc[3], rc[4]))
plot(r)
## see in context
plot(extent(rastersource), asp = 1)
plot(extent(rastersource, rc[1], rc[2], rc[3], rc[4]), add = TRUE)
plot(r, add = TRUE)
data("wrld_simpl", package = "maptools")
plot(sp::spTransform(wrld_simpl, projection(rastersource)), add = TRUE)
## use what we learned to drive rgdal more direclty
## so we can extract at a given level-of-detail, not just
## a slice subset
## note (y,x) from raster and rgdal::readGDAL orientation
dm <- c(nrow(rastersource), ncol(rastersource))
subfact <- 16
outdim <- trunc(dm/subfact)
offs <- c(0, 0)
rs <- rgdal::readGDAL(xml_specification, offset = offs, region.dim = dm, output.dim = outdim)
## <GDAL_WMS>
## <Service name="TMS">
## <ServerUrl>
## https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/MODIS_Aqua_CorrectedReflectance_Bands721/default/2018-02-08/250m/${z}/${y}/${x}.jpg</ServerUrl>
## </Service>
## <DataWindow>
## <UpperLeftX>-4194304</UpperLeftX>
## <UpperLeftY>4194304</UpperLeftY>
## <LowerRightX>4194304</LowerRightX>
## <LowerRightY>-4194304</LowerRightY>
## <TileLevel>3</TileLevel>
## <TileCountX>2</TileCountX>
## <TileCountY>2</TileCountY>
## <YOrigin>top</YOrigin>
## </DataWindow>
## <Projection>EPSG:3413</Projection>
## <BlockSizeX>512</BlockSizeX>
## <BlockSizeY>512</BlockSizeY>
## <BandsCount>3</BandsCount>
## </GDAL_WMS>
## has GDAL driver WMS
## and has 8192 rows and 8192 columns
## there we have it, a 512x512 level-of-detail read of the entire
## map
## beware, do not print the rgdal object
## just convert to raster and move on
rs <- raster(rs)
plot(rs, col = viridis::viridis(64))
plot(sp::spTransform(wrld_simpl, projection(rastersource)), add = TRUE)
## compare the lazy source at native resolution ...
rastersource
## class : RasterLayer
## band : 1 (of 3 bands)
## dimensions : 8192, 8192, 67108864 (nrow, ncol, ncell)
## resolution : 1024, 1024 (x, y)
## extent : -4194304, 4194304, -4194304, 4194304 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : <GDAL_WMS>
## <Service name="TMS">
## <ServerUrl>
## https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/MODIS_Aqua_CorrectedReflectance_Bands721/default/2018-02-08/250m/${z}/${y}/${x}.jpg</ServerUrl>
## </Service>
## <DataWindow>
## <UpperLeftX>-4194304</UpperLeftX>
## <UpperLeftY>4194304</UpperLeftY>
## <LowerRightX>4194304</LowerRightX>
## <LowerRightY>-4194304</LowerRightY>
## <TileLevel>3</TileLevel>
## <TileCountX>2</TileCountX>
## <TileCountY>2</TileCountY>
## <YOrigin>top</YOrigin>
## </DataWindow>
## <Projection>EPSG:3413</Projection>
## <BlockSizeX>512</BlockSizeX>
## <BlockSizeY>512</BlockSizeY>
## <BandsCount>3</BandsCount>
## </GDAL_WMS>
## names : GDAL_WMS.
## values : 0, 255 (min, max)
## ... with the subsampled, in-memory, collected object
rs
## class : RasterLayer
## dimensions : 512, 512, 262144 (nrow, ncol, ncell)
## resolution : 16384, 16384 (x, y)
## extent : -4194304, 4194304, -4194304, 4194304 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : band1
## values : 0, 235 (min, max)
## here's an experimental collect() method
#if (!exists("collect") && is.function(collect)) {
collect <- function(x, ...) {
UseMethod("collect")
}
#}
collect.BasicRaster <- function(x, ..., nrows = 512, ncols = nrows) {
dm <- c(raster::nrow(x), raster::ncol(x))
if (nrows != dm[1] || ncols != dm[2]) {
message(sprintf("%s is %ix%i, collecting output of %ix%i \n %s",
deparse(substitute(x)), dm[1], dm[2], nrows, ncols,
"use 'nrows = , ncols = ' to customize output size"))
}
outdim <- c(nrows, ncols)
outdim <- c(min(dm[1], nrows), min(dm[2], ncols))
offs <- c(0, 0)
if (raster::inMemory(x)) {
warning(sprintf("%s is already in-memory ... any resampling will be approximate", deparse(substitute(x))))
rs <- raster::aggregate(x, fact = dm/outdim)
} else {
rs <- rgdal::readGDAL(x@file@name, offset = offs, region.dim = dm, output.dim = outdim, silent = TRUE)
rs <- raster::raster(rs)
}
rs
}
collect(rastersource, nrows = 768, ncols = 768)
## rastersource is 8192x8192, collecting output of 768x768
## use 'nrows = , ncols = ' to customize output size
## class : RasterLayer
## dimensions : 768, 768, 589824 (nrow, ncol, ncell)
## resolution : 10922.67, 10922.67 (x, y)
## extent : -4194304, 4194304, -4194304, 4194304 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : band1
## values : 0, 243 (min, max)
## in-memory catch
#r <- raster(volcano)
#collect(r, nrows = 5, ncols = 16)
#collect(r)