Calculating Standardized Precipitation Index (SPI) in R

[Juweriya]

I am new to R and am trying to calculate the Standardized Precipitation Index (SPI) using the following code.

library(SPIGA)
library(readr)
SS <- read_csv("C:/Users/Aaisha/Desktop/Rain.csv", 
               col_types = cols(MONTH = col_integer(), 
                                year = col_integer())) # IMPORT DATA

newdata <- data.frame()   # MAKE A BLANK DATA FRAME
newdata <- data.frame(SS)   #PUT THE DATA IN THE NEW DATA FRAME
anyNA(newdata) #CHECK FOR MISSING DATA. TRUE = yes, data is missing. FALSE = no, data is not missing.
pob <-50 
iMax<-10
SPIML(newdata, scale =3)

when I run the code i get the error message

********** Calculation Station RF1 
Error in if ((Hx[i] > 0) & (Hx[i] <= 0.5)) { : 
  missing value where TRUE/FALSE needed
In addition: Warning message:
In pgamma(Pt, shape = alpha, scale = beta) : NaNs produced

I have ensured that there are no missing values.
Would appreciate any help to understand the error code.

[nirgrahamuk]

There is a relationship between you choice of scale and the minimum number of records/rows in your data.frame.
Using the SPIML documented example, which uses the data(Pm_Pisco); I saw if I reduced the number of entries below 26; while scaling at 3; it would error.
Do you have less than 26 entries?

scaling at 1 worked with 24 entries

[Juweriya]

Thank you for your reply.
Yes I have data from 1980 to 2019, for all 12 months. So I have 480 rows.
Also I want to calculate the SPI for multiple geographical locations, so I have arranged the data into seperate columns.

[nirgrahamuk]

Then I don't know what your problem is;
Could you make a reprex ?

[Juweriya]

I used reprex to produce the following

library(reprex)

library(SPIGA)
#> Loading required package: GA
#> Loading required package: foreach
#> Loading required package: iterators
#> Package 'GA' version 3.2.3
#> Type 'citation("GA")' for citing this R package in publications.
#> 
#> Attaching package: 'GA'
#> The following object is masked from 'package:utils':
#> 
#>     de
library(readr)
 SS <- read_csv("C:/Users/Aaisha/Desktop/Rain.csv", 
               col_types = cols(MONTH = col_integer(), 
                                year = col_integer())) # IMPORT DATA
newdata <- data.frame()   # MAKE A BLANK DATA FRAME
newdata <- data.frame(SS)   #PUT THE DATA IN THE NEW DATA FRAME
anyNA(newdata) #CHECK FOR MISSING DATA. TRUE = yes, data is missing. FALSE = no, data is not missing.
#> [1] FALSE
pob <-50 
iMax<-10
SPIML(newdata, scale =3)
#> ********** Calculation Station RF1
#> Warning in pgamma(Pt, shape = alpha, scale = beta): NaNs produced
#> Error in if ((Hx[i] > 0) & (Hx[i] <= 0.5)) {: missing value where TRUE/FALSE needed

reprex()
#> ℹ Non-interactive session, setting `html_preview = FALSE`.
#> CLIPR_ALLOW has not been set, so clipr will not run interactively
#> Error in switch(where, expr = stringify_expression(x_expr), clipboard = ingest_clipboard(), : EXPR must be a length 1 vector

[nirgrahamuk]

That didnt work...
we don't have your csv.
start your reprex from when you construct newdata; and share with us a representation of newdata

FAQ: How to do a minimal reproducible example ( reprex ) for beginners

A minimal reproducible example consists of the following items: A minimal dataset, necessary to reproduce the issue The minimal runnable code necessary to reproduce the issue, which can be run on the given dataset, and including the necessary information on the used packages. Let's quickly go over each one of these with examples: Minimal Dataset (Sample Data) You need to provide a data frame that is small enough to be (reasonably) pasted on a post, but big enough to reproduce your issue. Let's say, as an example, that you are working with the iris data frame head(iris) #> Sepal.Length Sepal.Width Petal.Length Petal.Width Species #> 1 5.1 3.5 1.4 0.…

[Juweriya]

Thank you for your reply.
This is my data and the code with the error.

dput(head(SS, 200)[c("year", "RF1", "RF2")])
structure(list(year = c(1980L, 1980L, 1980L, 1980L, 1980L, 1980L,
1980L, 1980L, 1980L, 1980L, 1980L, 1980L, 1981L, 1981L, 1981L,
1981L, 1981L, 1981L, 1981L, 1981L, 1981L, 1981L, 1981L, 1981L,
1982L, 1982L, 1982L, 1982L, 1982L, 1982L, 1982L, 1982L, 1982L,
1982L, 1982L, 1982L, 1983L, 1983L, 1983L, 1983L, 1983L, 1983L,
1983L, 1983L, 1983L, 1983L, 1983L, 1983L, 1984L, 1984L, 1984L,
1984L, 1984L, 1984L, 1984L, 1984L, 1984L, 1984L, 1984L, 1984L,
1985L, 1985L, 1985L, 1985L, 1985L, 1985L, 1985L, 1985L, 1985L,
1985L, 1985L, 1985L, 1986L, 1986L, 1986L, 1986L, 1986L, 1986L,
1986L, 1986L, 1986L, 1986L, 1986L, 1986L, 1987L, 1987L, 1987L,
1987L, 1987L, 1987L, 1987L, 1987L, 1987L, 1987L, 1987L, 1987L,
1988L, 1988L, 1988L, 1988L, 1988L, 1988L, 1988L, 1988L, 1988L,
1988L, 1988L, 1988L, 1989L, 1989L, 1989L, 1989L, 1989L, 1989L,
1989L, 1989L, 1989L, 1989L, 1989L, 1989L, 1990L, 1990L, 1990L,
1990L, 1990L, 1990L, 1990L, 1990L, 1990L, 1990L, 1990L, 1990L,
1991L, 1991L, 1991L, 1991L, 1991L, 1991L, 1991L, 1991L, 1991L,
1991L, 1991L, 1991L, 1992L, 1992L, 1992L, 1992L, 1992L, 1992L,
1992L, 1992L, 1992L, 1992L, 1992L, 1992L, 1993L, 1993L, 1993L,
1993L, 1993L, 1993L, 1993L, 1993L, 1993L, 1993L, 1993L, 1993L,
1994L, 1994L, 1994L, 1994L, 1994L, 1994L, 1994L, 1994L, 1994L,
1994L, 1994L, 1994L, 1995L, 1995L, 1995L, 1995L, 1995L, 1995L,
1995L, 1995L, 1995L, 1995L, 1995L, 1995L, 1996L, 1996L, 1996L,
1996L, 1996L, 1996L, 1996L, 1996L), RF1 = c(6.867742, 24.52051,
48.88, 112.8225, 208.11, 419.7175, 512.065, 443.455, 306.31,
191.9625, 61.86842, 14.29394, 6.867742, 24.52051, 48.88, 112.8225,
208.11, 419.7175, 512.065, 443.455, 306.31, 191.9625, 61.86842,
14.29394, 6.867742, 24.52051, 48.88, 112.8225, 208.11, 419.7175,
512.065, 443.455, 306.31, 191.9625, 61.86842, 14.29394, 6.867742,
24.52051, 48.88, 112.8225, 208.11, 419.7175, 512.065, 443.455,
306.31, 191.9625, 61.86842, 14.29394, 6.867742, 24.52051, 48.88,
112.8225, 208.11, 419.7175, 512.065, 443.455, 306.31, 191.9625,
61.86842, 14.29394, 6.867742, 24.52051, 48.88, 112.8225, 208.11,
419.7175, 512.065, 443.455, 306.31, 191.9625, 61.86842, 14.29394,
6.867742, 24.52051, 48.88, 112.8225, 208.11, 419.7175, 512.065,
443.455, 306.31, 191.9625, 61.86842, 14.29394, 6.867742, 24.52051,
48.88, 112.8225, 208.11, 419.7175, 512.065, 443.455, 306.31,
191.9625, 61.86842, 14.29394, 6.867742, 24.52051, 48.88, 112.8225,
208.11, 419.7175, 512.065, 443.455, 306.31, 191.9625, 61.86842,
14.29394, 6.867742, 24.52051, 48.88, 112.8225, 208.11, 419.7175,
512.065, 443.455, 306.31, 191.9625, 61.86842, 14.29394, 6.867742,
24.52051, 48.88, 112.8225, 208.11, 419.7175, 512.065, 443.455,
306.31, 191.9625, 61.86842, 14.29394, 6.867742, 24.52051, 48.88,
112.8225, 208.11, 419.7175, 512.065, 443.455, 306.31, 191.9625,
61.86842, 14.29394, 6.867742, 24.52051, 48.88, 112.8225, 208.11,
419.7175, 512.065, 443.455, 306.31, 191.9625, 61.86842, 14.29394,
6.867742, 24.52051, 48.88, 112.8225, 208.11, 419.7175, 512.065,
443.455, 306.31, 191.9625, 61.86842, 14.29394, 6.867742, 24.52051,
48.88, 112.8225, 208.11, 419.7175, 512.065, 443.455, 306.31,
191.9625, 61.86842, 14.29394, 6.867742, 24.52051, 48.88, 112.8225,
208.11, 419.7175, 512.065, 443.455, 306.31, 191.9625, 61.86842,
14.29394, 6.867742, 24.52051, 48.88, 112.8225, 208.11, 419.7175,
512.065, 443.455), RF2 = c(6.661835, 24.51724, 48.86134, 112.7645,
208.015, 419.5669, 511.7952, 443.2962, 306.2538, 191.8952, 61.8303,
14.29017, 6.661835, 24.51724, 48.86134, 112.7645, 208.015, 419.5669,
511.7952, 443.2962, 306.2538, 191.8952, 61.8303, 14.29017, 6.661835,
24.51724, 48.86134, 112.7645, 208.015, 419.5669, 511.7952, 443.2962,
306.2538, 191.8952, 61.8303, 14.29017, 6.661835, 24.51724, 48.86134,
112.7645, 208.015, 419.5669, 511.7952, 443.2962, 306.2538, 191.8952,
61.8303, 14.29017, 6.661835, 24.51724, 48.86134, 112.7645, 208.015,
419.5669, 511.7952, 443.2962, 306.2538, 191.8952, 61.8303, 14.29017,
6.661835, 24.51724, 48.86134, 112.7645, 208.015, 419.5669, 511.7952,
443.2962, 306.2538, 191.8952, 61.8303, 14.29017, 6.661835, 24.51724,
48.86134, 112.7645, 208.015, 419.5669, 511.7952, 443.2962, 306.2538,
191.8952, 61.8303, 14.29017, 6.661835, 24.51724, 48.86134, 112.7645,
208.015, 419.5669, 511.7952, 443.2962, 306.2538, 191.8952, 61.8303,
14.29017, 6.661835, 24.51724, 48.86134, 112.7645, 208.015, 419.5669,
511.7952, 443.2962, 306.2538, 191.8952, 61.8303, 14.29017, 6.661835,
24.51724, 48.86134, 112.7645, 208.015, 419.5669, 511.7952, 443.2962,
306.2538, 191.8952, 61.8303, 14.29017, 6.661835, 24.51724, 48.86134,
112.7645, 208.015, 419.5669, 511.7952, 443.2962, 306.2538, 191.8952,
61.8303, 14.29017, 6.661835, 24.51724, 48.86134, 112.7645, 208.015,
419.5669, 511.7952, 443.2962, 306.2538, 191.8952, 61.8303, 14.29017,
6.661835, 24.51724, 48.86134, 112.7645, 208.015, 419.5669, 511.7952,
443.2962, 306.2538, 191.8952, 61.8303, 14.29017, 6.661835, 24.51724,
48.86134, 112.7645, 208.015, 419.5669, 511.7952, 443.2962, 306.2538,
191.8952, 61.8303, 14.29017, 6.661835, 24.51724, 48.86134, 112.7645,
208.015, 419.5669, 511.7952, 443.2962, 306.2538, 191.8952, 61.8303,
14.29017, 6.661835, 24.51724, 48.86134, 112.7645, 208.015, 419.5669,
511.7952, 443.2962, 306.2538, 191.8952, 61.8303, 14.29017, 6.661835,
24.51724, 48.86134, 112.7645, 208.015, 419.5669, 511.7952, 443.2962
)), row.names = c(NA, -200L), class = c("tbl_df", "tbl", "data.frame"
))

library(reprex)

library(SPIGA)
#> Loading required package: GA
#> Loading required package: foreach
#> Loading required package: iterators
#> Package 'GA' version 3.2.3
#> Type 'citation("GA")' for citing this R package in publications.
#> 
#> Attaching package: 'GA'
#> The following object is masked from 'package:utils':
#> 
#>     de
library(readr)
 SS <- read_csv("C:/Users/Aaisha/Desktop/Rain.csv", 
               col_types = cols(MONTH = col_integer(), 
                                year = col_integer())) # IMPORT DATA
newdata <- data.frame()   # MAKE A BLANK DATA FRAME
newdata <- data.frame(SS)   #PUT THE DATA IN THE NEW DATA FRAME
anyNA(newdata) #CHECK FOR MISSING DATA. TRUE = yes, data is missing. FALSE = no, data is not missing.
#> [1] FALSE
pob <-50 
iMax<-10
SPIML(newdata, scale =3)
#> ********** Calculation Station RF1
#> Warning in pgamma(Pt, shape = alpha, scale = beta): NaNs produced
#> Error in if ((Hx[i] > 0) & (Hx[i] <= 0.5)) {: missing value where TRUE/FALSE needed

reprex()
#> ℹ Non-interactive session, setting `html_preview = FALSE`.
#> CLIPR_ALLOW has not been set, so clipr will not run interactively
#> Error in switch(where, expr = stringify_expression(x_expr), clipboard = ingest_clipboard(), : EXPR must be a length 1 vector