R 标准化降水指数 .nc 文件

Question

我正在尝试从 CHIRPS 月平均降水量数据中计算 SPI，因为它太大了，我将其缩减到我感兴趣的区域，这里是：https://www.dropbox.com/s/jpwcg8j5bdc5gq6/chirps_mensual_v1.nc?dl=0我这样做是为了打开它：

require(utils)
require(colorRamps)
require(RNetCDF)
require(rasterVis)
require(rgdal)
library(ncdf4)
library(raster)

datos2 <- nc_open("Datos/chirps_mensual_v1.nc")
ppt_array <- ncvar_get(datos2, "precip")
#I'm only taking complete years so I took out two months from 2018
ppt_mes <- ppt_array[ , ,1:444]

我知道有一个 SPI 库，但我不知道应该如何格式化数据才能使用它。所以我试图通过拟合伽马分布来在没有函数的情况下做到这一点，但我不知道如何为这个数据库做到这一点。

有谁知道如何使用函数或通过拟合分布来计算 SPI？

Answer 1

我不认为 SPI 包正在做你(或任何人(认为它正在做的事情。 如果使用debug(spi)并逐步执行代码，您将看到它在一个步骤中将经验累积分布函数(带 ecdf() (拟合到前两行和最后一行数据。 为什么是前两行和最后一行？ 我不知道，但是编写此包的人也使用for循环对矩阵执行t()。 更不用说我认为它应该使用伽马分布或皮尔逊III分布而不是ecdf()(根据Guttman，N.B.(1999(接受标准化降水指数：一种计算算法。JAWRA美国水资源协会杂志，35,311-322。

Answer 2

最后，我使用 SPI 库制作了它，结果将是每个网格点中每个月的值，如果您想计算特定区域的值，我也做了，但如果您也想要，我可以分享它：

另外，这个我使用 CRU 数据制作了它，但您可以调整它：

#spei cru 1x1
rm(list=ls(all=TRUE)); dev.off()
require(utils)
require(RNetCDF)
require(rasterVis)
require(rgdal)
library(ncdf4)
require(SPEI)
########################################################################################################

prec <- open.nc("pre_mensual.nc")
lon <- length(var.get.nc(prec, "lon"))
lat <- length(var.get.nc(prec, "lat"))
lon1 <- var.get.nc(prec, "lon")
lat1 <- var.get.nc(prec, "lat")
ppt  <- var.get.nc(prec, "pre") 
ppt  <- ppt[ , ,109:564] #31 18 456 (1980-2017)
anio = 456/12
###########################################################################################################
#Reshape data 
precip <- sapply(1:dim(ppt)[3], function(x)t(ppt[,,x]))
############################################################################################################
#This is for SPI-6, you can use either of them
spi_6 <- array(list(),(lon*lat))
for (i in 1:(lon*lat)) {
  spi_6[[i]] <- spi(precip[i,], scale=6, na.rm=TRUE)
}
#############################################################################################################
#Go back to an array form
sapply(spi_6, '[[',2 )->matriz_ppt 
ppt_6 <- array(aperm(matriz_ppt, c(2,1),c(37,63,456)));spi_c <- array(t(ppt_6), dim=c(37,63,456))
#############################################################################################################
    #Save to netcdf
for(i in 1:456) { 
  nam <- paste("SPI", i, sep = "")
  assign(nam,raster((spi_c[ , ,i]), xmn=min(lon1), xmx=max(lon1), ymn=min(lat1), ymx=max(lat1), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0")) )
}
gpcc_spi <- stack(mget(paste0("SPI", 1:456)))
outfile <- "spi6_cru_1980_2017.nc"
crs(gpcc_spi) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 
writeRaster(gpcc_spi, outfile, overwrite=TRUE, format="CDF", varname="SPEI", varunit="units",longname="SPEI CRU", xname="lon", yname="lat")

这不是最时尚的计算方法，但它确实有效。 :)

编辑：如果你想计算一个面积上的SPI/SPEI，这就是我所做的：

library(SPEI)
library(ncdf4)
library(raster)
#
pre_nc <- nc_open("pre_1971_2017_Vts4.nc")
pre <- ncvar_get(pre_nc, "pre")
pre <- pre[, , 109:564] #This is for the time I'm interested in
lats <- ncvar_get(pre_nc, "lat")
lons <- ncvar_get(pre_nc, "lon")
times <- 0:467
# Read mask
#This is a mask you need to create that adjusts to your region of interest
#It consist of a matrix of 0's and 1's, the 1's are placed in the area
#you are interested in
mask1 <- nc_open("cuenca_IV_CDO_05_final.nc")
m1 <- ncvar_get(mask1, "Band1")
m1[m1 == 0] <- NA
#
# Apply mask to data
#
pre1 <- array(NA, dim=dim(pre))
#
for(lon in 1:length(lons)){
  for(lat in 1:length(lats)){
    pre1[lon,lat,] <- pre[lon,lat,]*m1[lon,lat]
  } 
}
#
# Mean over the area of interest
#
mean_pre1 <- apply(pre1,c(3),mean, na.rm=TRUE)
# Calculate SPI/SPEI
spi1 <- matrix(data= NA, nrow = 456, ncol = 48)
for (i in 1:48) {
  spi1[,i] <- spi(data=ts(mean_pre1,freq=12),scale= i)$fitted
}
#This calculates SPI/SPEI-1 to SPI/SPEI-48, you can change it
# Save
#
write.table(spi1,'spi_1980_2017.csv',sep=';',row.names=FALSE)