[R-sig-Geo] accessing spatial imagery via rgdal's web map tile service driver

[Fischbach, Anthony]

The GDAL tile web map service appears to require an xml specification for
the TMS and a "window" specification indicating the bounds of the region to
be accessed, specified in the projected coordinate system using the
-projwin flag. I am uncertain how to pass this -projwin and its bounds to
the rgdal::readGDAL function.  The code snippet below stages the problem.
Please advise on how to acquire just the "window" of interest.

#Tile Web Mapping Service to earthdata imagery via rgdal
library(rgdal)  # r implementation of GDAL
require(raster) # r package for raster image manipulation, required for
raster capture of the readGDAL function result.
require(sp)     # r package for spatial class objects, including CRS
#
prj.geo<-CRS('+proj=longlat +datum=WGS84') # Geographic projection
prj.Polar_NSIDC <- CRS("+init=epsg:3413")
# (+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)
prj.StudyArea <- CRS("+proj=laea +lat_0=70 +lon_0=-170 +x_0=0 +y_0=0
+datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") # study area
projection
##
xLimits <- sort(c(170, -135)) # define the study area polygon using
geographic coordinates (eastings WGS84)
yLimits <- sort(c(52, 75)) # (northings WGS84)
pts<-expand.grid(xLimits, yLimits)
names(pts)<-c('x','y')
coordinates(pts) <- ~x+y
proj4string(pts) <- prj.geo
pts.Polar_NSIDC<-spTransform(pts, prj.Polar_NSIDC)
bb <- round(bbox(pts.Polar_NSIDC)) #
#
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>
')
#
cat(xml_specification)
GDALinfo(fname = xml_specification)
 (projWinParameters = paste('-projwin', bb[1, 1], bb[2,2], bb[1,2],
bb[2,1], sep = " "))
aqua <- readGDAL(fname = xml_specification) # , projWinParameters) ?? how
to implement the projwin constraints??
(aqua_stack <- stack(aqua))
#
plotRGB(aqua_stack) # plots the acquired image in RGB !! problem is that
the entire WMS domain is acquired.

	[[alternative HTML version deleted]]