6

Does anyone have any suggestions for burning a road network into a DEM? I have a Python script written up for this executing this in ArcMap, but I would love to hear if anyone knows of any freely available toolboxes that might have more complicated options/settings (such as removing a greater amount of elevation or a wider path due to the road type - i.e. different "burning options for small country roads versus wide city boulevards).

I'm planning on using this national Tiger/Line dataset: https://www.census.gov/geo/maps-data/data/tiger-line.html

If anyone has any general comments on using this dataset, or whether there may be a better one, I'm all ears. For instance, roads in this dataset appear like they are are broken down into different degrees of resolution - all roads at the county level shapefiles, Primary and Secondary Roads at the state level shapefiles, and Primary Roads at the national level shapefile. I'm thinking to get the coverage I want I would have to merge all that county level data together for the whole US.

  • 3
    Just be careful with tunnels. I once burned a very large valley into a DEM... – Richard Law Jan 14 '15 at 20:17
  • Ah jeez. That's going to be extremely difficult to deal with. Hope that's found in the attribute data. – traggatmot Jan 14 '15 at 20:34
  • @traggatmot Can I just ask why you would want to lower the elevations of a DEM along a road? I can understand burning streams into a DEM (though as you know I'm not keen on it) but I've never heard of road burning before. Just curious. – WhiteboxDev Jan 14 '15 at 21:17
  • AND BRIDGES! Roads go over watercourses on bridges which can cause massive hdrological issues when you raise. It might be more appropriate to do a CON(Roads < DEM, Roads ,DEM) than Mosaic to avoid raising the DEM to the bridge surface. – Michael Stimson Jan 14 '15 at 21:17
  • @WhiteboxDev, I'm looking to force surface flowpaths along road networks in areas with good road attribute information that may not have captured the presence of a road. It's just something I'm looking into testing out on a few watersheds before possibly constructing a national dataset. It's also something I've been told to do, whether or not it makes absolute sense. :( – traggatmot Jan 14 '15 at 21:21
2

The aspect of the question addressed by this answer concerns an efficient way to burn variable-width buffers into a DEM. Although it obviously could be done by extracting each road type, buffering it, and merging the resulting datasets, there's a better way.

The immediate objective is to create a 0-1 indicator grid of where to burn the DEM. After that, the calculations are straightforward and efficient.

To illustrate the idea, suppose there are just two road types, "primary" and "secondary," say, and you wish (a) to indicate those cells through which either a primary or a secondary road passes and (b) also to indicate those cells adjacent to a primary road. This is, in effect, a variable buffer whose typical width is about 1/2 cellsize for secondary roads and 3/2 cellsize for primary roads.

The solution consists of two parts. First, convert the roads data into raster format. Represent primary roads with values of 9, secondary roads with values of 1, and everything else with values of 0. (The origins of these values will become apparent as we go on.) Wherever two roads cross, use the larger of the two values.

Second, compute a weighted focal sum of this raster. (Weighted focal sums can be computed incredibly quickly even for large neighborhoods. As computation and parallel processing improve in the future, these will remain among the fastest possible operations on rasters.) The weights will be defined over a 3 by 3 neighborhood as given by this array:

1 1 1
1 9 1
1 1 1

Select all focal sums of 9 or greater: this is the variable buffer.

The effect, as you can readily check, is the following:

  • Any central cells with a value of 1 or greater will contribute at least 9 times their value, thereby ending in the output. Thus, all cells through which any road passes will be included.

  • Any neighborhoods that include any cell of value 9 will have focal sums of at least 9 and also end in the output. Thus, all cells adjacent to any primary road will be included.

  • However, when a cell is adjacent to (but does not cover) only secondary roads, the focal sum cannot exceed 1 * (1 + 1 + ... + 1) = 8. The first "1" is the value of a secondary road and the sum "(1 + 1 + ... + 1)" is the sum of all neighborhood weights that do not include the central square. Thus, such cells will not be included.

This is exactly as desired. You can see where the central weight of 9 and the primary value of 9 came from: it had to exceed the sum of weights in all edge and corner cells, which were arbitrarily given the value 1. Any value larger than 8 would have worked.

As another example, let there be three levels of roads: primary, secondary, and tertiary. Suppose you want a 5/2 cell buffer of the primary, a 3/2 cell buffer of the secondary, and a 1/2 cell buffer of the tertiary roads. A reasonable neighborhood, and its weights, is

0   1   1   1   0
1 109 109 109   1
1 109 981 109   1
1 109 109 109   1
0   1   1   1   0

The previous example was multiplied by 109 = 12*9 + 1 and bordered by 1's and 0's (the 0's are too far from the center to be of interest). The values to give to the roads are now

tertiary:    1
secondary:   9 =  9*1 = (8)*1 + 1
primary:   981 = 12*9 + 8*(12*9+1) + 1

This time, compare the weighted focal sum to 981. As before, you can see that the places where the sum equals or exceeds 981 are precisely those cells that are either (a) on any kind of road (because the central weight is 981) or (b) next to a primary or secondary road (because 9 times any single one of the middle weights is 981) or (c) within two cells of a primary road (because 981 times any of the nonzero weights is at least 981). No combination of secondary roads along the border with tertiary roads not in the center can exceed 980 = 8*109 + 12*9.

The same technique can be extended to apply to long tiers of road types and arbitrary buffer radii, provided only that larger radii correspond to the higher classes of roads (which more or less will define the road class in the first place). In this fashion a potentially large number of extract-buffer-recombine operations are replaced by a single focal sum-comparison operation. This could be advantageous when working at a national scale: for the variable buffers to make much sense, the cellsize would be on the order of tens of meters, leading to a raster of tens or hundreds of trillions of cells (in the continental US, anyway). Computational efficiency will matter.

  • I do currently use an indicator grid (0/1) for all road types. I don't currently buffer the roads, and I am confused by the width of the buffer for primary/secondary you mention in the 1st example - 1/2 vs 3/2 - but I am still digesting your answer so it may come to me... – traggatmot Jan 14 '15 at 21:43
  • This answer responds to your request to "remove a greater amount of elevation or a wider path due to the road type." Upon specifying how wide each road type should be, you have to convert the width to cell sizes. For instance, with a cellsize of 10m, if you want widths of 10m, 30m, and 65m, they will convert to 1, 3, and 6.5 cellsizes. Since these are widths, the corresponding buffer radii are half that: 1/2, 3/2, and 6.5/2 cellsizes. You would likely round the latter up, giving 0.5, 1.5, and 3.5 cellsizes. This would require a 7 by 7 neighborhood. – whuber Jan 14 '15 at 23:31
5

I would first assign an elevation to your road network by creating a new field in your vector data, and populating it with the desired elevation value of your roads.

Then, convert it to raster using your existing DEM as the 'master raster' - ie. use the same cell size, extent, snap raster, etc.

Now you have your road network in raster form, with the values as the elevation you set earlier.

Next, use the Mosaic tool, which has an option for 'mosaic type'. This basically says, in the list of raster you're going to mosaic (ie. the layers in the dialogue box), which raster do you want to set the value of the cell being written.

So - if you put your Road raster first in the list, then your DEM second, choose the 'first' option in the 'mosaic type', and since all your rasters are the same cell size, extent, etc., you will burn your road raster into your DEM!

  • With this process you will want to consider road width, as roads will vary according to the number of lanes. Unless you have a high resolution DEM, "burning" roads into a DEM could create more issues that its worth. What resolution DEM do have access to? – Ryan Garnett Jan 14 '15 at 21:22
  • My Python script currently converts the road data into a raster using the protocol you describe, and then lowers elevation in any overlapping cells by a fixed amount, which I may vary based on the road attribute data. – traggatmot Jan 14 '15 at 21:23
  • @RyanGarnett, I have 30m and 10m, and i'm considering only doing it for the 10m. – traggatmot Jan 14 '15 at 21:34
  • @RyanGarnett good point! – DPSSpatial Jan 14 '15 at 21:36
  • @taggatmot a 10m DEM will not allow for single lane roads. Just keep that in mind when you are doing your "burning" – Ryan Garnett Jan 14 '15 at 23:29
0

Here are some alternative approaches to the same problem:

The road enforcement algorithm (REA) presented in this paper manipulates flow direction matrices alongside linear landscape features (i.e., roads) by converging the flow patterns towards depressions.

http://www.spatialhydrology.net/index.php/JOSH/article/view/19

And this paper introduces an open source software, GeoNET, that solves the problem only on a small scale:

http://www.sciencedirect.com/science/article/pii/S1364815216301219

Hope this helps.

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.