Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

How can I get an elevation profile for a band of terrain?

The highest elevation within 10 km (on each side of the defined line) should be taken into account.

I hope my question is clear. Thank you very much in advance.

share|improve this question
    
Is the line defining your profile a simple straight line, or does it consist of several segments with corners? –  Jake Feb 6 '13 at 8:50
    
The line consist of several segments. But all segment are straight line. :) I mean there are no curves. –  Kara Feb 6 '13 at 8:56
    
Just...as they say..spitballing...but could you buffer the line with a 10km buffer. then select all features within the buffer...then select the highest values? –  Gerald Feb 6 '13 at 9:22
1  
Could you provide an image of what you want to achive? –  Alexandre Neto Feb 6 '13 at 9:23
    
@ Alex: the result I want is a regular Elevation Graph. But with a 10km buffer in order that highest value 10km each side of the selected path is shown on the graph. –  Kara Feb 6 '13 at 10:57

3 Answers 3

up vote 10 down vote accepted

Following on from the comments, here's a version that works with perpendicular line segments. Please use with caution as I haven't tested it thoroughly!

This method is much more clunky than @whuber's answer - partly because I'm not a very good programmer, and partly because the vector processing is a bit of a faff. I hope it'll at least get you started if perpendicular line segments are what you need.

You'll need to have the Shapely, Fiona and Numpy Python packages installed (along with their dependencies) to run this.

#-------------------------------------------------------------------------------
# Name:        perp_lines.py
# Purpose:     Generates multiple profile lines perpendicular to an input line
#
# Author:      JamesS
#
# Created:     13/02/2013
#-------------------------------------------------------------------------------
""" Takes a shapefile containing a single line as input. Generates lines
    perpendicular to the original with the specified length and spacing and
    writes them to a new shapefile.

    The data should be in a projected co-ordinate system.
"""

import numpy as np
from fiona import collection
from shapely.geometry import LineString, MultiLineString

# ##############################################################################
# User input

# Input shapefile. Must be a single, simple line, in projected co-ordinates
in_shp = r'D:\Perp_Lines\Centre_Line.shp'

# The shapefile to which the perpendicular lines will be written
out_shp = r'D:\Perp_Lines\Output.shp'

# Profile spacing. The distance at which to space the perpendicular profiles
# In the same units as the original shapefile (e.g. metres)
spc = 100

# Length of cross-sections to calculate either side of central line
# i.e. the total length will be twice the value entered here.
# In the same co-ordinates as the original shapefile
sect_len = 1000
# ##############################################################################

# Open the shapefile and get the data
source = collection(in_shp, "r")
data = source.next()['geometry']
line = LineString(data['coordinates'])

# Define a schema for the output features. Add a new field called 'Dist'
# to uniquely identify each profile
schema = source.schema.copy()
schema['properties']['Dist'] = 'float'

# Open a new sink for the output features, using the same format driver
# and coordinate reference system as the source.
sink = collection(out_shp, "w", driver=source.driver, schema=schema,
                  crs=source.crs)

# Calculate the number of profiles to generate
n_prof = int(line.length/spc)

# Start iterating along the line
for prof in range(1, n_prof+1):
    # Get the start, mid and end points for this segment
    seg_st = line.interpolate((prof-1)*spc)
    seg_mid = line.interpolate((prof-0.5)*spc)
    seg_end = line.interpolate(prof*spc)

    # Get a displacement vector for this segment
    vec = np.array([[seg_end.x - seg_st.x,], [seg_end.y - seg_st.y,]])

    # Rotate the vector 90 deg clockwise and 90 deg counter clockwise
    rot_anti = np.array([[0, -1], [1, 0]])
    rot_clock = np.array([[0, 1], [-1, 0]])
    vec_anti = np.dot(rot_anti, vec)
    vec_clock = np.dot(rot_clock, vec)

    # Normalise the perpendicular vectors
    len_anti = ((vec_anti**2).sum())**0.5
    vec_anti = vec_anti/len_anti
    len_clock = ((vec_clock**2).sum())**0.5
    vec_clock = vec_clock/len_clock

    # Scale them up to the profile length
    vec_anti = vec_anti*sect_len
    vec_clock = vec_clock*sect_len

    # Calculate displacements from midpoint
    prof_st = (seg_mid.x + float(vec_anti[0]), seg_mid.y + float(vec_anti[1]))
    prof_end = (seg_mid.x + float(vec_clock[0]), seg_mid.y + float(vec_clock[1]))

    # Write to output
    rec = {'geometry':{'type':'LineString', 'coordinates':(prof_st, prof_end)},
           'properties':{'Id':0, 'Dist':(prof-0.5)*spc}}
    sink.write(rec)

# Tidy up
source.close()
sink.close()

The image below shows an example of the output from the script. You feed in a shapefile representing your centre-line, and specify the length of the perpendicular lines and their spacing. The output is a new shapefile containing the red lines in this image, each of which has an associated attribute specifying its distance from the start of the profile.

