**Home**»

**Software**»

**GRASS GIS: raster, vector, and imagery analysis**»

**Vector Operations**» Raster profile along arbitrary line segments

# Raster profile along arbitrary line segments

**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 lengthcalculations 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

**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
```

```
# 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')
```

```
#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)
```

## Software

- General Purpose Programming with Scripting Languages
- LaTeX Tips and Tricks
- PostGIS: Spatially enabled Relational Database Sytem
- PROJ: forward and reverse geographic projections
- GDAL and OGR: geodata conversion and re-projection tools
- R: advanced statistical package
- GRASS GIS: raster, vector, and imagery analysis
- Generic Mapping Tools: high quality map production