0

I am a beginner in R, I need to use gstat or idw function to get a rainfall raster, but I have an error in the formula. How can I change it? This is my code:

setwd("C:/Users/A/Desktop/DatiPrecipitazioni/dati")
mydata<-read.csv(file="prova.csv",header=TRUE,sep=";",dec=".",na="")

dataset<-read.csv(file="Dataset.csv",header=TRUE,sep=";",dec=".",na="")

DTM=raster("DTM.tif")

vals<-mydata$HPrec

simplerasterPREC = function(dataset,DTM,vals){

 coordinates(dataset) = ~x+y

 proj4string(dataset)<- proj4string(DTM)

 idw_formula = gstat(id="vals",formula = as.formula(paste(vals,"~1")), data = dataset,nmax=10, maxdist=50000, set = list(idp=2)) 

 rast_int <- interpolate(DTM, idw_formula)

 rast_int <- mask(rast_int,DTM)

 return(rast_int)
}

If I see what happens line by line I have this error:

coordinates(dataset) = ~x+y

proj4string(dataset)<- proj4string(DTM)

idw_formula = gstat(id="vals",formula = as.formula(paste(vals,"~1")), data = dataset,nmax=10, maxdist=50000, set = list(idp=2))

Error in terms.formula(formula) : invalid term in model formula
  • Where comes mydata? – aldo_tapia Nov 14 '17 at 15:42
  • mydata is a csv file that contains the precipitation data and coordinates of the points where they are measured. It is in a folder on my desktop. – Fabiana Nov 14 '17 at 15:49
  • If mydata is the same one than dataset, change vals with HPrec in the function and try it. (gstat(id="HPrec",formula = HPrec ~ 1, data = dataset,nmax=10, maxdist=50000, set = list(idp=2))) – aldo_tapia Nov 14 '17 at 15:53
  • In dataset I have only coordinates, in mydata I have the coordinates and precipitation points. If I change vals with Hprec anyway it doesn't return rast_int. If I do (plot(rast_int)), I have this error: – Fabiana Nov 14 '17 at 16:06
  • I have an error: rast_int not found – Fabiana Nov 14 '17 at 16:07
1

I've done it like this: first make a grid to which you'll interpolate your met station values,

library(gstat)
library(ggplot2)
library(rgdal)

grd <- expand.grid(x=seq(from=bbox(met_stations)[1,1], to=bbox(met_stations)[1,2], by=0.05), 
                   y=seq(from=bbox(met_stations)[2,1], to=bbox(met_stations)[2,2], by=0.05))
coordinates(grd) <- c("x", "y")
gridded(grd) <- T

then you run the idw:

idw <- idw(formula =  your_rainfall ~ 1, locations = met_stations, 
           newdata = grd)
idw.output = as.data.frame(idw)  
names(idw.output)[1:3] <- c("long", "lat", "var1.pred") 


ggplot() + geom_tile(data = idw.output, aes(x = long, y = lat, fill = var1.pred))
  • If you look at the OP's code they are using the raster::interpolate function to make an estimate to a prediction grid. This is accomplishing exactly the same thing that you are outlining, in a simpler way using an existing raster extent and resolution. You really did not address the underlying issue at hand, which is not estimation to an array. The error is indicating an error in the formula call to gstat. – Jeffrey Evans Nov 14 '17 at 17:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.