## Introduction

Labeled contour lines are a default visual feature on topographic maps. Geographic Information Systems (GIS) such as ArcGIS and QGIS typically allow to display such rotated labels as part of the contour’s styling properties. However, at times it can be useful to have such labels explicitly as point geometries with attached attributes that contain information about height and rotation. In this post, we will show how such point labels can be created automatically with Python.

## Data

For this tutorial, we are using a freely available digital elevation model (DEM) based on data from the Shuttle Radar Topography Mission (SRTM) that can be obtained from CGIAR. Specifically, we are picking an area covering Switzerland and use QGIS to clip it to a manageable extent, project it to the Swiss coordinate reference system (CRS) CH1903/+LV95 (EPSG:2056) and generate contour lines with an interval of 50 meters. Detailed explanations on how contour lines can be derived from a DEM can be found in other places, for example here and here. Figure 1 shows the result.

## Creating Points along Contour Lines

The first step in achieving this goal is to generate points in equal intervals along each of the contour lines. This has already has been discussed several times elsewhere. However, because it is a prerequisite for rotated point labels, we are recapitulating the relevant steps here as well. First, we need to load the contour lines and define the desired properties of the point labels to be generated:

```
import fiona
# open input contours
in_contours_path = "data/contours.shp"
contours = fiona.open(in_contours_path)
# define output contours
out_point_path = "data/points.shp"
schema = { 'geometry': 'Point', 'properties': {'elev': 'int', 'angle': 'float', 'angle_d': 'float', 'angle_d90': 'float'}}
points = fiona.open(out_point_path, 'w', driver="ESRI Shapefile", schema=schema, crs=contours.crs)
```

The first thing required is to import fiona (line 1). This can either be achieved with PIP or with Anaconda, depending on which package management system you prefer. Next, we are defining the path to our contours shapefile (line 4) and open it (line 5). Furthermore, we define the paths where the point labels should be written to (line 9) and define their schema. The schema holds information about the geometry type and the attributes (line 10). In our case, the geometry type is points, and they will hold attribute values about the elevation (“elev”), the angle in radians (“angle”), the angle in degrees (“angle_d”) and the angle in degrees rotated by 90° (“angle_d90”). The latter is necessary to display the labels at the correct angle in QGIS. All these information are used to create the output shapefile (line 11). Note that the CRS is taken directly from from the contours shapefile.

Next, we will iterate over the contours and create points in equal intervals:

```
from shapely.geometry import shape, mapping
height_offset = 100
point_distance = 1000
height_attribute_name = "ELEV"
for contour in contours:
contour_height = int(contour["properties"][height_attribute_name])
# only write points along contour if height is a multiple of height_offset
if not contour_height % height_offset == 0:
continue
# determine the number of points to place for each contour
contour_geometry = shape(contour["geometry"])
point_count = int(contour_geometry.length / point_distance)
for point_number in range(point_count):
point_1 = contour_geometry.interpolate(float(point_number) / point_count, normalized=True)
points.write({ 'geometry': mapping(point_1), 'properties': {'elev': contour_height, 'angle': None, 'angle_d': None, 'angle_d90': None}})
```

The first step we need to take is to import the *shape *and *mapping *functions of shapely (line 1). *shape *allows us to convert a geometry into a shapely object that allows us to use advanced geospatial functions on them. *mapping *converts those objects back to conventional geometries. Just like fiona, shapely can be installed via PIP or Anaconda.

Then we define three variables (lines 3-5). *height_attribute_name *holds the attribute name that holds the elevation information of the contours. *height_offset *is being used to only generate points for those contour lines whose height is a multiple of *height_offset* to avoid clutter. *point_distance* defines the spacing between the point labels.

Then, we iterate over all the single contours (line 7). For each contour, the elevation data of a contour is accessed via *height_attribute_name *(line 8). However, we skip all of those contours which are not a multiple of *height_offset *(lines 11 and 12). Then, we access the contour geometry and convert it to a shapely object (line 15). Now, we can access its *length *property and devide it by the desired points distance to obtain the number of point labels to be placed on that contour (line 16).

Using this information, we create a loop iterates over the interger values from 0 to *point_count* – 1 in each iteration storing the value in *point_number* (line 18). By dividing *point_number *by *point_count*, we obtain a value between 0 and 1. This represents the relative position along that contour starting from its first vertex to its last vertex. Note that we cast *vertex_number *to float to avoid integer divisions. This value can be fed into the *interpolate* method to convert this relative position to an absolute geographical position with x and y coordinates (line 19). We just need to make sure that the *normalized *argument is set to True. Else, a value between 0 and the actual length of the contour would be required. To give a concrete example, if *float(vertex_number) / vertices_count* evaluates to 0.25 and the length of the contour line is 120m, it will return the point coordinates 30m from the first vertex of the geometry along the contour. The resulting point geometry along with the height of the contour for which the point was created are then written to the output file. Note that we use the *mapping* function of shapely to convert the shape object back to an ordinary geometry and that we define all angles as *None *since we are not calculating them yet (line 20). he result is depicted in Figure 2.

## Calculating Local Contour Line Orientation

Great, the point labels have been placed along the contour lines according to our specification. However, although they store the height values, they are not oriented correctly along the contour lines. Hence, we need to adjust our script a bit:

