0

I have a shapefile with fields (TIME, LATITUDE, LONGITUDE, RAINFALL). I am trying to generate IDW interpolation based on TIME(per day).

library(raster)
library(gstat)
library(sp)
library(sf)
library(rgdal)

shp <- readOGR(('*******\\1.shp'))
extent <- extent(shp)
##### Here I generate grid for n elements of time
grd <- list()
for (i in 1:length(shp$TIME)) {
    grd[[i]] <- expand.grid(x = seq(from = round(extent@xmin),
                                  to = round(extent@xmax),
                                  by = 1),
                          y = seq(from = round(extent@ymin),
                                  to = round(extent@ymax),
                                  by = 1))
  coordinates(grd[[i]]) <- ~x + y
  gridded(grd[[i]]) <- TRUE
  proj4string(grd[[i]]) <- proj4string(shp)
}
####Upto above everything works well

idw.list <- list()
for (i in 1:length(grd)){
  idw.list[[i]] <- gstat::idw(RAINFALL ~ 1, shp, newdata=grd[[i]], idp=2.0) 
}

#### At last the error comes as:
**[1] "+proj=longlat +datum=WGS84 +no_defs"
[1] "+proj=longlat +datum=WGS84 +no_defs"
Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  : 
  var1 : data item in gstat object and newdata have different coordinate reference systems**

Even though the CRS of both objects is same.

Can someone suggest me a solution?

Sonali Sharma
  • 1
  • 1
  • Welcome to SO! Please see [How do I ask a good question?](https://stackoverflow.com/help/how-to-ask) and [How to make a great R reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). As a side note, I'd check out here: https://swilke-geoscience.net/post/spatial_interpolation/. It's super helpful. Once you figure out for 1 layer, it's very easy to apply to the raster stack. – AndrewGB Jan 19 '22 at 06:09
  • 1
    Hello @AndrewGillreath-Brown thank you for your guide. I am editing my question to add with a new query. – Sonali Sharma Jan 19 '22 at 08:48
  • What happens if you just run `a <- gstat::idw(RAINFALL ~ 1, shp, newdata=grd, idp=2.0)`. May we know of `class(grd)`? Just reading the help it seems, that `idw` expect spatial object as a new data, not a list element. Also might be helpful: https://www.geo.fu-berlin.de/en/v/soga/Geodata-analysis/geostatistics/Inverse-Distance-Weighting/IDW-interpolation-of-weather-data/index.html – Grzegorz Sapijaszko Jan 27 '22 at 19:48

0 Answers0