[R-sig-Geo] Masking interpolations
Wesley Roberts
wroberts at csir.co.za
Mon Feb 16 10:58:09 CET 2009
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/
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