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

[Tomislav Hengl]

FYI

I have just finishing preparing an R script that accompanies my paper in Computers in Geosciences on
finding the grid cell size for various applications.

The script and the data sets can be obtained from:

http://spatial-analyst.net/pixel.php 

This gives a step-by-step guide to estimate suitable grid cell sizes given an input data with its
inherent properties: scale, positional accuracy, inspection density, spatial correlation and
complexity of terrain. The bottom line is that there are objective ways to calculate a suitable cell
size (or at least a range of suitable cell sizes), and that the selection of grid cell size can be
automated. For example, using the akima package we can automate estimation of the grid:

> lbrary(maptools)
> elevations <- readShapePoints("elevations.shp", proj4string=CRS(as.character(NA)))
> library(akima)
> elevations.interp <- interp(x=elevations at coords[,1],y=elevations at coords[,2],z=elevations$VALUE,
extrap=T, linear=F)
> image(elevations.interp, col=bpy.colors(), asp=1)
> library(spatstat)
> dem.interp <- as.SpatialGridDataFrame.im(as.im(elevations.interp))
> dem.interp at grid

But I would rather first fit a variogram and then select the cell size based on the number of point
pairs and given the effective range of spatial autocorrelation.


For more info see also:

Hengl T., 2006. Finding the right pixel size. Computers and Geosciences, 32(9): 1283-1298.
http://dx.doi.org/10.1016/j.cageo.2005.11.008 


Tom Hengl
http://spatial-analyst.net 



-----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: dinsdag 27 mei 2008 15:27
To: Edzer Pebesma; Paul Hiemstra
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Defining a grid for interpolations ?

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) ) ?

Best regards

Mauricio

2008/5/26 Edzer Pebesma <edzer.pebesma at uni-muenster.de>:
> 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
>>>
>>>
>>
>>
>
>



-- 
hzambran

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