accessing spatial imagery via rgdal's web map tile service driver

Thu, Feb 8, 2018 3:56 PM

I've put a fuller example together, thanks for the prompt, it's nice to
know how to do this:

http://rpubs.com/cyclemumner/358029

Note that raster is using underlying rgdal offset/output.dim functionality
to drive GDAL itself here, something that's not exposed well - raster makes
it seamless, but the underlying tools are a bit obscure and could be used
more widely.

The example uses row/col indexing via the extent(RasterLayer, numeric, ...)
method, you can otherwise use any extent in the native (north Polar
Stereographic) coordinate system. I don't know if raster interprets
aggregate(, fact = ) calls into translations for rgdal's output.dim, but it
would be powerful to be easy to do that.

Cheers, Mike.

On Fri, 9 Feb 2018 at 10:32 Michael Sumner <mdsumner at gmail.com> wrote:

Oh, I had a typo in the xml_specification due to blithe copy/paste from
email. It works fine as intended.

I'll update and report back, this is a good example for fleshing out some
issues.

Cheers!


On Fri, 9 Feb 2018 at 08:46 Fischbach, Anthony <afischbach at usgs.gov>
wrote:

Mike, thank you for digging in on this.  I have submitted my
sessionInfo() and comments to your git script regarding an error on the
raster function call. - Tony

Anthony Fischbach, Wildlife Biologist
USGS Alaska Science Center Walrus Research Program
4210 University Drive
<https://maps.google.com/?q=4210+University+DriveAnchorage,+AK+99508&entry=gmail&source=g>
Anchorage, AK 99508
<https://maps.google.com/?q=4210+University+DriveAnchorage,+AK+99508&entry=gmail&source=g>

voice: (907) 786-7145

http://alaska.usgs.gov/science/biology/walrus/

On Thu, Feb 8, 2018 at 12:26 PM, Michael Sumner <mdsumner at gmail.com>
wrote:

Hi Anthony, what OS are you on? Can you share sessionInfo() or
devtools::session_info()?

But, it could be simpler, but  it doesn't work on my system, I'll find
out more:

https://gist.github.com/mdsumner/ca332f9c5112e00edcc40fdcf2a4b968

Cheers, Mike.

On Thu, 8 Feb 2018 at 06:28 Fischbach, Anthony <afischbach at usgs.gov>
wrote:

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]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

accessing spatial imagery via rgdal's web map tile service driver

Thread (5 messages)