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

[Fischbach, Anthony]

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