Example script output

As @whuber has said in the comments, once you've got to this stage the rest is fairly easy. The image below shows another example with the output added to ArcMap.

enter image description here

Use the Feature to Raster tool to convert the perpendicular lines into a categorical raster. Set the raster VALUE to be the Dist field in the output shapefile. Also remember to set the tool Environments so that Extent, Cell size and Snap raster are the same as for your underlying DEM. You should end up with a raster representation of your lines, something like this:

enter image description here

Finally, convert this raster to an integer grid (using the Int tool or the raster calculator), and use it as the input zones for the Zonal Statistics as Table tool. You should end up with an output table like this:

enter image description here

The VALUE field in this table gives the distance from the start of the original profile line. The other columns give various statistics (maximum, mean etc.) for the values in each transect. You can use this table to plot your summary profile.

NB: One obvious problem with this method is that, if your original line is very wiggly, some of the transect lines may overlap. The zonal statistics tools in ArcGIS cannot deal with overlapping zones, so when this happens one of your transect lines will take precedence over the other. This may or may not be a problem for what you're doing.

Good luck!

share|improve this answer
3  
+1 That's a nice start to a great contribution! If you look closely at the second figure, you will notice some shorter transects: they are the ones crossing near bends. This is because your algorithm to compute the transects incorrectly assumes the displacement of each segment will equal spc, but bends shorten the displacements. Instead, you should normalize the transverse direction vector (divide its components by the vector's length) and then multiply that by the desired radius of the transect. –  whuber Feb 13 '13 at 23:37
    
You're quite right - thanks for the feedback @whuber! Hopefully fixed now... –  JamesS Feb 14 '13 at 0:04
    
Dear James, I'll try that thank you very much. This solution suits perfectly. –  Kara Feb 14 '13 at 10:20

The highest elevation within 10 km is the neighborhood maximum value computed with a circular 10 km radius, so just extract a profile of this neighborhood maximum grid along the trajectory.

Example

Here is a hillshaded DEM with a trajectory (black line running from bottom to top):

DEM

This image is approximately 17 by 10 kilometers. I chose a radius of just 1 km rather than 10 km to illustrate the method. Its 1 km buffer is shown outlined in yellow.

The neighborhood maximum of a DEM will always look a little strange, because it will tend to jump in value at points where one maximum (a hill top, perhaps) falls just beyond 10 km and another maximum at a different elevation comes just within 10 km. In particular, hilltops that dominate their surroundings will contribute perfect circles of values centered at the point of local maximum elevation:

Neighborhood max

Darker is higher on this map.

Here is a plot of the profiles of the original DEM (blue) and the neighborhood maximum (Red):

Profiles

It was computed by dividing the trajectory into regularly spaced points at 0.1 km apart (starting at the southern tip), extracting the elevations at those points, and making a joined scatterplot of the resulting triples (distance from beginning, elevation, maximum elevation). The point spacing of 0.1 km was chosen to be substantially smaller than the buffer radius but large enough to make the computation go quickly (it was instantaneous).

share|improve this answer
    
That's not entirely accurate though, right? Instead of a circular buffer around each point, shouldn't an orthogonal line of length 20 km be used to sample the underlying raster? At least that's how I would interpret Kara's requirement of "the highest value within 10km on each side of the line" being taken into account. –  Jake Feb 6 '13 at 19:52
4  
@jake I wouldn't say "inaccurate": you merely offer an alternative interpretation. "On each side of" is a vague term that could use better qualification. I can propose solutions for interpretations like yours; one method uses a zonal maximum. However, it's more complicated and much slower in execution. Why don't we first see what the OP thinks of this simple solution? –  whuber Feb 6 '13 at 19:58
    
Bad choice of words, I shouldn't have used "accurate" -- sorry about that –  Jake Feb 6 '13 at 20:20
1  
Because you know how to use the Profile tool, you're almost done. QGIS has interfaces to GRASS which include neighborhood operations. Just apply the neighborhood maximum operation using r.neighbors and profile its result. –  whuber Feb 7 '13 at 16:34
1  
@JamesS You don't want to do a parallel shift, you want to make each cross-profile perpendicular to the center line. (The parallel shift approach can be implemented exactly as I describe here by using an appropriate long and skinny neighborhood for the neighborhood max calculation.) I'm pretty sure you can find code on this site for constructing sets of equally-spaced perpendicular line segments along a polyline; that's the hard part. Everything else is just a matter of extracting the DEM values along those segments and summarizing them. –  whuber Feb 8 '13 at 16:23

I had the same problem and tried James S' solution, but couldn't get the GDAL to work with Fiona.

Then I discovered the SAGA algorithm "Cross Profiles" in QGIS 2.4, and got exactly the result I wanted and that I presume you are looking for too (see below).

enter image description here

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.