ETCCDI индекс продолжительности теплого заклинания (wsdi) в RStudio с использованием оператора климатических данных (cdo)
Я пытаюсь вычислить ETCCDI индекс продолжительности теплого заклинания (wsdi) в RStudio, используя оператор климатических данных (cdo). Я использовал следующий скрипт и проверил результаты с помощью значений, доступных на веб-странице KNMI. Результаты не совпадают с доступными в KNMI. Все файлы от clippedtasmax_1981.nc до clippedtasmax_2100.nc (с единицами измерения в dC) находятся в папке input_folder. Не могли бы вы сказать мне, что не так с этим сценарием.
library(raster)
library(ncdf4)
startyear <- 1981
endyear <- 2100
for (i in startyear:endyear)
{
file_tmax <- paste(input_folder,"\\clippedtasmax_ ",i,".nc",sep="")
command <- paste("cdo timmin ",file_tmax," ",workdir,"tmax_min.nc",sep="")
system(command)
command <- paste("cdo timmax ",file_tmax," ",workdir,"tmax_max.nc",sep="")
system(command)
command <- paste("cdo timpctl,90 ",file_tmax," ",workdir,"tmax_min.nc ",workdir,"tmax_max.nc ",workdir,"temp_tmax_90.nc",sep="")
system(command)
command <- paste("cdo -add ",workdir,"temp_tmax_90.nc -sub ",file_tmax," ",file_tmax," ",workdir,"temp_tmax_90_timesteps.nc",sep = "")
system(command)
command <- paste("cdo eca_hwfi ",file_tmax," ",workdir,"temp_tmax_90_timesteps.nc ",outdir,"\\Warmspell_",i,".nc",sep = "")
system(command)
}
Затем я пересэмплировал выходные файлы "Warmspell_i.nc".
rm(list=ls())
library(raster)
#create new reference raster
referenceraster<-raster(xmn=80, xmx=90, ymn=25, ymx=30,
crs=CRS('+proj=longlat +datum=WGS84'))
res(referenceraster)<-2.5
referenceraster
#files to be resampled
wsdi<-"X:\\outdir\\warmspell_1981.nc"
wsdi<-raster(wsdi)
wsdi
resample_raster<-resample(wsdi, referenceraster, method="bilinear", nas.rm=TRUE)
values(resample_raster)
extent<-extent(80, 90, 25, 30)
e <- extract(resample_raster, extent, weights=TRUE, fun=mean)
e
Заранее спасибо.