This example requires the devtools, raster, and lubridate packages, which can be installed from CRAN. And finally the CubeR package can be installed from Eurac’s Gitlab, see also https://gitlab.inf.unibz.it/earth_observation_public/CubeR We use the development branch.

# devtools::install_git("https://gitlab.inf.unibz.it/earth_observation_public/CubeR", ref = "dev")

Load the required packages:

library(CubeR)

## Warning: replacing previous import 'magrittr::extract' by 'raster::extract'
## when loading 'CubeR'

## Warning: replacing previous import 'lubridate::origin' by 'raster::origin'
## when loading 'CubeR'

## Warning: replacing previous import 'urltools::url_parse' by
## 'xml2::url_parse' when loading 'CubeR'

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following object is masked from 'package:base':
## 
##     date

library(raster)

## Loading required package: sp

library(magrittr)

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:raster':
## 
##     extract

List all available datasets (the first ones):

getCapability() %>% head

## [1] "ERA5_Air_Temperature_2m_Europe_Hourly"              
## [2] "ERA5_Total_Sky_Direct_Solar_Radiation_Europe_Hourly"
## [3] "EURAC_SNOW_CLOUDREMOVAL_MODIS_ALPS_LAEA"            
## [4] "EURAC_SNOW_MODIS_ALPS_LAEA"                         
## [5] "MOSAIC_ALPS_NIR_R_G"                                
## [6] "MOSAIC_TEST"

Set to MODIS snow cover (choose original or cloudremoval)

coverage <- "EURAC_SNOW_MODIS_ALPS_LAEA"

Get the whole extent

bbox <- coverage_get_bounding_box(coverage = coverage)
print(bbox)

## [1] "3847097.79599"                 "4982445.60783"                
## [3] "2208288.965829979065247953306" "2872447.74497"

Show the temporal availability

coverage_get_temporal_extent(coverage = coverage)

## [1] "2002-07-03T12:00:00.000Z" "2019-06-09T12:00:00.000Z"

List bands

coverage_get_bands(coverage = coverage)

## [1] "SnowMap"     "SnowQuality"

Now, let’s download the SnowMap for the whole extent for an arbitrary date (e.g. 12 January 2018). The function will return a list of rasters (even if only one band and one date is selected). This might take a while, depending on the internet connection.

raster_list <- image_from_coverage(coverage = coverage,
                                   slice_E = bbox[1:2],
                                   slice_N = bbox[3:4],
                                   date = ymd("2018-01-12"),
                                   bands = "SnowMap")
raster_list

## [[1]]
## class      : RasterLayer 
## dimensions : 2867, 4901, 14051167  (nrow, ncol, ncell)
## resolution : 231.6564, 231.6564  (x, y)
## extent     : 3847098, 4982446, 2208289, 2872448  (xmin, xmax, ymin, ymax)
## crs        : +init=epsg:3035 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 5  (min, max)

And a simple plot

rr <- raster_list[[1]]
plot(rr)

Now let’s download a cloudremoval image. Since it has the same specificiation as the original maps, we only need to change the coverage. Specification of the band is optional, since the coverage has only one band.

coverage <- "EURAC_SNOW_CLOUDREMOVAL_MODIS_ALPS_LAEA"
bbox <- coverage_get_bounding_box(coverage = coverage)
rr_cloudremoval <- image_from_coverage(coverage = coverage,
                                       slice_E = bbox[1:2],
                                       slice_N = bbox[3:4],
                                       date = ymd("2018-01-12"))[[1]]

plot(rr_cloudremoval)

The function can also be used to access only parts of the data, by specifying different slices. E.g. let’s take only the central part between 4.2e6 and 4.6e6 for X, and 2.5e6 and 2.7e6 for Y.

coverage <- "EURAC_SNOW_CLOUDREMOVAL_MODIS_ALPS_LAEA"
bbox <- coverage_get_bounding_box(coverage = coverage)
rr_cloudremoval_sub <- image_from_coverage(coverage = coverage,
                                       slice_E = c(4.2, 4.6) * 10^6,
                                       slice_N = c(2.5, 2.7) * 10^6,
                                       date = ymd("2018-01-12"))[[1]]

plot(rr_cloudremoval_sub)