17

I want to perform a line density analysis for my line data, based on an attribute.

The heatmap option (kernel density analysis) seems to be only possible for point data.

I have tried to convert my line data into points, and create a heatmap, but the attributes for the weight were lost on the conversion.

This is a line density analysis performed with ArcGIS. I was looking to do the same with QGIS.

enter image description here

Is it possible to perform a similar line density analysis in QGIS? and this is the data on QGIS

Improve this question
0

5 Answers 5

Reset to default
10

Since QGIS 3.12 there is an algorithm called Line Density listed in the interpolation algorithm group in the Processing Toolbox (accessible via Processing -> Toolbox).

It is a direct port of the ArcGIS Line Density algorithm.

In this example I turn track lines of marine traffic

track lines

into a heat map

marine traffic heat map

Improve this answer
7

I can't find a way to do this in QGIS via the DB Manager / SQLite functions, but in PostGIS, we do the following:

select
(ST_DumpPoints(l.geom)).geom
from linedata as l

This, exported to a new table, allows the points of the lines (for our case from PGRouting students to their school) so we can visualize the hotspot of routes the students 'could' take:

enter image description here

Improve this answer
7

Using a similar approach to the answer from @DPSSpatial, I would use the QChainage plugin to create evenly spaced points along each line.

This plugin generates points along the line, but it doesn't copy across any attributes over from the line.

Here I'm assuming WGS84 (4326, Degrees Lat/Lon). You probably want to do this using an appropriate projection in meters.

  • use QChainage to generate points along the line, e.g. every 0.05 degrees
  • Buffer this layer e.g. by 0.001 degrees.
  • Do a Spatial Join (Vector > Data Management Tools > Join Attributes by Location) (using "Intersects"), from the buffered points to your original line layer. This will let the points pick up the influence value from the lines beneath them.
  • Convert this layer back to points using Vector > Geometry Tools > Polygon Centroids

Finally, I would use the standard heatmap renderer. Make sure you use the weight points by setting and use the influence field.

You should now have something like this. The numbers along each line show the "influence" value.

enter image description here

Improve this answer
7

In fact, it is possible to do line density in QGIS using the standard menus sending parameters to the GDAL executables.

Using gdal_rasterize using the -burn -add switch would create a raster of densities of lines or points, with no search radius as kernel densities and similar. Just how many points or how much lines within a raster cell.

gdal_rasterize is found in the menu Raster -> Conversion -> Rasterize. In the last text field, you press the pencil button and alter the text to include "-burn 1 -add" similar to below. Naturally, you change the -tr (resolution) parameters and the input and output name according to your situation.

gdal_rasterize -burn 1 -add -tr 200.0 200.0 -l lines C:/lines.shp C:/out_lines.tif

Jukka Rahkonen gave me the hints on the gdal_dev mailing list. I've been looking for this for literally years.

Concerning the raster cell values for rasterized lines, I am unsure of what the values represent. Apparently it is not number of lines or lenghts of lines.

Line density with gdal_rasterize

Improve this answer
1

It is possible to perform a heat map from a vector line using the following approach

Step 1: create a Numpy array of zeros

Step 2: from line vector create nested buffers vectors, the first 2000 meters, 1980 meters, and so on. (The step must be selected by the user)

Step 3: Next each buffer vector is transformed to raster, with a value of 0.1 (the value must be selected by the user)

Step 4: Each raster is added to a Zero Numpy array, beginning for the outer (2000 meters buffer)

The Script generates 100 buffers (2000 to 0 meters) from the vector line and generates the following raster. For this work, I used a lineal ponderation function, since in each iteration the script sums 0.1 to the original raster. Heat map raster generated from the code using a lineal ponderation function

from qgis.PyQt.QtCore import QVariant
from PyQt5 import *
from osgeo import gdal, osr
import processing
from qgis.core import *
import numpy as np
gdal.UseExceptions()

## ***************************************************************************************** ##
##              *** Function of the SCRIPT ***
# The aim of this script is to create a Heat Map from a line vector
#
##
## INPUT: - Line Vector ShapeFile
##        - Buffer distance: Must be expressed in meters and must be multiple of raster resolution  
##        - Cte: value to calculate heat map values  
##
## OUTPUT: - Raster: raster with heat map with values 0 for pixels out of buffer's limit
##                   
##
## **************************************************************************************** ##
#
##  **INPUTH PATHS**  ##
RasterPath = 'C:/Users/Marcos/Documents/Pyqgis/Aves Argentinas/Raster'

## **UPLOAD REFERENCE RASTER** ##
qfdr = QFileDialog()
title='RASTER Open'
path ='C:/Users/Marcos/Documents/Pyqgis/Aves Argentinas'
fRast = QFileDialog.getOpenFileNames(qfdr, title, path)
print("Este es F:", fRast)
ListaRast = fRast[0]
PathRast = str(ListaRast[0])
#PathRast = str(PathRast)
ds_Raster = gdal.Open(PathRast)
if ds_Raster is None:
    print('Could not open', + ds_Raster)
    sys.exit(1)
    
