5

I'm trying to generate a profile along a track (bike trail) using the Profile Tool Plugin from QGIS, but it consumes more than 14GB's of RAM so it can't finish it's work.

  • DEM: 200m GeoTIFF
  • Profile-Line (the trail): 600 nodes, 300km length

Both are already generalized.

As you can see in the screenshot below, this task consumes more then 14GB's of RAM. I tried that on different PC's and operating systems. enter image description here

I need an vector graphic (SVG) as output. Is there an alternative to the Profile Tool or a workaround?

2
  • Add gdalinfo report about the structure of the DEM file. Most interesting info is if the GeoTIFF is tiled or not. Do you want a profile line with about 300x5=1500 vertices?
    – user30184
    Jan 7, 2020 at 16:10
  • 2
    Please harmonize your question title and your question text. Do you want to solve it specifically with the Profile Tool or are you looking for alternatives to it. If the latter, please remove the bit about RAM consumption as it would be irrelevant then.
    – inc42
    Jan 11, 2020 at 8:57

6 Answers 6

4
+100

Instead of using a plugin, try using optimised and native tools.

In processing toolbox, try to SAGA > Terrain Analysis > Profiles from lines.

enter image description here

1

Not a solution, more like a workaround if you are using Linux. You could assign more storage to swap. this doesn't make the progress faster, but ensures it doesn't fail.

I had a similar problem while mosaicking xyz grids. My RAM (16GB) was way too small. I gave the swap memory additional 100GB via gparted. While the progress could use the RAM for active calculations, everything that was more than the ram could handle was temporarily stored in the swap area.

1
  • For memory intensive processes, this will make the progress a lot (like, a lot; I cannot stress a lot enough...) slower (i.e. several orders of magnitude, corresponding to disk & memory I/O speed difference). But it does work...
    – geozelot
    Jan 9, 2020 at 13:46
1

I guess the following solution will solve your problem. It uses the QGIS 3.x Atlas Export feature to build terrain profiles together with the corresponding maps:

enter image description hereThere's no need to install any Plugin, because HTML/Javascript (Highcharts) will generate the terrain profiles. Python GDAL is reponsible for fetching the elevation data from the DEM file. I tested the solution with a VRT file, which was a merge of > 5000 TIF files. During my tests I produced terrain profiles from more than 10.000 points as well, without any memory problems.

I use QGIS Project variables for the profile input parameters ("dhmFile","elevationDistance") and for the elevation calculation result ("elevationData") :

enter image description here

All the needed Python Expression Functions are directly stored in the QGIS project "Python Macros".

from qgis.core import qgsfunction,QgsExpressionContextUtils,QgsExpression,QgsProject,QgsCoordinateReferenceSystem,QgsCoordinateTransform
import os,tempfile,re,gdal,osr
from gdalconst import GA_ReadOnly
from qgis.PyQt.QtCore import QTimer,QEventLoop

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  rtnX = geoMatrix[2]
  rtnY = geoMatrix[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / xDist)
  return (pixel, line)

def getElevation(geoMatrix,band,x,y,fname,fileDict):
    if fname[-4:].lower() == '.vrt':
        locInfo = band.GetMetadataItem('Pixel_%s_%s' % world2Pixel(geoMatrix,x,y), 'LocationInfo')
        fname = re.sub(r'.*<File>(.*)</File>.*', r'\1', locInfo)

    if fname not in fileDict.keys():
        dataset = gdal.Open(fname, GA_ReadOnly)
        if dataset:
            matrix = dataset.GetGeoTransform()
            fileDict[fname] = (dataset,matrix)
        else:
            return 0
    else:
        (dataset,matrix) = fileDict[fname]

    bnd = dataset.GetRasterBand(1) # 1-based index
    (pixel,line) = world2Pixel(matrix, x, y)
    elev = (bnd.ReadAsArray(pixel, line, 1, 1))[0][0]
    return elev

