[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
area=readShapePoly('/.../my_area')
# 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
levels(area$data1)
area[area$data1=='catchment',] -> catchment
# set up a fitting grid
grid <- expand.grid(x = seq(674400,730700,by=100), y =
seq(4615100,4744400,by=100))
coordinates(grid) <- c("x", "y")
catch.coor <- overlay(grid,catchment)
str(catch.coor)
catch.grid <- grid[!is.na(catch.coor),]
summary(catch.grid)
plot(catch.grid)
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
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
--
Mauricio
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