```
import math
for point_number in range(point_count):
dx, dy, point_1 = calculate_dx_dy(contour_geometry, point_number, point_count, small_number)
# calculate the angles for QGIS
angle = math.atan2(dx, dy)
angle_d = round(math.degrees(angle), 4)
angle_d90 = round((math.degrees(angle) + 90) % 360, 4)
points.write({ 'geometry': mapping(point_1), 'properties': {'elev': contour_height, 'angle': angle, 'angle_d': angle_d, 'angle_d90': angle_d90}})
```

To obtain the orientation along a contour, we need the local change in terms of the x-direction and the y-direction. These information will be provided by the function *calculate_dx_dy* (line 5) that is explained later on. Once, these values are obtained, we can use the atan2() function to calculate the angle of the resulting vector in radians (line 8). Radians can be easily converted to degrees with the math.degrees() function (line 9). However, this value needs to be rotated by 90° degrees to conform to the format required by QGIS (line 10).

Finally, we need to define a very small number which can be used to calculate an offset along a contour line. One way to achieve this is to calculate the inverted length of contour line, which, when used in the *interpolate* method, will always evaluate to 1m for that contour line. Multiplying it by 0.01 will yield 1mm. Hence, this will read:

`small_number = (1 / contour_geometry.length * 0.001)`

Now, we have everything in place to come up with the function that actually calculates the offsets in the x-direction and the y-direction:

```
def calculate_dx_dy(contour_geometry, point_number, point_count, small_number):
# get points along the linestring, one exactly at the interpolated position, the second with an offset according to small_number
point_1 = contour_geometry.interpolate(float(point_number) / point_count, normalized=True)
point_2 = contour_geometry.interpolate((float(point_number) / point_count) + small_number, normalized=True)
# calculate the differences between the two coordinates
dx = point_2.coords[0][0] - point_1.coords[0][0]
dy = point_2.coords[0][1] - point_1.coords[0][1]
return dx, dy, point_1
```

This function first calculates the two points along the contour line, *point_1* at the actual location of the point label (line 3) and *point_2* one 1mm further along the contour line starting from the point label (line 4). The coordinates of the resulting points can be accessed with the *coords* property. Points only hold a single coordinate tuple, which can be accessed with the index 0. To obtain the x coordinate for that tuple, one can use the index 0 and for the y coordinate one can use the index 1. By subtracting the respective x coordinates of the two points and the respective y coordinates, the changes along the contour line at this position can be calculated for both directions (lines 7 and 8). These differences along with the position of the point label are returned by the function (line 10). Executing the full script, hence yields the following result:

The labels are rotated correctly along the local direction of the contour lines. The whole script can be found here:

```
# author: Magnus Heitzler
# website: https://www.corundum-digital.com
import fiona, math
from shapely.geometry import shape, mapping
def calculate_dx_dy(contour_geometry, point_number, point_count, small_number):
# get points along the linestring, one exactly at the interpolated position, the second with an offset according to small_number
point_1 = contour_geometry.interpolate(float(point_number) / point_count, normalized=True)
point_2 = contour_geometry.interpolate((float(point_number) / point_count) + small_number, normalized=True)
# calculate the differences between the two coordinates
dx = point_2.coords[0][0] - point_1.coords[0][0]
dy = point_2.coords[0][1] - point_1.coords[0][1]
return dx, dy, point_1
height_offset = 100
point_distance = 1000
height_attribute_name = "ELEV"
in_contours_path = "data/contours.shp"
contours = fiona.open(in_contours_path)
out_point_path = "data/points.shp"
schema = { 'geometry': 'Point', 'properties': {'elev': 'int', 'angle': 'float', 'angle_d': 'float', 'angle_d90': 'float'}}
points = fiona.open(out_point_path, 'w', driver="ESRI Shapefile", schema=schema, crs=contours.crs)
for contour in contours:
contour_height = int(contour["properties"][height_attribute_name])
# only write points along contour if height is a multiple of height_offset
if not contour_height % height_offset == 0:
continue
# determine the number of points to place for each contour
contour_geometry = shape(contour["geometry"])
point_count = int(contour_geometry.length / point_distance)
# offset will be 1 mm
small_number = (1 / contour_geometry.length * 0.001)
for point_number in range(point_count):
dx, dy, point_1 = calculate_dx_dy(contour_geometry, point_number, point_count, small_number)
# calculate the angles for QGIS
angle = math.atan2(dx, dy)
angle_d = round(math.degrees(angle), 4)
angle_d90 = round((math.degrees(angle) + 90) % 360, 4)
points.write({ 'geometry': mapping(point_1), 'properties': {'elev': contour_height, 'angle': angle, 'angle_d': angle_d, 'angle_d90': angle_d90}})
contours.close()
points.close()
```

## Conclusion

In this tutorial, we have discussed a way to calculate point labels along contour lines with locally adjusted rotation, a property that is often desired when creating cartography products. Along the way, you should now have an idea about how to use the libraries fiona and shapely to modify vector data on a low level (i.e., working with single coordinates) and with high level functions (i.e., the *interpolate* method). The method could be optimized in several ways, for example by adding a small offset along a contour line for the first point label as it often is placed at a discontinuity where the first and the last vertices of a contour line meet, thus looking somewhat off. An example for this issue can be found in Figure 4.

If you have any questions about the workflow, feel free to place a comment below. If you like to stay informed about upcoming content from Corundum Digital, consider subscribing to our newsletter here!