[R-sig-Geo] Defining a grid for interpolations ?

Oehler, Friderike (AGPP) Friderike.Oehler at fao.org
Fri May 23 15:12:48 CEST 2008

Dear Mauricio,

with some help of Dr. Rossiter, I once created a grid covering my study area
by importing a vector file delimiting the borders of the area and overlaying
it over the expanded grid:

# import the vector shape-file

# in case that the vector file includes various polygons: extract from these
only the polygon you want to use for your grid, for instance: 
area at data
area[area$data1=='catchment',] -> catchment

# set up a fitting grid
grid <- expand.grid(x = seq(674400,730700,by=100), y =
coordinates(grid) <- c("x", "y")
catch.coor <- overlay(grid,catchment)
catch.grid <- grid[!is.na(catch.coor),]

In case that you have a vector file of your catchment available, or can
create one, maybe this helps. However, the more experienced gstat users or
developers certainly know better advice.

Cheers, Friderike

-----Original Message-----
From: r-sig-geo-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of Mauricio Zambrano
Sent: 23 May 2008 14:13
To: r-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] Defining a grid for interpolations ?

Dear List,

My question is about how to define the grid to be used for the
interpolations, using R 2.7.0 and gstat 0.9-47

I'm working in a catchment of ~3000 km2, with daily rainfall data of several
stations (more than 40), and I would like to interpolate daily values within
all the catchment using ordinary Kriging.

For defining the grid in which the interpolation will be carried out, at the
beginning, I tried with

#setting a grid each 100m vertical and horizontal
dx <- seq(674400,730700,by=100)
dy <- seq(4615100,4744400,by=100)
catch.grid <- expand.grid(dx,dy)

#setting the names of the columns of the grid
names(catch.grid) <- c("x","y")

#setting the coordinates of the grid
coordinates(catch.grid) <- ~x+y

#interpolating with the inverse distance
pp.idw <- krige(PP_DAILY_MEAN_MM~1, meteo_catch_nNA, catch.grid)

#plotting the interpolated values
spplot(pp.idw["var1.pred"], main="Daily Mean PP, [mm]")

but looking at the figure I realized that the interpolations were carried out
considering all the cells within the squared grid, and not only within the

After, I tried reading a raster file (each 100m) and using it as the grid,

DEMM100m <- read.asciigrid("Catch_DEM_c100m")

but the results that I got seems to be wrong, because I got high values where
low values were expected and viceversa. (I really appreciate any help
clarifying me what I'm doing wrong )

Is there any way to define a grid, starting from a shapefile of the catchment
boundaries ?. For example, I would like to define something similar to the
"meuse.grid" dataset.

Thanks in advance and best regards


Linux user #454569 -- Ubuntu user #17469

R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo

More information about the R-sig-Geo mailing list