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/ 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 -- This message is subject to the CSIR's copyright terms and conditions, e-mail legal notice, and implemented Open Document Format (ODF) standard. The full disclaimer details can be found at http://www.csir.co.za/disclaimer.html. and is believed to be clean. MailScanner thanks Transtec Computers for their support. _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
