ECMWF land-cover ncdf examples

ecmwf-gdal.R

library(gdalcubes)
library(dplyr)
library(sf)
library(spData)
aoi <- spData::world |> filter(name_long=="United States") 

gdalcubes_options(parallel = 100) # can be way more than "real" cpus.

# netcdf is an inconsistent format where auto-detection of layers and CRS doesn't always work. We specify this crap in the URL:
urls <- 'vrt://NETCDF:"/vsicurl/https://minio.carlboettiger.info/public-ecmwf/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7cds.nc":lccs_class?a_srs=OGC:CRS84'

img <- gdalcubes::create_image_collection(urls, date_time = "2015-01-01")


# We can create "views" with arbitrary resolution in space and time and in arbitrary projections.
# We can do most band math operations in streaming fashion.  Can also do zonal stats. 

#box <- st_bbox(world)
box <- st_bbox(aoi)
view <- cube_view(srs = "EPSG:4326",
                  extent = list(t0 = "2015-01-01", t1 = "2015-12-31",
                                left = box[1], right = box[3],
                                top = box[4], bottom = box[2]),
                  #dx = .00277, dy=.00277, # original resolution
                  dx = .0277, dy=.0277,    # lower resolution for quick pass
                  dt = "P1Y")

raster_cube(img, view) |> plot(col = viridisLite::viridis)

raster_cube(img, view) |> write_ncdf("cube.ncdf")


# We use gdal_translate (gdal_warp if we need to reproject) to sample the 
sf::gdal_utils("translate", 
               'NETCDF:"/vsicurl/https://minio.carlboettiger.info/public-ecmwf/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7cds.nc":lccs_class',
               "lc-2015.tif",
               options = c("-a_srs", "EPSG:4326", "-of", "COG"))
urls <- "https://minio.carlboettiger.info/public-ecmwf/lc-2015.tif"



r <- terra::rast("/vsicurl/https://minio.carlboettiger.info/public-ecmwf/lc-2015.tif")
r # metadata is fine

## Do not try to plot without downscale!
# r |> plot()

# stars will auto-downscale:
library(stars)
r <- stars::read_stars("/vsicurl/https://minio.carlboettiger.info/public-ecmwf/lc-2015.tif")
r
r |> plot()