[R-sig-Geo] Kriging in gstat - "number of rows in data.frame and SpatialPoints don't match" error

[Jon Skoien]

Hi Anna,

The problem seems to be this line:
 > grd <- as(grid, "SpatialPixelsDataFrame")
as the grid doesn't have any data except for the coordinates.
Try instead with:
 > grd <- as(grid, "SpatialPixels")

Cheers,
Jon


On 11/8/2016 3:50 PM, Anna Szyniszewska wrote:
> Hi,
>
> I have done kriging using gstat a few times before, but this time I keep getting an error and I am unsure what is that the krige() function doesn't like - my spatial pixels or point data frames.
> Error in validObject(.Object) :
>   invalid class "SpatialPointsDataFrame" object: number of rows in data.frame and SpatialPoints don't match
>
> I tried to modify my code in several ways, but I was not successful in finding solution. Help of an experienced eye, would be very much appreciated - many thanks in advance.
>
> Below is my fully reconstructible code:
>
> X <- c(32.04, 32.07, 32.11, 32.12, 32.16, 32.17, 32.2, 32.22, 32.24, 32.27,
>        32.3, 32.37, 32.42, 32.45, 32.49, 32.49, 32.51, 32.52, 32.56, 32.57,
>        32.58, 32.59, 32.59, 32.6, 32.61, 32.63, 32.67, 32.7, 32.7, 32.72, 32.74,
>        32.75, 32.78, 32.82, 32.83, 32.88, 32.9, 32.95, 32.96, 32.97, 33.01,
>        33.06, 33.07, 33.08, 33.09, 33.1, 33.13, 33.16, 33.17, 33.17, 33.17,
>        33.17, 33.21, 33.22, 33.25, 33.27, 33.28, 33.29, 33.3, 33.34, 33.36,
>        33.41, 33.41, 33.41, 33.42, 33.43, 33.45, 33.48, 33.49, 33.51, 33.54,
>        33.58, 33.64, 33.64, 33.65, 33.71, 33.74, 33.8, 33.84, 33.9)
> Y <- c(-3.51, -4.79, -4.59, -3.7, -3.02, -4.63, -4.99, -4.42, -3.67, -3.18,
>        -3.27, -3.91, -3.49, -4.82, -4.4, -3.77, -3.73, -4.28, -4.07, -3.81,
>        -3.96, -4.92, -4.42, -3.41, -3.27, -3.01, -4.13, -4.98, -4.58, -3.11,
>        -3.86, -4.45, -4.91, -4.18, -4.33, -4.22, -3.9, -4.96, -3.1, -3.19,
>        -4.21, -3.75, -3.64, -4.06, -3.85, -3.35, -4.24, -4.59, -4.74, -4.44,
>        -4.31, -3.22, -4.2, -3.76, -3.17, -3.03, -3.74, -4.74, -4.16, -3.09,
>        -3.73, -4.23, -3.69, -3.55, -4.33, -4.38, -4.54, -4.45, -3.41, -4.7,
>        -3.65, -4.43, -4.39, -3.11, -4.14, -4.78, -3.61, -3.21, -3.25, -3.36)
> Z <- c(0.87, 0.73, 0.83, 0.2, 0.23, 0.93, 0.1, 0.77, 0.8, 0.27, 0.5, 0.13, 0.33,
>        0.1, 0.13, 0.03, 0.13, 0.87, 0.07, 0.17, 0, 0.17, 0.1, 0.2, 0.13, 0.1,
>        0, 0.07, 0.73, 0.5, 0, 0.17, 0, 0, 0, 0, 0, 0.07, 0.1, 0.07, 0, 0.1,
>        0.03, 0.1, 0, 0.13, 0.13, 0.03, 0, 0, 0, 0.1, 0.23, 0, 0.17, 0.17, 0.03,
>        0, 0.1, 0.4, 0.07, 0.17, 0.13, 0.1, 0.27, 0, 0, 0, 0.03, 0, 0.03, 0.17,
>        0.1, 0, 0.03, 0, 0.03, 0, 0.03, 0)
>
> Data1 <- data.frame(X,Y,Z)
> coordinates(Data1) <- cbind(Data1$X,Data1$Y)
>
> # Create spatial grid for kriging
> x <- rep(seq(33,35,0.01),each=201)
> y <- rep(seq(-5,-3,0.01),times=201)
>
> grid <- data.frame(x,y)
>
> coordinates(grid) <- c("x","y")
> gridded(grid)=TRUE
>
> grd <- as(grid, "SpatialPixelsDataFrame")
>
> # Project SpatialPointsDataFrame and SpatialPixelsDataFrame:
> proj4string(Data1) <-  CRS("+proj=longlat +datum=WGS84")
> proj4string(grd) <-  CRS("+proj=longlat +datum=WGS84")
>
> # Calculate variogram:
> vgm1 <- variogram(Z~1,Data1,width = 10,cutoff=130)
>
> plot(vgm1)
>
> # Set initial variogram function values
> psill = 0.05
> distance = 80
> nugget = 0.01 # constant
>
>
> # Fit and plot variogram model:
> null.vgm <- vgm(psill, "Sph", distance, nugget) # variogram function
> vgm.fit <- fit.variogram(vgm1, model=null.vgm, fit.ranges=TRUE,
>                          fit.method=1) # Raw WLS
> # Plot semivariogram function:
> plot(vgm1,vgm.fit)
> Inc.kriged <- krige(Z~1,Data1,grd,model=vgm.fit)
>
>
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Jon Olav Skøien
Joint Research Centre - European Commission
Institute for Space, Security & Migration
Disaster Risk Management Unit

Via E. Fermi 2749, TP 122,  I-21027 Ispra (VA), ITALY

jon.skoien at jrc.ec.europa.eu
Tel:  +39 0332 789205

Disclaimer: Views expressed in this email are those of the individual 
and do not necessarily represent official views of the European Commission.