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 <afischb...@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 <mdsum...@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 <afischb...@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@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 >> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g> >> Kingston Tasmania 7050 Australia >> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g> >> >> > -- Dr. Michael Sumner Software and Database Engineer Australian Antarctic Division 203 Channel Highway Kingston Tasmania 7050 Australia [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo