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