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')) } }哥白尼陆地监测服务（Copernicus Land Monitoring Services）提供2014年至今的欧洲区域每日地表土壤湿度栅格数据（分辨率1 km，版本1）。该数据集每日仅覆盖欧洲陆地区域的5%至10%，且存在场景线条状明显伪影。 本数据集为基于两类聚合方式生成的土壤湿度（取值范围0–100%）长期聚合产品，原始数据为每日土壤湿度影像： 1. 长期季度聚合（其中qr.1对应冬季，qr.2对应春季，qr.3对应夏季，qr.4对应秋季）； 2. 年度分位数聚合，包含P.05、P.50与P.95三个分位数。 该土壤湿度栅格数据基于哨兵1号（Sentinel 1）卫星数据生成，详细技术说明参见以下文献： Bauer-Marschallinger, B. ; Freeman, V. ; Cao, S. ; Paulik, C. ; Schaufler, S. ; Stachl, T. ; Modanesi, S. ; Massari, C. ; Ciabatta, L. ; Brocca, L. ; Wagner, W. 基于哨兵1号的全球土壤湿度监测：充分利用优势与克服挑战. 《IEEE地球科学与遥感汇刊》, 2019, 1-20. DOI: 10.1109/TGRS.2018.2858004 用户可通过以下链接以.netCDF格式获取并下载原始数据：https://globalland.vito.be/download/manifest/ssm_1km_v1_daily_netcdf/ 文件名匹配格式"soil.moisture_s1.clms.qr.*.p0.*.gf_m_1km_20140101_20241231_eu_epsg4326_v20250211.tif"的文件为经间隙填充的季度土壤湿度估算产品。本间隙填充模型通过约25万个随机训练点构建，关联了CHELSA气候生物气候变量、欧洲空间局气候变化倡议（ESA CCI）积雪覆盖概率、ESA CCI森林与裸地百分比覆盖率以及全球水产品（Global Water Pack）长期地表水占比等辅助数据。该间隙填充模型的决定系数（R²）为0.96，均方根误差（RMSE）为土壤湿度幅值的6.5%。 聚合过程通过R语言的terra包结合matrixStats::rowQuantiles函数实现。泛欧洲区域的瓦片系统与陆地区域掩码同样可供获取。 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=" ", 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=" ", 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')) } }