Dear Alessandro,
If I understand correctly, you have densely collected LiDAR data and would like
to see if the
population/variogram parameters of your target variable (elevation) differ
drastically in different
part of the area. In the case of elevation data and larger areas - I am sure
they do.
In order to do something like this, you would need a tool that can estimate
these parameters locally
i.e. in a moving window (Haas, 1990). This would then allow you to map the
parameters over the whole
area of interest and then do statistical comparison for the purpose of
stratification (Lloyd and
Atkinson, 1998) or localized prediction (Walter et al., 2001).
As far as I know, such analysis is not available in R, but you might try to try
playing with a small
tool Budiman Minasny developed few years ago (Bishop et al. 2006):
http://www.usyd.edu.au/su/agric/acpa/software/digeman.htm
I imagine that something like this could be easily either run from R (it is
distributed as windows
*.exe) or even implemented as a package (Budi works mainly in Matlab). I find
Digeman extremely
interesting (it actually produces maps of variogram parameters!), but have not
much time to test it.
See also: http://www.usyd.edu.au/su/agric/acpa/vesper/vesper.html
all the best,
T. Hengl
http://spatial-analyst.net
References:
Bishop, T.F.A., Minasny, B., McBratney, A.B., 2006. Uncertainty analysis for
soil-terrain models.
International Journal of Geographical Information Science 20, 117-134.
Lloyd, C.D., Atkinson, P.M., 1998. Scale and the spatial structure of landform:
optimizing sampling
strategies with geostatistics. In: Proceedings of the Third International
Conference on
GeoComputation, University of Bristol, UK, 17–19 September 1998. University of
Bristol, Bristol, UK,
16pp.
Haas, T. C., 1990. Kriging and automated semivariogram modelling within a
moving window. Atmospheric
Environment 24A: 1759–1769.
Walter, C., McBratney, A. B., Donuaoui, A., Minasny, B., 2001. Spatial
prediction of topsoil
salinity in the Chelif valley, Algeria, using local ordinary kriging with local
variograms versus
whole-area variogram. Australian Journal of Soil Research 39: 259–272.
-----Original Message-----
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of
Alessandro
Sent: donderdag 7 augustus 2008 19:51
To: [email protected]
Subject: [R-sig-Geo] stationary test
Hi All,
is there un code to do a stationary test in R? because I have a data-base
and I wish to know there is or not a stationary in my data-set.
Alessandro
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
User Guide for DigEMAn (Digital Elevation Model Uncertainty Analysis)
** Theory **
Program for uncertainty analysis & multi windows size
Multi windows based on Jo Wood's thesis Chapter 4
http://www.soi.city.ac.uk/~jwo/phd/
LHS: Latin Hypercube Sampling
McKay, M.D., Beckman, R.J. & Conover, W.J. 1979.
A comparison of three methods for selecting values of input variables
in the analysis of output from a computer code.
Technometrics, 21, 239-245.
Iman, R.L. & Conover, W.J. 1982.
A distribution-free approach to inducing rank correlation among input
variables.
Communications in Statistics Part B-Simulation and Computation, 11,
311-334.
The program works as follows:
-Input DEM (with its uncertainty/ std_dev)
-Determine windows size, 3x3 or 5x5 or 7x7 etc...
nside=1, 3x3 windows cells will look like, centre no.5:
7 8 9
4 5 6
1 2 3
nside=2, 5x5 windows will look like, centre no. 13:
21 22 23 24 25
16 17 18 19 20
11 12 13 14 15
6 7 8 9 10
1 2 3 4 5
- For each point in DEM:
* Extract the neighbourhood
* calculate its local variogram
* Fit a model through the semivariogram data (gamma)
* LHS elevation of the windows cells : elev_1, elev_2....elev_ns1
For elev=1 to ns1
* Fit a quadratic response surface, use 1/gamma as weight
get the parameters, a,b,c,d,e,f & their std. dev
* LHS parameters a,b,c,d,e,f : par_1, par_2, ... , par_ns2
For par=1 to ns2
calculate terrain attributes
next par
next elev
* average the terrain attributes & calc the std dev
The uncertainty of the terrain attributes is the sum of the uncertainty in
elevation & uncertainty in the DTM model
Max. 4000 rows & 4000 columns of dem
** How to use it **
Modify the dem.ini file with your parameter values
the description on left hand side of the '=' sign should not be modified
infile='dem_test.txt'
Input data file, ascii text file with x,y,z,(std_z)
name should be enclosed by quotation mark ' '
output file name
outfile='demout.txt'
dis=10
distance/resolution of dem cell
nside=2
number of half windows size for terrain attributes calculation
nside=1 will give 3x3 neighbourhood, 2 will give 5x5, etc...
use 2 for this instance
ntype=4
type of variogram,
(1) sperical, (2) exponential, (3) gaussian, (4) stable model, (5) Matern, (6)
generalised gauchy
reccomended for dem: (4) or (5) avoid using (3) and (6)
if ntype=0
no variogram modelling, the weight is unity
ntype=-1, weight is inverse proportional to distance from the centre point
ntype=-2, weight is inverse proportional to squared distance from the centre
point
ns1=20
number of LHS for uncertainty in elevation,
for DEM without uncertainty set ns=1, data file will only have 3 columns: x, y,
z
if ns>1 data file should have 4 columns: x,y,z,std_z
ns2=20
number of LHS for uncertainty in quadratic modelling, for terrain attributes
calculation
iprint=0
iprint=1 will print all the realisations of elevation & calculated attributes,
and parameters of the variogram, large file..
iplot=1
attributes to be displayed
1,2, 3, 4.. correspond to the order of the output attributes, e.g. 1= slope, 2=
aspect.. see 'output' below
if iplot=0: no graphic display
pmin=0
min value for the graphics scale of attributes defined above
pmax=10
max value
islop=0
islop=0 use elevation for semivariance calculation,
islop=1 slope is used for semivariance,
islop=2 use std dev of elevation for semivariance calculation.
** Output
The output files will give x, y, z, attributes:
>From Wilson & Gallant
slope : in %, perform arctan(slope/100) to convert into degree
aspect : in degree
profc : profile curvature (in 100/m2)
planc : plan curvature
gradc : gradient/tangential curvature
>From Shary:
meanc : mean curvature
unsph : unsphericity
diffc : differential curvature
minc : minimum curvature
maxc : maximum curvature
>From Jo Wood:
longc : longitudinal curvature
crosc : cross sectional curvature
The std dev of the attributes will be given in the same name as output file
with 'sd_' prefix
The realisations of the attributes will be given in the same name as output
file with 'sim_' prefix
The parameters of variogram will be given in the name with 'vario_' prefix
** Running batch
Execute:
Digeman.exe dem.ini
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo