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 catchment.

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

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 catchment.

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

  

Hi Mauricio,

To get a grid based on a shapefile you can use the command "spsample".

library(sp)
library(maptools)
source_poly = readShapePoly("/path/to/poly")
# cellsize is in map units (e.g. km), also see "?spsample"
grd = spsample(source_poly, cellsize = c(10e3,10e3), type = "regular")
gridded(grd) = TRUE # Make it a grid
summary(grd)
# Visualize the grid
spplot(source_poly, sp.layout = list("sp.points", grd))

"grd" can now be used for interpolations.

hth and cheers,
Paul

Mauricio Zambrano wrote:

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 catchment.

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

  

Please note that this "regular" grid is still random, as the first point is
sampled randomly from the first cell. For a fixed, reproducable grid, add
the argument offset = c(0.5, 0.5)
--
Edzer

Paul Hiemstra wrote:

Hi Mauricio,

To get a grid based on a shapefile you can use the command "spsample".

library(sp)
library(maptools)
source_poly = readShapePoly("/path/to/poly")
# cellsize is in map units (e.g. km), also see "?spsample"
grd = spsample(source_poly, cellsize = c(10e3,10e3), type = "regular")
gridded(grd) = TRUE # Make it a grid
summary(grd)
# Visualize the grid
spplot(source_poly, sp.layout = list("sp.points", grd))

"grd" can now be used for interpolations.

hth and cheers,
Paul

Mauricio Zambrano wrote:

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 catchment.

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

Dear Paul and Edzer,

Thanks a lot for your answers.

Before reading the solution proposed by Paul, I had tried with:

# reading the boundary of the catchment
library(maptools)
catchment <- readShapePoly("only_catchment.shp")
catchment.grid <- spsample(catchment, cellsize=100, offset = c(0.5,
0.5)) # reading the boundary of the catchment

 and then I did an IDW interpolation with:

pp.idw <- krige(PP_DAILY_MEAN_MM~1, meteo, catchment.grid)

and it works.

However, after reading your solution, I realised that I didn't use the command:

gridded(catchment.grid)

but I got an interpolation anyway.

Where is performed the interpolation when I don't use a gridded grid (
gridded(catchment.grid) ) ?
  

Mauricio Zambrano wrote:

Dear Paul and Edzer,

Thanks a lot for your answers.

Before reading the solution proposed by Paul, I had tried with:

# reading the boundary of the catchment
library(maptools)
catchment <- readShapePoly("only_catchment.shp")
catchment.grid <- spsample(catchment, cellsize=100, offset = c(0.5,
0.5)) # reading the boundary of the catchment

 and then I did an IDW interpolation with:

pp.idw <- krige(PP_DAILY_MEAN_MM~1, meteo, catchment.grid)

and it works.

However, after reading your solution, I realised that I didn't use the
command:

gridded(catchment.grid)

but I got an interpolation anyway.

Where is performed the interpolation when I don't use a gridded grid (
gridded(catchment.grid) ) ?

it is actually

gridded(catchment.grid) <- TRUE

that changes the object, if it not gridded.

For interpolation, it doesn't matter whether you loop over a set of points,
or a set of points on a grid.

--
Edzer

But if the object (catchment.grid) is not gridded, what points are
used for carrying out the interpolation "pp.idw <-
krige(PP_DAILY_MEAN_MM~1, meteo, catchment.grid "?