[R-sig-Geo] Masking interpolations

Wesley Roberts wroberts at csir.co.za
Mon Feb 16 11:34:11 CET 2009 

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


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

More information about the R-sig-Geo mailing list