ds_RasterRows = ds_Raster.RasterXSize
print("Este es el Row", ds_RasterRows)
ds_RasterCols = ds_Raster.RasterYSize
ds_RasterBands = ds_Raster.RasterCount
# Raster Georreference #
geotransform = ds_Raster.GetGeoTransform()
wkt = ds_Raster.GetProjection()
ZeroRaster = ds_Raster.GetRasterBand(1).ReadAsArray()
width,height = ZeroRaster.shape
ZeroRaster = ZeroRaster * 0 
print("Este es 3000 3000", ZeroRaster[3000,3000])

#
#
## **UPLOAD LINE VECTOR** ##
qfdv = QFileDialog()
title='VECTOR Open'
path ='C:/Users/Marcos/Documents/Pyqgis/Aves Argentinas/Red de Caminos Digitalizada UTM 20S.shp'
fVect = QFileDialog.getOpenFileNames(qfdv, title, path)
print("Este es F:", fVect)
ListaVect = fVect[0]
PathVect = str(ListaVect[0])
BufferPath = 'C:/Users/Marcos/Documents/Pyqgis/Aves Argentinas/Buffer' + '.shp'
print("Este es el Path:", PathVect);
#
n = 1
Cte = 1/100
## BEGING ITERATION: MAKIE BUFFERS VECTOR ##
while n < 101:
    Bondad = Cte
    RadioBuffer = 2000-( 2000*((n-1)/100))
    ## **BUFFER GENERATION** ##
    processing.run("native:buffer", {'INPUT': 'C:/Users/Marcos/Documents/Pyqgis/Aves Argentinas/Red de Caminos Digitalizada UTM 20S.shp',
                'DISTANCE': RadioBuffer,
                'SEGMENTS': 10,
                'DISSOLVE': True,
                'END_CAP_STYLE': 0,
                'JOIN_STYLE': 0,
                'MITER_LIMIT': 10,
                'OUTPUT': BufferPath})
    #
    BufferLyr = QgsVectorLayer(BufferPath, "", "ogr")
    vpr = BufferLyr.dataProvider()
    vpr.addAttributes([QgsField("Value", QVariant.Double)])
    BufferLyr.updateFields()
    #
    fid = 0
    attrs = { 0: "Bufer", 1:0, 2:3888, 3:Bondad }
    BufferLyr.dataProvider().changeAttributeValues({ fid : attrs })
    #
    ##   **RASTER TO VECTOR**   ##
    processing.run("grass7:v.to.rast", {'input':BufferPath,
                    'type':[3],
                    'where':'',
                    'use':0,
                    'attribute_column':'Value',
                    'rgb_column':None,
                    'label_column':None,
                    'value':1,
                    'memory':300,
                    'output':'C:/Users/Marcos/Documents/Pyqgis/Aves Argentinas/Raster/Buffer'+ str(RadioBuffer) +'.tiff',
                   'GRASS_REGION_PARAMETER':'452836.47857953,592833.71854762,6572938.39776585,6713618.64114623 [EPSG:32720]',
                    'GRASS_REGION_CELLSIZE_PARAMETER':19.99960570972714,
                    'GRASS_RASTER_FORMAT_OPT':'',
                    'GRASS_RASTER_FORMAT_META':'',
                    'GRASS_SNAP_TOLERANCE_PARAMETER':-1,
                    'GRASS_MIN_AREA_PARAMETER':0.0001})
    BufferLyr = None
    ##  **DIFFERENCE  ZERORASTER - BUFFER**  ##
    ds_Buffer = gdal.Open('C:/Users/Marcos/Documents/Pyqgis/Aves Argentinas/Raster/Buffer'+ str(RadioBuffer) +'.tiff')
    geotransform = ds_Buffer.GetGeoTransform()
    wkt = ds_Buffer.GetProjection()
    Buffer = ds_Buffer.GetRasterBand(1).ReadAsArray()
    Buffer[np.isnan(Buffer)] = 0
    ZeroRaster[np.isnan(ZeroRaster)] = 0
    ZeroRaster = ZeroRaster + Buffer
    ds_Buffer = None
    Buffer = None
    n += 1

##  **SAVUNG RASTER RESULT**  ##
driver = gdal.GetDriverByName("GTiff")
output_file = "C:/Users/Marcos/Documents/Pyqgis/Aves Argentinas/Raster/Resultado.tiff"

dst_ds = driver.Create(output_file,
                       height,
                       width,
                       1,
                       gdal.GDT_Float32)
ZeroRaster = np.array(ZeroRaster)



#writting output raster
dst_ds.GetRasterBand(1).WriteArray( ZeroRaster )
#setting nodata value
dst_ds.GetRasterBand(1).SetNoDataValue(-999)
#setting extension of output raster
dst_ds.SetGeoTransform(geotransform)
# setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
#Close output raster dataset

ds = None
dst_ds = None
print("Script END")
Improve this answer

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.