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

Michael Sumner mdsumner at gmail.com
Fri Feb 9 00:56:08 CET 2018


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

	[[alternative HTML version deleted]]



More information about the R-sig-Geo mailing list