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 Anchorage, AK 99508 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> > > [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo