Surface soil moisture for Europe 2014-2024 at 1 km annual and quarterly aggregates

Copernicus Land Monitoring Services provides Surface Soil Moisture 2014-present (raster 1 km), Europe, daily – version 1. Each day covers only 5 to 10% of European land mask and shows lines of scenes (obvious artifacts). This is the long-term aggregates of daily images of soil moisture (0–100%) based on two types of aggregation:

• Long-term quarterly (qr.1 - winter, qr.2 - spring, qr.3 - summer and qr.4 - autumn),
• Annual quantiles P.05, P.50 and P.95,

The soil moisture rasters are based on Sentinel 1 and described in detail in: 

• Bauer-Marschallinger, B. ; Freeman, V. ; Cao, S. ; Paulik, C. ; Schaufler, S. ; Stachl, T. ; Modanesi, S. ; Massari, C. ; Ciabatta, L. ; Brocca, L. ; Wagner, W. Toward Global Soil Moisture Monitoring With Sentinel-1: Harnessing Assets and Overcoming Obstacles. IEEE Transactions on Geoscience and Remote Sensing 2019, 1 - 20. DOI 10.1109/TGRS.2018.2858004

You can access and download the original data as .nc files from: https://globalland.vito.be/download/manifest/ssm_1km_v1_daily_netcdf/.

Aggregation has been generated using the terra package in R in combination with the matrixStats::rowQuantiles function. Tiling system and land mask for pan-EU is also available.

library(terra)
library(matrixStats)
g1 = terra::vect("/mnt/inca/EU_landmask/tilling_filter/eu_ard2_final_status.gpkg")
## 1254 tiles
tile = g1[534]
nc.lst = list.files('/mnt/landmark/SM1km/ssm_1km_v1_daily_netcdf/', pattern = glob2rx("*.nc$"), full.names=TRUE)
## 3726
## test it
#r = terra::rast(nc.lst[100:210])

agg_tile = function(r, tile, pv=c(0.05,0.5,0.95), out.year="2015.annual"){
  bb = paste(as.vector(ext(tile)), collapse = ".")
  out.tif = paste0("./eu_tmp/", out.year, "/sm1km_", pv, "_", out.year, "_", bb, ".tif")
  if(any(!file.exists(out.tif))){
    r.t = terra::crop(r, ext(tile))
    ## each tile is 100x100 pixels 365 days
    r.t = as.data.frame(r.t, xy=TRUE, na.rm=FALSE)
    sel.c = grep(glob2rx("ssm$"), colnames(r.t))
    ## remove everything outside the range
    t1s = cbind(data.frame(matrixStats::rowQuantiles(as.matrix(r.t[,sel.c]), probs = pv, na.rm=TRUE)),  data.frame(x=r.t$x,  y=r.t$y))
    #str(t1s)
    ## write to GeoTIFFs
    r.o = terra::rast(t1s[,c("x","y","X5.","X50.","X95.")], type="xyz", crs="+proj=longlat +datum=WGS84 +no_defs")
    for(k in 1:length(pv)){ 
      terra::writeRaster(r.o[[k]], filename=out.tif[k], gdal=c("COMPRESS=DEFLATE"), datatype='INT2U', NAflag=32768, overwrite=FALSE)
    }
    rm(r.t); gc()
    tmpFiles(remove=TRUE)
  }
}

## quarterly values:
lA = data.frame(filename=nc.lst)
library(lubridate)
lA$Date = ymd(sapply(lA$filename, function(i){substr(strsplit(basename(i), "_")[[1]][4], 1, 8)}))
#summary(is.na(lA$Date))
#hist(lA$Date, breaks=60)
lA$quarter = quarter(lA$Date, fiscal_start = 11)
summary(as.factor(lA$quarter))

for(qr in 1:4){
  #qr=1
  pth = paste0("A.q", qr)
  rs = terra::rast(lA$filename[lA$quarter==qr])
  #agg_tile(rs, tile, out.year=pth)
  x = parallel::mclapply(sample(1:length(g1)), function(i){try( agg_tile(rs, tile=g1[i], out.year=pth) )}, mc.cores=20)
  for(type in c(0.05,0.5,0.95)){
    x <- list.files(path=paste0("./eu_tmp/", pth), pattern=glob2rx(paste0("sm1km_", type, "_*.tif$")), full.names=TRUE)
    out.tmp <- paste0(pth, ".", type, ".sm1km_eu.txt")
    vrt.tmp <- paste0(pth, ".", type, ".sm1km_eu.vrt")
    cat(x, sep="\n", file=out.tmp)
    system(paste0('gdalbuildvrt -input_file_list ', out.tmp, ' ', vrt.tmp))
    system(paste0('gdal_translate ', vrt.tmp, ' ./cogs/soil.moisture_s1.clms.qr.', qr, '.p', type, '_m_1km_20140101_20241231_eu_epsg4326_v20250206.tif -ot "Byte" -r "near" --config GDAL_CACHEMAX 9216 -co BIGTIFF=YES -co NUM_THREADS=80 -co COMPRESS=DEFLATE -of COG -projwin -32 72 45 27'))
  }
}

## per year ----
for(year in 2015:2023){
  l.lst = nc.lst[grep(year, basename(nc.lst))]
  r = terra::rast(l.lst)
  ## test it:
  pth = paste0(year, ".annual")
  x = parallel::mclapply(sample(1:length(g1)), function(i){try( agg_tile(r, tile=g1[i], out.year=pth) )}, mc.cores=40)
  ## Mosaics:
  for(type in c(0.05,0.5,0.95)){
    x <- list.files(path=paste0("./eu_tmp/", pth), pattern=glob2rx(paste0("sm1km_", type, "_*.tif$")), full.names=TRUE)
    out.tmp <- paste0(pth, ".", type, ".sm1km_eu.txt")
    vrt.tmp <- paste0(pth, ".", type, ".sm1km_eu.vrt")
    cat(x, sep="\n", file=out.tmp)
    system(paste0('gdalbuildvrt -input_file_list ', out.tmp, ' ', vrt.tmp))
    system(paste0('gdal_translate ', vrt.tmp, ' ./cogs/soil.moisture_s1.clms.annual.', type, '_m_1km_', year, '0101_', year, '1231_eu_epsg4326_v20250206.tif -ot "Byte" -r "near" --config GDAL_CACHEMAX 9216 -co BIGTIFF=YES -co NUM_THREADS=80 -co COMPRESS=DEFLATE -of COG -projwin -32 72 45 27'))
  }
}