[R-sig-Geo] Masking interpolations

[Paul Hiemstra]

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 at coords[,1]))
> y.range <- as.integer(range(a at 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