Masking interpolations

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.hiemstra at 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 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
    
      

************************************************************************************************************************************************************************************************************

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/35273872 at 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/35273872 at 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/35273872 at 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.hiemstra at 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 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