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([EMAIL PROTECTED],1],[EMAIL > PROTECTED],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)) > [EMAIL PROTECTED] 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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Mauricio Zambrano Sent: dinsdag 27 mei 2008 15:27 To: Edzer Pebesma; Paul Hiemstra Cc: r-sig-geo@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 <[EMAIL PROTECTED]>: > 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@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo