Skip to content

creating irregularly-shaped grids for kriging

5 messages · Scanlan, Craig, Paul Hiemstra, Pierre Roudier +1 more

#
Hi Craig,

The easiest way is to get a polygon of you study area and than use 
spsample() to create a grid. You could also make a convex hull using the 
chull command and spsample from that. Depending on your configuration of 
points this might be an option. You could try autoKrige from the automap 
pacakge to get an idea of how this last option looks, it uses a convex 
hull to create a new_data object if it is not given.

library(automap)
data(meuse)
coordinates(meuse) = ~x+y

kr = autoKrige(meuse)
plot(kr)

cheers,
Paul
Scanlan, Craig wrote:

  
    
#
Hi Craig,

I have a workaround for this using the getpoly() inout()

2009/7/20, Scanlan, Craig <craig.scanlan at agric.wa.gov.au>:
#
Hi again, and apologises to the list for the posting problem,

I am using a workaround for this using inout() and getpoly() (or
locator()) on agricultural fields that have irregular boundaries:

1- Make your regular grid using expand.grid. Then you got your
rectangular grid under a variable called grd.

2a- If you already have the coordinates of the boundaries of the
polygon, have a look to ?inout from the splancs package to sort the
coordinates of your grid.

bound <- read.csv('myboundary.txt',header=TRUE)
grd <- grd[inout(grd,bound),]

2b- If you want to enter manually those boundaries, plot the
coordinates of your grid and enter the polygon using getpoly() from
the same package :

plot(grd,pch='+')
bound <- getpoly()
grd <- grd[inout(grd,bound),]

Here is a simple example that works on my box:

require(splancs)
x <- seq(0,100,by=2)
y <- seq(0,100,by=5)
grd <- expand.grid(x=x,y=y)
plot(grd,pch='+')
bound <- getpoly()
grd.irregular <- grd[inout(grd,bound),]
plot(grd,pch='+')
points(grd.irregular,pch=16,col='red')

Then you can convert it to a sp object:
names(grd.irregular)[1] <- 'x'
names(grd.irregular)[2] <- 'y'
coordinates(grd.irregular) <- c('x','y')
gridded(grd.irregular) <- TRUE

All the best,

Pierre

2009/7/20, Pierre Roudier <pierre.roudier at gmail.com>:
#
Hi Craig,

If you have a shapefile of your agricultural fields, you can use the
function "readShapePoly" (maptools package) for getting the polygons
of your agricultural fields, and then use the "spsample" command to
create the grid within the polygons, as Paul mentioned. Finally, you
can use the grid generated by spsample as the object with prediction
locations.

Something like:

1) require(maptools)

2) shp <- readShapePoly(yourshapefilename, proj4string=yourp4s)

3) agr.grid <- spsample(shp, type="regular", cellsize=yourcell.size,
offset = c(0.5, 0.5))

4) gridded(agr.grid) <- TRUE

5) x <- krige(yourvalue~1, locations=yourlocations, newdata=agr.grid,
model=yourmodel) (or the autoKrige version of this)


I hope this helps.

Kinds,

Mauricio