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

[Michael Sumner]

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 <afischbach at 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 <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>
>>
>>
> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]