...in addition:
You could also use AcrGIS, QGis or another GIS package to make a polygon
that delineates the areas in which you want to interpolate. You can then
use spsample() like I showed you in the previous mail. Making this
polygon can also be done in R using spplot and grid.locator, if you are
interested I can send you some example code.
Another option is to calculate the square grid you make using
expand.grid, perform the interpolation and throw away any points that
are above a certain kriging variance (that you have to choose). This
would be a good solution if you need to automatically interpolate a lot
of lidar sets that differ a lot in configuration of the points. it is
not the most elegant and efficient solution though :).
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
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
