[GIS] How to count the number of vector lines that pass though a given raster cell

Question

I am working on a project (currently in GRASS7 though I'm not against exporting to Arc or R or other software if necessary) where I am generating thousands of shortest cost paths (using the r.walk and r.drain modules) to various points around my map of Papua New Guinea. What I am interested in is figuring out how to take the generated vector lines and create a raster that counts the number of lines that pass through each cell of that raster. This will give me an estimate of which locations on the island are most likely to be located on a shortest path from one point to another.

One way to do it would be to (using python)

1. Create a raster initalized to 0
2. use a python script to run though and rasterize each line individually and
3. add each line to the raster using raster calculator to create a count.

The problem with this approach is that for the tens or even hundreds of thousands of lines I'm using, it takes a very very long time. I would like a way to directly count the number of lines. Something like http://grasswiki.osgeo.org/wiki/Count_points_in_raster_cells
except with lines.

Answer 1

Here's a Python solution that will run 100,000 simple linestrings in 56 seconds. My code could probably more efficient, but it's just a quick hack. This particular solution just generates n random lines, computes the slope and iterates over the column or row, depending on whether the line's slope is >= abs(1). If the slope if infinite it just increments the row. The script loads an existing raster to grab the projection/transform/size and creates a new uint32 raster for the output.

#!/usr/bin/env python

import gdal
import sys
import numpy as np
from random import random
from shapely.geometry import LineString, Point

# Arguments are, input tif, output tif, number of lines (default 100)
infile, outfile = sys.argv[1:3]
if infile == outfile:
    raise Exception('Input and output must differ.')

lineCount = 1000
try:
    lineCount = int(sys.argv[3])
except: pass
if lineCount < 1:
    raise Exception('Line count must be > 0')

# Open the source raster and get some info.
src = gdal.Open(infile, gdal.GA_ReadOnly)
trans = src.GetGeoTransform()
proj = src.GetProjection()
cols = src.RasterXSize
rows = src.RasterYSize
src = None

# Resolution
rx = trans[1]
ry = trans[5] # Negative in UTMN.

# Determine the boundaries.
minx = trans[0]
maxy = trans[3]
maxx = minx + cols * rx
miny = maxy + rows * ry

# Create a grid for counting.
grid = [0] * cols * rows

# Generate some random lines. You'd load your own here.
# For this code to work as-is, you'd have to break them up
# into simple linestrings.
lines = []
for i in range(lineCount):
    x1 = minx + (maxx - minx) * random()
    y1 = miny + (maxy - miny) * random()
    x2 = minx + (maxx - minx) * random()
    y2 = miny + (maxy - miny) * random()
    lines.append(LineString([(x1, y1), (x2, y2)]))

# Infinity
inf = float('Inf')

# Snap to grid so the col/row calcs work correctly
minx = int(minx / rx) * rx
miny = int(miny / -ry) * -ry
maxx = minx + (cols + 1) * rx
maxy = miny + (rows + 1) * -ry

# Iterate over the lines...
for line in lines:
    # Get start and end points.
    x1, y1 = line.coords[0]
    x2, y2 = line.coords[1]
    # Get the slope. Use + or - inf to indicate the direction of a vertical.
    slope = (y2 - y1) / (x2 - x1) if (x2 - x1) != 0. else inf * (1 if (y2 - y1) > 0 else -1)
    # Rasterize the line...
    x, y = (int(x1 / rx) * rx, int(y1 / -ry) * -ry)
    while x <= x2:
        # Compute the row/col indices and increment the grid.
        c = int(((x - minx) / (maxx - minx)) * cols)
        r = int(((y - miny) / (maxy - miny)) * rows)

        # Increment the grid count
        grid[r * cols + c] += 1

        # Compute the value of y, given x or
        if abs(slope) == inf:
            y += -ry * (1 if slope > 0. else -1)
        elif abs(slope) >= 1.:
            x = (1. / slope) * (y - y1) + x1
            y += -ry * (1 if slope > 0. else -1)
        else:
            y = slope * (x - x1) + y1
            x += rx

# Convert grid to np array for writing
grid = np.array(grid).reshape(rows, cols)

# Create a new tiff.
drv = gdal.GetDriverByName('GTiff')
dst = drv.Create(outfile, cols, rows, 1, gdal.GDT_UInt32)
dst.SetGeoTransform(trans)
dst.SetProjection(proj)
band = dst.GetRasterBand(1)
band.WriteArray(grid)

Here's what the output looks like. Red is high count (79), blue is low (1). (Remember, these are just random lines, so they concentrate in the centre.)