How to fill missing Elevations to empty Contour Lines with PyQGIS

FilledCompound.png

For decades, spatial data was considered only as simple data, it was stored on different formats, with different standards with a sense of storing a location, a place or just creating a map without taking in mind how the data will be useful for future evaluations. From this underrepresentation of spatial data, we get today some extrange cases of vector/raster files without system of reference, spatial data in local coordinates without datum, and sometimes contour lines without elevation or any other attribute.

Elevation contour lines without the elevation attribute is common when we import contour lines from Autocad DXF files, but it has also happened that contour data was stored on a hard disk drive without elevation attribute and years later data is found and there is no one to ask to restore the missing information. This tutorial shows a practical procedure to fill missing elevations on contour lines with the use of PyQGIS on QGIS 3. The procedure uses an intersection line that crosses the contour lines where the base elevation and interval is known. There are some specific instructions to run the script that are well described on the video.

It is strongly recommended to follow all the steps with the data provided on the input part before applying the procedure to your specific case with your own data.

Video

Python code

This is the PyQGIS code used to fill the missing Elevations

from PyQt5.QtCore import QVariant

mapcanvas = iface.mapCanvas()
layers = mapcanvas.layers()

intLineObject = layers[0]
topoLines = layers[1] 
print('Layer Names: '+intLineObject.name()+' '+topoLines.name())

intLine = intLineObject.getFeature(1)
print("Intersection Line: "+intLine.geometry().asWkt())

originPoint = intLine.geometry().vertexAt(0)
print("Origin Point: "+originPoint.asWkt())

topoLines.startEditing()

for elevLine in topoLines.getFeatures():
    intPointObject = elevLine.geometry().intersection(intLine.geometry())
    
    if intPointObject:
        intPoint = intPointObject.asPoint()
        refPoint = QgsPoint(intPoint.x(),intPoint.y())
        distance = originPoint.distance(refPoint)
        print(distance)
        
        topoLines.changeAttributeValue(elevLine.id(), 1, distance)        

topoLines.commitChanges()

##### Segunda Parte ######

topoLines.startEditing()

topoLines.setSubsetString("Dist > 0")

baseElev = 350
interInterval = 50

def get_name(f):
    return f['Dist']

topoFeatures = sorted(topoLines.getFeatures(),key=get_name)

i=0

for elevLine in topoFeatures:
    print(elevLine["Dist"])
    topoLines.changeAttributeValue(elevLine.id(), 0, baseElev+i*interInterval)
    topoLines.changeAttributeValue(elevLine.id(), 1, 0)
    i+=1
    
topoLines.commitChanges()

topoLines.setSubsetString("")

Input data

You can download the input files for this tutorial from this link.

Comment

Saul Montoya

Saul Montoya es Ingeniero Civil graduado de la Pontificia Universidad Católica del Perú en Lima con estudios de postgrado en Manejo e Ingeniería de Recursos Hídricos (Programa WAREM) de la Universidad de Stuttgart con mención en Ingeniería de Aguas Subterráneas y Hidroinformática.

Smiley face

Subscribe to our free e-newsletter for tutorials, articles, webminars, courses and more.