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

[Michael Sumner]

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

	[[alternative HTML version deleted]]