[R] LIDAR Problem in R (THANKS for HELP)

[Dylan Beaudette]

On Tuesday 05 August 2008, Alessandro wrote:
> Hi All,
>
>
>
> I am a PhD student in forestry science and I am working with LiDAR data set
> (huge data set). I am a brand-new in R and geostatistic (SORRY, my
> background it’s in forestry) but I wish improve my skill for improve
> myself. I wish to develop a methodology to processing a large data-set of
> points (typical in LiDAR) but there is a problem with memory. I had created
> a subsample data-base but the semivariogram is periodic shape and I am not
> to able to try a fit the model. This is a maximum of two weeks of work (day
> bay day) SORRY. Is there a geostatistical user I am very happy to listen
> his suggests. Data format is X, Y and Z (height to create the DEM) in txt
> format
>

Hi,

Remeber, R is memory-bound. Why exactly are you trying to develop a method 
that will inherently demand support for massive datasets (LiDAR) within a 
framework designed to work with samples (R)?

There are a wealth of gridding techniques suitable for LiDAR data that are not 
going to fail on massive datasets. A couple come to mind:

1. the blockmean / blockmode functions in GMT
2. r.in.xyz in GRASS
3. RST interpolation
4. natural neighbor interpolation
5 ... lots more ...

i.e. why do you want to use something as memory/processor/model dependent as 
kriging to grid such a massive and regularly-spaced set of data?

I would suggest that the time spent trying to accomplish what you propose 
might be better spent on thinking about the underlying purpose of the 
experiment. Please don't take my comments as anything other that one 
researcher trying to save another researcher's valuable time and effort.

Cheers,

Dylan


> I have this questions:
>
>
>
> 1.       After the random selection (10000 points) and fit.semivariogram
> model is it possible to use all LiDAR points? Because the new LiDAR power
> is to use huge number of points (X;Y; Z value) to create a very high
> resolution map of DEM and VEGETATION. The problem is the memory, but we can
> use a cluster-linux network to improve the capacity of R
>
>
>
> 2.       Is it possible to improve the memory capacity?
>
>
>
> 3.       The data has a trend and the qqplot shows a Gaussian trend. Is it
> possible to normalize the data (i.e. with log)?
>
>
>
> 4.       When I use this R code “subground.uk = krige(log(Z)~X+Y,
> subground, new.grid, v.fit, nmax=40)” to appear an Error massage: Error in
> eval(expr, envir, enclos) : oggetto "X" non trovato
>
>
>
> I send you a report and attach the image to explain better.
>
>
>
> all procedure is write in R-software and to improve in gstat . The number
> of points of GROUND data-set (4x2 km) is 5,459,916.00. The random sub- set
> from original data-set is 10000 (R code is:  > samplerows
> <-sample(1:nrow(testground),size=10000,replace=FALSE) > subground
> <-testground[samplerows,])
>
>
>
> 1.       Fig 1 - Original data-set Histogram: show two populations;
>
> 2.       Fig 2 – original data-set density plot: show again two populations
> of data;
>
> 3.        Fig 3 –Original data-set Boxplot: show there aren’t outlayers un
> the data-set (the classification with terrascan is good and therefore there
> isn’t a problem with original data)
>
> 4.       Fig 4  - plot random data-set in the space
>
> 5.        Fig 5 – ordinary kriging: show a trend in the space (hypothesis:
> the points are very close in the space)
>
> 6.       Fig 6 –de-trend dataset with:  v <- variogram (log(Z)~X+Y,
> subground, cutoff=1800, width=100))
>
> 7.       Fig 7 – map of semi-variogram: show an anisotropy in the space (0°
> is Nord= 135° major radius 45° minus radius)
>
> 8.       Fig 8- semi-variogram with anisotropy (0°, 45°, 90°, 135°).  Best
> shape is from 135°
>
> 9.       Fig 9-  semi-variogram fit with Gaussian Model. R code is:
> > v = variogram(Z~X+Y, subground, cutoff=1800, width=200, alpha=c(135))
> >
> > v.fit = fit.variogram(v, vgm(psill = 1, model="Gau", range=1800, nugget=
>
> 0, anis=c(135, 0.5)))
>
>
>
> R code:
>
>
>
> testground2 <-
> read.table(file="c:/work_LIDAR_USA/R_kriging_new_set/ground_26841492694149_
>x yz.txt", header=T)
>
> class (testground2)
>
> coordinates (testground2)=~X+Y # this makes testground a
> SpatialPointsDataFrame
>
> class (testground2)
>
> str(as.data.frame(testground))
>
>
>
> hist(testground$Z,nclass=20) #this makes a histogram
>
> plot(density(testground$Z)) #this makes a plot density
>
> boxplot(testground$Z)#this makes a boxplot
>
>
>
> samplerows<-sample(1:nrow(testground),size=10000,replace=FALSE) #select n.
> points from all data-base
>
> subground <-testground[samplerows,]
>
> hist(subground$Z,nclass=20) #this makes a histogram
>
> plot(density(subground$Z)) #this makes a plot density
>
> boxplot(subground$Z)#this makes a boxplot
>
> spplot(subground["Z"], col.regions=bpy.colors(), at = seq(850,1170,10))
>
>
>
> library(maptools)
>
> library(gstat)
>
> plot(variogram(Z~1, subground)) #Ordinary Kriging (without detrend)
>
> # if there is a trend we must use a detrend fuction Z~X+Y
>
> x11(); plot(variogram(log(Z)~X+Y, subground, cutoff=1800, width=80))
> #Universal Kriging (with detrend)
>
> x11(); plot(variogram(log(Z)~X+Y, subground, cutoff=1800, width=80, map=T))
>
> x11(); plot(variogram(log(Z)~X+Y, subground, cutoff=1800, width=80,
> alpha=c(0, 45, 90, 135)))
>
> v = variogram(log(Z)~X+Y, subground, cutoff=1800, width=80, alpha=c(135,
> 45))
>
> v.fit = fit.variogram(v, vgm(psill = 1, model="Gau", range=1800, nugget= 0,
> anis=c(135, 0.5)))
>
> plot(v, v.fit, plot.nu=F, pch="+")
>
> # create the new grid
>
> new.grid <- spsample(subground, type="regular", cellsize=c(1,1))
>
> gridded(new.grid) <- TRUE
>
> fullgrid(new.grid) <- TRUE
>
> new.grid at grid
>
> #using Universal Kriging
>
> subground.uk = krige(log(Z)~X+Y, subground, new.grid, v.fit, nmax=40)
> #ERROR



-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341