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/.

The files with pattern "soil.moisture_s1.clms.qr.*.p0.*.gf_m_1km_20140101_20241231_eu_epsg4326_v20250211.tif" are the gap-filled soil moisture quarterly estimates. For gap filling I build a model using cca 250k random training points and relationship with CHELSA climate bioclimatic variables, ESA CCI snow cover probability, ESA CCI forest and bare areas percent cover and Global Water Pack long-term surface water fraction. The gap-filling model had an R-square of 0.96 and RMSE of 6.5% of soil moisture.

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))
    r.t = as.data.frame(r.t, xy=TRUE, na.rm=FALSE)
    sel.c = grep(glob2rx("ssm$"), colnames(r.t))
    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))
    ## 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])
  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)
  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'))
  }
}