Hi,
If you want to use the standard convex hull of automap, just omit the
new_data argument. That will stuff 5000 points into the convexhull on a
regular grid. If you want to specify the pixelsize you need to adapt the
create_new_data function like this:
create_new_data = function(obj) {
# Function that creates a new_data object if one is missing
convex_hull = chull(coordinates(obj)[,1],coordinates(obj)[,2])
convex_hull = c(convex_hull, convex_hull[1]) # Close the polygon
d = Polygon(obj[convex_hull,])
# This line is different, in addition I assume that your coordinates
are in meters
new_data = spsample(d, type = "regular", cellsize = 10e-2)
gridded(new_data) = TRUE
return(new_data)
}
# perform the interpolation
kr = autoKrige(I~1, a, create_new_data(a), nmax=100)
cheers,
Paul
Wesley Roberts wrote:
Excellent, that will help alot. I have the newest version of automap and am running an experiment now.
Just to confirm I need to omit the statement 'grd(new_data argument)' which in my case
is 'gridded(grd) <- TRUE'. That way the resolution is defined in x.range <- .....
and y.range <- ..... statements but the entire grid (square) is not defined and thus
interpolation only occurs where data exists?
You need to either omit the grd part from the call to autoKrige or use
the new function I showed above.
Wesley
Wesley Roberts MSc.
Researcher: Earth Observation (Ecosystems)
Natural Resources and the Environment
CSIR
Tel: +27 (21) 888-2490
Fax: +27 (21) 888-2693
"To know the road ahead, ask those coming back."
- Chinese proverb
Paul Hiemstra <p.hiems...@geo.uu.nl> 02/16/09 12:08 PM >>>
Hi,
A convex hull should do the trick, if you have the latest automap
version (>0.5-2) you don't have do to anything to get a convex hull.
Just omit grd (new_data argument) and automap will automatically use a
convex hull.
cheers,
Paul
Wesley Roberts wrote:
Hi Paul,
Many thanks for the reply. I am not sure about the convex hull approach as I am
not sure how to implement it as part of my program. Is the code you wrote below
replacing the following statements? Would I then pass grd to the autoKrige
statement as shown in the original code below?
************************************************************************************************************************************************************************************************************
x.range <- as.integer(range(a...@coords[,1]))
y.range <- as.integer(range(a...@coords[,2]))
grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=0.1),
y=seq(from=y.range[1], to=y.range[2], by=0.1))
coordinates(grd) <-~ x+y
gridded(grd) <- TRUE
************************************************************************************************************************************************************************************************************
I have uploaded some pics of my points and interpolations to flickr.
Link 1 shows the timber compartment with the lidar points overlayed and the bounding polygon used to subset the point data set. You cant really see the irregularity of the points but trust me, they are all over the place. Average distance between points is about 17cm so a 10cm interpolation resolution should be okay.
Link 1
http://www.flickr.com/photos/35273...@n07/3283623165/
Link 2, is the height interpolation using the parameters you suggested on
Friday. As you can see the lack of points outside the polygon results in a type
of edge effect at the boundaries. I can mask the rest out but would prefer to
limit the interpolation to minimize errors at the boundaries.
Link 2
http://www.flickr.com/photos/35273...@n07/3283623941/
Finally link 3 shows what I think is the kriging variance although I cant be
sure. When I import the tiff written by writeGDAL there are three bands (using
GRASS). The first is the interpolated variable, the next two are a mystery to
me. If this is indeed the krig variance then limiting the interpolation based
on kriging variance seems like a good idea? What do you think?
Link 3
http://www.flickr.com/photos/35273...@n07/3284446466/
This is indeed the kriging variance.
Many thanks,
Wesley
Wesley Roberts MSc.
Researcher: Earth Observation (Ecosystems)
Natural Resources and the Environment
CSIR
Tel: +27 (21) 888-2490
Fax: +27 (21) 888-2693
"To know the road ahead, ask those coming back."
- Chinese proverb
Paul Hiemstra <p.hiems...@geo.uu.nl> 02/16/09 10:29 AM >>>
Hi Wesley,
You could take a look at using a convex hull. I'm not sure if this will
fix your problem as we cannot see how exactly your points are irregular.
The latest version on my website (0.5-2) uses a convex hull off the data
if you don't pass a new_data object. You could try this. The function
making the convex hull is:
create_new_data = function(obj) {
# Function that creates a new_data object if one is missing
convex_hull = chull(coordinates(obj)[,1],coordinates(obj)[,2])
convex_hull = c(convex_hull, convex_hull[1]) # Close the polygon
d = Polygon(obj[convex_hull,])
new_data = spsample(d, 5000, type = "regular")
gridded(new_data) = TRUE
return(new_data)
}
If you want to call it directly from the package use
automap:::create_new_data.
cheers,
Paul
Wesley Roberts wrote:
Dear R-sig-geo'ers
I am currently running some interpolations using automap written by Paul Hiemstra. So far my interpolations have been producing suitable results except for one problem. From the code you will see that the boundaries of the spatial grid are determined using the range of the X and Y coordinates creating a square grid. My point data do not cover the entire grid and I would only like to interpolate in areas where data exists otherwise I get a significant edge effect. Is it possible to limit / mask my interpolation to only predict where data exists?
The point data are lidar canopy returns for an irregular shaped timber compartment and number around 10 000 irregular spaced points.
Any help on this matter would be greatly appreciated.
Kind regards,
Wesley
library(automap)
library(gstat)
a <- read.csv("AreaOne_4pts.csv", header=TRUE)
coordinates(a) <-~ x+y
x.range <- as.integer(range(a...@coords[,1]))
y.range <- as.integer(range(a...@coords[,2]))
grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=0.1),
y=seq(from=y.range[1], to=y.range[2], by=0.1))
coordinates(grd) <-~ x+y
gridded(grd) <- TRUE
height = autoKrige(H~1, a, grd, nmax=100)
writeGDAL(height$krige_output, fname="test.tiff", drivername ="GTiff", type =
"Float32")
intensity = autoKrige(I~1, a, grd, nmax=100)
writeGDAL(intensity$krige_output, fname="test.tiff", drivername ="GTiff", type =
"Float32")
Wesley Roberts MSc.
Researcher: Earth Observation (Ecosystems)
Natural Resources and the Environment
CSIR
Tel: +27 (21) 888-2490
Fax: +27 (21) 888-2693
"To know the road ahead, ask those coming back."
- Chinese proverb
--
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone: +31302535773
Fax: +31302531145
http://intamap.geo.uu.nl/~paul
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
