Submitted by dylan on Wed, 2006-12-20 23:45.

**Overview**

A simple experiment illustrating the calculation of elevation profiles along arbitrary line segments using a combination of GRASS and R. Fractal dimension was calculated for each profile by spectral decomposition (FFT), as suggested in (Davis, 2002). Note that the elevation profiles are calculated in a rather non-standard way, as elevation is plotted against **3D** distance along the original line segment. This means that resulting transect **length** calculations are sensitive to the resolution of the raster (or distance between sample points), across a rough surface: i.e. smaller grid spacing will result in a longer effective (3D) transect length.

**References**

- Davis, J.C. Statistics and Data Analysis in Geology John Wiley and Sons, 2002

**Sample profile creation I: **Vector representation of 2 transects performed by an adventurous soil scientist near McCabe Canyon, Pinnacles National Monument.

**Create 2 profiles in GRASS**

# digitize some new lines in an area of interest

v.digit -n map=valley bgcmd="d.rast elev_meters"

v.digit -n map=mtns bgcmd="d.rast elev_meters"

# convert to points

v.to.points -i -v -t in=valley out=valley_pts type=line dmax=50

v.to.points -i -v -t in=mtns out=mtns_pts type=line dmax=50

# set resolution to match DEM

g.region -pa res=4

# extract elevation at each point

v.drape in=valley_pts out=valley_pts_3d type=point rast=elev_4m method=cubic

v.drape in=mtns_pts out=mtns_pts_3d type=point rast=elev_4m method=cubic

# export to text files

v.out.ascii valley_pts_3d > valley_pts_3d.txt

v.out.ascii mtns_pts_3d > mtns_pts_3d.txt

# start R

R

**Sample profile creation II: **Elevation profiles for the paths taken by the soil scientist here.

**Process profile data, and plot in R**

# read in the data, note that there is no header, and the columns are pipe delimited

valley <- read.table('valley_pts_3d.txt', header=F, sep="|", col.names=c('easting','northing','elev','f'))

mtns <- read.table('mtns_pts_3d.txt', header=F, sep="|", col.names=c('easting','northing','elev','f'))

# compute the length minus 1 ; i.e. the number of points - 1

# use for the lagged distance calculation

s.valley <- length(valley$easting) - 1

s.mtns <- length(mtns$easting) - 1

# calculate the lagged distance along each component of the 3D vector

dx.valley <- valley$easting[2:(s.valley+1)] - valley$easting[1:s.valley]

dy.valley <- valley$northing[2:(s.valley+1)] - valley$northing[1:s.valley]

dz.valley <- valley$elev[2:(s.valley+1)] - valley$elev[1:s.valley]

dx.mtns <- mtns$easting[2:(s.mtns+1)] - mtns$easting[1:s.mtns]

dy.mtns <- mtns$northing[2:(s.mtns+1)] - mtns$northing[1:s.mtns]

dz.mtns <- mtns$elev[2:(s.mtns+1)] - mtns$elev[1:s.mtns]

# combine the vector components

# 3D distance formula

c.valley <- sqrt( dx.valley^2 + dy.valley^2 + dz.valley^2)

c.mtns <- sqrt( dx.mtns^2 + dy.mtns^2 + dz.mtns^2)

# create the cumulative distance along line:

dist.valley <- cumsum(c.valley[1:s.valley])

dist.mtns <- cumsum(c.mtns[1:s.mtns])

# plot the results

par(mfrow=c(2,1))

plot(dist.valley, valley$elev[1:s.valley], type='l', xlab='Distance from start (m)', ylab='Elevation (m)', main='Profile')

plot(dist.mtns, mtns$elev[1:s.mtns], type='l', xlab='Distance from start (m)', ylab='Elevation (m)', main='Profile')

**Sample profile creation III: **Spectral decomposition (by FFT) for each profile.

**Compute fractal dimension of each transect by FFT**

#compute the fractal dimension via [1]

s.valley <- spectrum(valley$elev[1:s.valley])

s.mtns <- spectrum(mtns$elev[1:s.mtns])

# compute linear model of log(power) ~ log(1/frequency)

# use the slope estimate in [2]

# first define the model formuals: for plotting and for lm()

fm.valley <- log(s.valley$spec) ~ log(1/s.valley$freq)

fm.mtns <- log(s.mtns$spec) ~ log(1/s.mtns$freq)

# next compute the least-squares regression lines via lm()

lm.valley <- lm(fm.valley)

lm.mtns <- lm(fm.mtns)

# extract the slope of each regression line:

slope.valley <- coef(lm.valley)[2]

slope.mtns <- coef(lm.mtns)[2]

# compute Df - the fractal dimension from [3]

Df.valley <- (5 - slope.valley) / 2

Df.mtns <- (5 - slope.mtns) / 2

# plot the figures

plot(fm.valley, main=paste("Valley Profile\nFractal Dimension: ", round(Df.valley,3), sep=''), xlab='Log(1/frequency)', ylab='Log(Power)' )

abline(lm.valley)

plot(fm.mtns, main=paste("Mountain Profile\nFractal Dimension: ", round(Df.mtns,3), sep=''), xlab='Log(1/frequency)', ylab='Log(Power)')

abline(lm.mtns)