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 <mdsum...@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 <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 > > -- 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