[R-sig-Geo] Problem overlying Kriging grid with measurement points

[Roger Bivand]

On Wed, 28 May 2008, Mauricio Zambrano wrote:

> Dear List,
>
> I carried out a IDW interpolation and I got a plot of the
> interpolations values over the grid that I defined.
>
> However, for a better interpretation of the results, I would like to
> plot, over the Kriging grid, the points that were used as data for the
> interpolation, which are stored in a shapefile.
>
> For doing this, I tried with:
> spplot(pp.idw["var1.pred"], main="Daily Mean PP, [mm]", sp.layout =
> list("SpatialPointsDataFrame", meteo_pt[,"PP_DAILY_MEAN_MM"] ) )
>
> but I only get the plot of the interpolated value over the grid, and
> the plot shows the following error message:
>
> "error using packet 1"
> "Missing 'data' argument, no omission value"
>
> Could you give some hint about what I'm doing wrong ?.

sp.layout=list("sp.points", meteo_pt)

See Note in ?spplot

Roger

>
> Additional information is:
>
> class(pp.idw)
> [1] "SpatialPixelsDataFrame"
> attr(,"package")
> [1] "sp"
>
> class(meteo_pt)
> [1] "SpatialPointsDataFrame"
> attr(,"package")
> [1] "sp"
>
>
> And complete tried procedure:
>
> p4s <- CRS("+proj=utm +zone=30 +ellps=intl +units=m +no_defs")
>
> meteo <- read.csv2("Meteorological_by_basin.csv")
>
> # Selecting only those stations that have a NON "NA" value on the
> field "PP_DAILY_MEAN_MM":
> meteo_nna <- meteo[!is.na(meteo[,"PP_DAILY_MEAN_MM"]),]
>
> # Settting the COORDINATES
> coordinates(meteo_nna) <- ~EAST_ED50 + NORTH_ED50
>
> # Projecting the coordinates of the meteorological stations into the
> right system
> proj4string(meteo_nna) = p4s
> #proj4string(meteo_nna) <- CRS("+init=epsg:28992") ; # another way to
> set up the coordinates
>
> # Reading the boundary of the catchment
> library(maptools)
> catchment <- readShapePoly("only_catchment.shp", proj4string=p4s)
>
> # Defining a sampling GRID of 1000mx1000m with regular spacing
> # For avoiding a random grid (sampled randomly from the first cell),
> and getting a fixed,
> # and reproducable grid, it is neccessary to add the argument "offset
> = c(0.5, 0.5)"
> catchment.grid <- spsample(catchment, type="regular", cellsize=1000,
> offset = c(0.5, 0.5))
>
> # Makin possible that the grid can be used in the interpolations:
> gridded(catchment.grid) <- TRUE
>
> # Plotting the grid
> spplot(catchment, sp.layout = list("sp.points", catchment.grid))
>
> library(gstat)
> # Interpolating with the INVERSE DISTANCE WEIGHTED
> pp.idw <- krige(PP_DAILY_MEAN_MM~1, meteo_nna, catchment.grid)
>
> # Plotting the interpolated values
> spplot(pp.idw["var1.pred"], main="Daily Mean PP, [mm]")
>
> # Reading the Meteorologicla Stations whitin the cachment
> meteo_pt <- readShapePoints("meteo_by_basin.shp", proj4string=p4s)
>
> #Attempting to get the overlay:
> spplot(pp.idw["var1.pred"], main="Daily Mean PP, [mm]", sp.layout =
> list("SpatialPointsDataFrame", meteo_pt[,"PP_DAILY_MEAN_MM"] ) )
>
>
> Thanks in advance,
>
> Mauricio
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no