@qgsfunction(args=1, group='Custom', usesgeometry=True)
def GetElevationData(values, feature, parent):
    layerId = values[0]
    layer_srs = QgsProject.instance().layerTreeRoot().findLayer(layerId).layer().crs().authid()
    dhm_srs = ''

    id = feature.id()

    dist = float(QgsExpressionContextUtils.projectScope(QgsProject.instance()).variable('elevationDistance'))
    elevationData = QgsExpressionContextUtils.projectScope(QgsProject.instance()).variable('elevationData')
    if not elevationData or (id != int(elevationData.split(',')[0]) or dist != float(elevationData.split(',')[1])):
        dhmFile = QgsExpressionContextUtils.projectScope(QgsProject.instance()).variable('dhmFile')
        ds = gdal.Open(dhmFile, GA_ReadOnly)
        prj = ds.GetProjection()
        srs=osr.SpatialReference(wkt=prj)
        if srs.IsProjected:
            dhm_srs = 'EPSG:' + srs.GetAttrValue('AUTHORITY',1)

        band = ds.GetRasterBand(1)
        geoMatrix = ds.GetGeoTransform()
        fileDict = {}
        geom = feature.geometry()

        if layer_srs != dhm_srs:
            # transform geometry to dhm_srs
            sourceCrs = QgsCoordinateReferenceSystem(layer_srs)
            destCrs = QgsCoordinateReferenceSystem(dhm_srs)
            tr = QgsCoordinateTransform(sourceCrs, destCrs, QgsProject.instance())
            geom.transform(tr)

        data = []

        if dist != 0:
            l = geom.length()
            sum = 0
            distances = []
            while sum+dist < l:
                sum += dist
                distances.append(sum)

            # interpolate points on linestring
            points2d = [(lambda g: (g.x(), g.y()))(geom.interpolate(d).asPoint()) for d in distances]
            vertices = geom.asPolyline()
            start = (vertices[0].x(),vertices[0].y())
            end = (vertices[-1].x(),vertices[-1].y())
            points2d.insert(0,start) # prepend start point

            stat = 0
            for (x,y) in points2d:
                elev = getElevation(geoMatrix,band,x,y,dhmFile,fileDict) # get elevation from VRT
                if stat < geom.length():
                    data.append([stat,round(elev,2)])
                stat = stat + dist

            elev = getElevation(geoMatrix,band,end[0],end[1],dhmFile,fileDict) # get elevation from last vertex
            data.append([round(geom.length(),2),round(elev,2)])
        else:
            # if dist = 0 collect all linestring vertices
            points2d = [(v.x(),v.y()) for v in geom.asPolyline()]

            i = 0
            for (x,y) in points2d:
                elev = getElevation(geoMatrix,band,x,y,dhmFile,fileDict)
                dist = geom.distanceToVertex(i)
                data.append([round(dist,2),round(elev,2)])
                i = i+1

        elevationData = '%s,%s,%s' % (id,dist,str(data))
        QgsExpressionContextUtils.setProjectVariable(QgsProject.instance(),'elevationData',elevationData)
        del ds
        for (dataset,matrix) in fileDict.values():
            del dataset


    return re.search('(\[.+)',elevationData).group(1)


@qgsfunction(args=0, group='Custom')
def wait1000(values, feature, parent):
    loop = QEventLoop()
    QTimer.singleShot(1000,loop.quit)
    loop.exec_()
    return 0


QgsExpression.registerFunction(wait1000)
QgsExpression.registerFunction(GetElevationData)

def openProject():
    pass

def saveProject():
    pass

def closeProject():
    QgsExpression.unregisterFunction('wait1000')
    QgsExpression.unregisterFunction('GetElevationData')
    pass


And here is the HTML/Javascript code for the profile (it needs to be put into a QGIS HTML frame):

<!DOCTYPE html>
<html>
<head>
<style type="text/css">
    #container {
    max-width: 1000px;
    height: 250px;
    margin: 1em auto;
}
</style>

<title>Highcharts Demo</title>
</head>

<body>
<script src="http://code.highcharts.com/highcharts.js"></script>
<script src="http://code.highcharts.com/modules/annotations.js"></script>

<div id="container" style="height: 250px; min-width: 380px;"></div>

<script type="text/javascript">//<![CDATA[

var elevationData = [% GetElevationData( @atlas_layerid ) %];

