6

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.

  • You can't do raster/vector directly. What you can do is either create a polygon from the raster cells (1 cell = 1 polygon) then overlay. Tools exist for raster=>polygon but depending on the size of the raster this might not be the best option. You can create a polygon on the fly using QGIS python but if you don't script that option isn't the best for you. Try clipping the raster out into a few zones, convert to polygons and then intersect and summarize. – Michael Stimson Jul 18 '15 at 0:08
  • @MichaelMiles-Stimson Thanks for the response! I was afraid that would be the answer. As you've guessed, my raster is quite large, but splitting it up is a good idea. I could write a script to break it into 100 different blocks and process each of those individually without too much trouble I think. I've done some basic python scripting before. Could you expand a little on what you mean by creating a polygon "on the fly" in QGIS? Thanks again! – James ZD Jul 19 '15 at 15:25
  • In the QGIS API for python you can create a polygon from known points, from the raster you get upper left (origin) Cell size X and Y, rows and columns then iterate for each row and each column and create a polygon from points which can then be used as a boolean while iterating your line layer with Intersect qgis.org/api/… – Michael Stimson Jul 19 '15 at 21:04
7

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.)

Line-in-cell count raster

1
  1. Convert Raster to Polygons. One cell per polygon.
  2. Use a Spatial Join from the Polygon layer to the to Lines. Much like this How can all polylines within a polygon be selected and then operated on to determine the average value from an attribute field? (do not check the sum box, this will give you count I hope).
  3. Convert the polygons back to raster.

Seems convoluted but I see no reason it will not be rapid unless the raster has hundreds-of-thousands of cells and it is easily conducted in a model or similar. I am sure you could use the same process in QGIS or Grass.

0

If the least cost paths were extracted from a raster map, then I assume they would have line vertexes that fall on the raster cells (unless you used the knight's step but it may still work well enough).

You could then:

  1. run v.patch on all your vector LCPs (I think this would be pretty fast even on a large number of vector maps)
  2. run v.to.points, which will convert all of the line vertices to a map of points, which will fall on your raster cells
  3. run v.out.ascii on those points and pipe it to r.in.xyz method=n to get a count, which is described in the link you gave

All those algorithms are relatively fast in my experience.

Here's a toy example. It's in grass 6 but I imagine it'll still work:

g.region w=0 e=10 s=0 n=10 res=1
r.mapcalc cost_surf=1.0 

r.cost input=cost_surf output=cumulative_cost coordinate="2.0,2.0" stop_coordinate="8.0,8.0" --o
r.drain input=cumulative_cost output=lcp1 voutput=lcp1 coordinate="8.0,8.0" --o

r.cost input=cost_surf output=cumulative_cost coordinate="2.0,8.0" stop_coordinate="8.0,2.0" --o
r.drain input=cumulative_cost output=lcp2 voutput=lcp2 coordinate="8.0,2.0" --o

v.patch input=lcp1,lcp2 output=lcp_all --o

v.to.points -v input=lcp_all output=lcp_all_pts type=line --o
v.out.ascii input=lcp_all_pts output=- format=point | r.in.xyz input=- output=path_cnt method=n

GRASS 6.4.3 (ncspm):~ > r.out.ascii input=path_cnt output=-
north: 10
south: 0
east: 10
west: 0
rows: 10
cols: 10
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 1 0 0 0 0 0 1 0 
0 0 0 1 0 0 0 1 0 0 
0 0 0 0 1 0 1 0 0 0 
0 0 0 0 0 2 0 0 0 0 
0 0 0 0 1 0 1 0 0 0 
0 0 0 1 0 0 0 1 0 0 
0 0 1 0 0 0 0 0 1 0 
0 0 0 0 0 0 0 0 0 0 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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