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.
Answer
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.
No comments:
Post a Comment