Highcharts.chart('container', {
    chart: {
        type: 'area'
    },

    title: {
         style: { color: "#2b435d" },
        text: "[%  attribute(  @atlas_feature ,'name' ) %]"
    },

    subtitle: {
         style: { color: "#2b435d" },
        text: 'Total Length: [%round($length)%]m'
    },
    xAxis: {
        labels: {
              style: { color: "#2b435d" },
            format: '{value} m'
        },
        minTickInterval: 250,
        title: {
            text: ' '
        }
    },

    yAxis: {
        startOnTick: true,
        endOnTick: false,
        maxPadding: 0.35,
        title: {
            text: null
        },
        labels: {
              style: { color: "#2b435d" },
            format: '{value} m'
        }
    },

    legend: {
        enabled: false
    },
    plotOptions: {
                area: {
                    fillColor: {
                        linearGradient: {
                            x1: 0,
                            y1: 0,
                            x2: 0,
                            y2: 1
                        },
                        stops: [
                            [0, Highcharts.getOptions().colors[7]],
                            [1, Highcharts.Color(Highcharts.getOptions().colors[2]).setOpacity(0).get('rgba')]
                        ]
                    },
                    marker: {
                        radius: 2
                    },
                    lineWidth: 1,
                    states: {
                        hover: {
                            lineWidth: 1
                        }
                    },
                    threshold: null
                }
            },

    series: [{
        data: elevationData,
        lineColor: Highcharts.getOptions().colors[1],
        color: Highcharts.getOptions().colors[2],
        fillOpacity: 0.5,
        name: 'Elevation',
        marker: {
            enabled: false
        },
        threshold: null
    }]

});
//]]> 

</script>
</body>
</html>


Because of timing problems (QGIS won't wait for our Javascript to finish), we have to add a "wait" Expression function to the HTML frame (i.e. use "Exclude item from export" to attach the function): enter image description here

To display the name of the track/trail, the Atlas feature needs an attribute called "name". To speed up the process, it's best to use the same SRS for the Linestrings as for the DEM.

If you like to output every vertex of your track, simply set "elevationDistance" to 0.

You can copy the distance/z values from "elevationData" if you need it for further processing. But if there are too many values for the text box, you have to fetch the values with Python.

QgsExpressionContextUtils.projectScope(QgsProject.instance()).variable('elevationData')

If you change the "elevationDistance" variable or jump to the next Atlas feature, a new profile will be generated. If you only reopen the Layout, the data from the QGIS Project variable "elevationData" will be read for faster processing.

P.S. don't forget to enable Python Macros in your QGIS project!

1
  • Just a short note to my answer: if the SRS differs between DEM and Linestring, the Python routine tries to determine the EPSG Code of the DEM, which sometimes could fail. In that case you can set the right code using "gdal_edit.py".
    – christoph
    Jan 17, 2020 at 7:53
1

Using ProfileTool plugin, try unchecking the 'Interpolated profile' checkbox under the profile graph before selecting your profile line. When this option is unchecked, the profile is only evaluated at the nodes of the input polyline, otherwise the polyline is interpolated along a greater number of points.

This may help reduce your memory usage, although I am not sure why you are getting such a high value in first place. Feel free to open an issue on https://github.com/PANOimagen/profiletool describing your problem, if you can provide sample files we can take a look at it.

1
  • This would ignore any DEM pixels that lie between track points and I find it highly unlikely that such result would be any use for the questioner.
    – inc42
    Jan 11, 2020 at 8:56
0

Here's a workaround that will let you create a profile graph in different software.

  1. Run the Extract Vertices tool on the bike path to get its vertices as points. The output layer will have new fields, including vertex index (beginning at 0) and distance along original geometry.

  2. Use the 'Sample raster values` tool to copy the elevation values from the DEM to the vertex points.

  3. Export the vertex point layer in CSV or XLSX format.

  4. Create a line graph in your graphing software of choice. Any spreadsheet editing program should be able to create the graph, eg Excel. Use the "distance along original geometry" field for the X axis, and the "elevation" field for the Y axis.

  5. Export the graph in SVG format if your graphing software has that option. Otherwise, export it as an image and convert it to SVG using a different software.

Note: As inc42 pointed out, this method ignores any DEM pixels that lie between track points. So for best results, you will actually want to densify your path line before running this method.

2
  • This would ignore any DEM pixels that lie between track points and I find it highly unlikely that such result would be any use for the questioner.
    – inc42
    Jan 11, 2020 at 8:53
  • 1
    @inc42 That's a good point. I added a note about that to my answer. I still think this answer could be helpful as a workaround. Obviously it's not a full solution, but sometimes you just need a result now. Until someone truly resolves the underlying issue, a workaround might be the best available at the moment.
    – csk
    Jan 14, 2020 at 16:09
-4

A quick Google search leaves me with this possible solution/workaround:

[Extract raster data along shapefile in QGIS

I've quickly tried it in QGIS, which leaves me with the a generated shapefile of points along a line and their extracted height. The sampling distance is a little bit less than 7 meters on a 5 meter input grid, I don't know if this can be changed. Make sure the CRS of the line and the DEM are the same.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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