Line extraction with GRASS GIS

This post presents another way how to use Landsat imagery to map features for OpenStreetMap in an easy way with free and open source software only. After remapping the Nam Ngum 1, Nam Leuk reservoir and other lakes using a similar approach, I wanted to improve the Nam Ngum river. The second largest river in Laos and an important tributary to the Mekong was in OpenStreetMap only mapped as a single line instead of an area, although it's more than hundred metre wide downstream of the Nam Ngum 1 reservoir.

As approach I chose to apply edge detection algorithms to the satellite images, namely the Sobel filter. To get used to the rather new Python scripting in GRASS GIS, I wrote a Python script to preprocess the images. The Python scripting ability is a great addition to GRASS GIS, since it facilitates considerably the entry point to scripting.

First all necessary modules needed to be imported.

#!/usr/bin/env python

from grass.script import core as gcore
from grass.script import raster as graster
from string import join
import sys

In the main function I defined a list of regions to reduce the computational area and thus processing time.

def main():

    # Create a list of regions along the river to reduce the computation time
    regions = []
    regions.append({'w': 279850, 'e': 300648, 's': 2003246, 'n': 2017059})
    regions.append({'w': 260146, 'e': 280944, 's': 2001343, 'n': 2015156})
    regions.append({'w': 242189, 'e': 262987, 's': 2005578, 'n': 2019392})
    regions.append({'w': 234104, 'e': 243982, 's': 2017410, 'n': 2036064})
    regions.append({'w': 231550, 'e': 245964, 's': 2033931, 'n': 2052559})

As in my previous approach I used again the near infrared band from Landsat images to distinct water bodies from land.

    
    # The input map
    inputMap = "20011130_B.4"
    # Make sure the region especially the resolution is correctly set
    gcore.run_command("g.region", rast=inputMap)

A loop through the five regions applied the Sobel filter to the subimages. Alternatively it would also possible to use the GRASS GIS i.zc module, that also serves the purpose to detect edges. The results of the Sobel filter and the i.zc module are very similar.

    
    # A list to store the b/w maps    
    mapList = []
    # Loop all regions
    for region in regions:
        # Temporary name for the current map
        currentMap = "tmp.%s" % len(mapList)
        gcore.debug("Current map: %s" % currentMap)
        
        # Set the (sub-)region
        gcore.run_command("g.region",
                             w = region['w'],
                             e = region['e'],
                             s = region['s'],
                             n = region['n'])
                             
        # The Sobel edge detection filter
        sobelCmd = "%s = sqrt(pow((-1)*20011130_B.4[-1,-1] + 20011130_B.4[-1,1] - 2*20011130_B.4[0,-1] + 2*20011130_B.4[0,1] - 20011130_B.4[1,-1] + 20011130_B.4[1,1],2) + pow((-1)*20011130_B.4[-1,-1] - 2*20011130_B.4[-1,0] - 20011130_B.4[-1,1] + 20011130_B.4[1,-1] + 2*20011130_B.4[1,0] + 20011130_B.4[1,1],2))" % currentMap
        
        graster.mapcalc(sobelCmd, overwrite = True)

        #gcore.run_command("i.zc", input = "20011130_B.4", output= currentMap )

In the next step I separated the image to 1 and null values using the map calculator.

      
        bwCmd = "{m}.bw = if({m} >= 120, 1, null())".format(m=currentMap)
        graster.mapcalc(bwCmd)

        # Append the current map to patchInput
        mapList.append("%s.bw" % currentMap)

Then I set the region again to the full extent, patched the images together and applied the thinning to the raster file to prepare the vectorization.

        
    # Set the region to the full extent
    gcore.run_command("g.region", rast=join(mapList, ','))
    
    # Merge the raster maps       
    patchedMap = "namngum_river"
    gcore.run_command("r.patch", input = join(mapList, ','), output = patchedMap, overwrite=True)
    
    # Prepare the map to vectorize with r.thin
    thinMap = "namngum_river.thin"
    gcore.run_command("r.thin", input = patchedMap, output = thinMap, overwrite=True)
    
    # Vectorize the map
    vectorMap = "namngum_river_raw"
    gcore.run_command("r.to.vect", input = thinMap, output = vectorMap, type = 'line', overwrite=True)

The final step in the main function was to clean up all temporary raster maps.

    # Clean up all temporary raster maps
    for i in range(len(mapList)):
        maps = "tmp.%s,tmp.%s.bw" % (i, i)
        gcore.run_command("g.remove", rast = maps)

The main entry point for the whole script.

if __name__ == "__main__":
    options, flags = gcore.parser()
    sys.exit(main())

As next step I had to manually clean the resulting vector file. This can be done in different ways, I edited it interactively in GRASS GIS by deleting the false line scraps.

The interactive editing in GRASS GIS

Last but not least I smoothed the lines with v.generalize, the module that advances to one of my favourite in GRASS GIS.

v.generalize input=namngum_river output=namngum_river_smooth method=snakes alpha=0.5 beta=0.5 threshold=0

Nam Ngum river before and after the smoothing

Finally I exported the vector map, converted it to OpenStreetMap's XML format and tagged it correctly in JOSM.

v.out.ogr input=namngum_river_smooth type=line dsn=~/namngum_river.shp
python ogr2osm.py -e 32648 -o ~/namngum_river.osm -a ~/namngum_river.shp

Overlaying the extracted river on Bing imagery in JOSM shows that the result is very accurate.

The resulting river verified with Bing imagery