"Debugging" lazyraster

[baerbock]

Dear hypertidy,

how is the gdal API exactly called at the TMS example ( http://rpubs.com/cyclemumner/358029 )?

I would like to understand this a little bit deeper.

[mdsumner]

That actually uses rgdal, I can't remember really the link I was thinking to lazyraster ... I will have a look and flesh it out some more 😉

[mdsumner]

In terms of lazyraster you should ignore that Rpubs post, the lazyraster readme has an example but here I use the same source as on Rpubs:

  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>
         ')

library(lazyraster)
lr <- lazyraster(xml_specification)
## we can see the overall metadata
lr
#> class         : LazyRaster
#> dimensions    : 8192, 8192 (nrow, ncol)
#> resolution    : 1024.0000, 1024.0000 (x, y)
#> extent        : -4194304.0000,  4194304.0000, -4194304.0000,  4194304.0000 (xmin, xmax, ymin, ymax)
#> crs           : <placeholder>
#> values        :   0.0000, 230.0000 (min, max - range from entire extent)
#> window extent : <whole extent>
#> window index  : <->

## and crop to a smaller extent 
(xr <- lazycrop(lr, raster::extent(0,  4194304.0000, 0,  4194304.0000)))
#> class         : LazyRaster
#> dimensions    : 8192, 8192 (nrow, ncol)
#> resolution    : 1024.0000, 1024.0000 (x, y)
#> extent        : -4194304.0000,  4194304.0000, -4194304.0000,  4194304.0000 (xmin, xmax, ymin, ymax)
#> crs           : <placeholder>
#> values        :   0.0000, 230.0000 (min, max - range from entire extent)
#> window extent :       0.0000, 4194304.0000,       0.0000, 4194304.0000
#> window index  : 4096, 8191, 0, 4096

## but no data is pulled until we ask for it with plot() or as_raster()
as_raster(lr, dim = c(300, 300))
#> class       : RasterLayer 
#> dimensions  : 300, 300, 90000  (nrow, ncol, ncell)
#> resolution  : 27962.03, 27962.03  (x, y)
#> extent      : -4194304, 4194304, -4194304, 4194304  (xmin, xmax, ymin, ymax)
#> coord. ref. : NA 
#> data source : in memory
#> names       : layer 
#> values      : 0, 225  (min, max)

^{Created on 2019-04-06 by the reprex package (v0.2.1)}

The actual GDAL stuff is driven by vapour, down in this internal function where it has figured out the window specification and vapour gets the raw data from GDAL, GDAL handles all the rest with its RasterIO facility down in the C++.

https://github.com/hypertidy/lazyraster/blob/d3818f3d825ab02f01f063f003047d75668ddc49/R/lazyraster.R#L148-L149

Does that clarify? Lazyraster isn't very function still, but I use it a lot to get data quickly. stars can do this now too, with read_stars(, proxy = TRUE) by using internal GDAL code in sf:

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>
       ')

stars::read_stars(xml_specification, proxy = TRUE)
#> stars_proxy object with 1 attribute in file:
#> $`GDAL_WMS>\n       `
#> [1] "<GDAL_WMS>\n       <Service name=\"TMS\">\n       <ServerUrl>\n       https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/MODIS_Aqua_CorrectedReflectance_Bands721/default/2019-04-05/250m/${z}/${y}/${x}.jpg</ServerUrl>\n       </Service>\n       <DataWindow>\n       <UpperLeftX>-4194304</UpperLeftX>\n       <UpperLeftY>4194304</UpperLeftY>\n       <LowerRightX>4194304</LowerRightX>\n       <LowerRightY>-4194304</LowerRightY>\n       <TileLevel>3</TileLevel>\n       <TileCountX>2</TileCountX>\n       <TileCountY>2</TileCountY>\n       <YOrigin>top</YOrigin>\n       </DataWindow>\n       <Projection>EPSG:3413</Projection>\n       <BlockSizeX>512</BlockSizeX>\n       <BlockSizeY>512</BlockSizeY>\n       <BandsCount>3</BandsCount>\n       </GDAL_WMS>\n       "
#> 
#> dimension(s):
#>      from   to   offset delta                       refsys point values
#> x       1 8192 -4194304  1024 +proj=stere +lat_0=90 +la...    NA   NULL
#> y       1 8192  4194304 -1024 +proj=stere +lat_0=90 +la...    NA   NULL
#> band    1    3       NA    NA                           NA    NA   NULL
#>         
#> x    [x]
#> y    [y]
#> band

^{Created on 2019-04-06 by the reprex package (v0.2.1)}