Armchair mapping using GRASS GIS

After reading my previous post a friend asked me why I didn't use just GRASS GIS. That's really a valid question. Since the GDAL/OGR utilities are some of my daily tools, I'm more familiar with them than with GRASS. Out of interest I developed a workflow how to extract and vectorize waterbodies from Landsat images using (almost solely) GRASS GIS. The resulting features I wanted to upload to OpenStreetMap.
This time the Nam Ngum 1 hydro power reservoir in Laos was my test area. Like Nam Leuk last time it was only roughly digitized and some of its islands needed to be remapped.

Since GRASS uses own raster and vector data formats, it was necessary to import the images:

GRASS 6.4.1 (Landsat_128047):~ > r.in.gdal input=Laos/Landsat7/LE71280472011042PFS00/L71128047_04720110211_B40.TIF output=20110211_B40
GRASS 6.4.1 (Landsat_128047):~ > r.in.gdal input=Laos/Landsat7/LE71280472011330EDC00/L71128047_04720111126_B40.TIF output=20111126_B40
GRASS 6.4.1 (Landsat_128047):~ > r.in.gdal input=Laos/Landsat7/LE71280472011026PFS00/L71128047_04720110126_B40.TIF output=20110126_B40

Set the computational region to the area of interest:

GRASS 6.4.1 (Landsat_128047):~ > g.region w=231000 n=2082000 e=280000 s=2038000

I merged the three images to remove the stripes. It was necessary to indicate the -z option to set zero values transparent.

GRASS 6.4.1 (Landsat_128047):~ > r.patch -z input=20111126_B40,20110211_B40,20110126_B40 output=NamNgum_B40

Merged Landsat images in GRASS

Then I separated the image values with the map calculator. Alternatively the r.classify module is also an option.

GRASS 6.4.1 (Landsat_128047):~ > r.mapcalc 'NamNgum.bw = if(NamNgum_B40 <= 30, 1, null())'

I took value 30 as threshold, but this depends on the image and needs to be evaluated. I set the threshold rather low to make sure that no land was classified as water. To assign mixed pixel along the lakeshore to water I expanded the waterbodies for one pixel with r.grow.

GRASS 6.4.1 (Landsat_128047):~ > r.grow input=NamNgum.bw output=NamNgum.bw.grow

Like in my previous workflow there were noisy pixels I wanted to get rid of. I had to export the image to edit it in GIMP. Because I set the nodata option to zero, no data pixels are shown transparent in GIMP.

GRASS 6.4.1 (Landsat_128047):~ > r.out.gdal input=NamNgum.bw.grow format=PNG output=~/NamNgum.bw.png nodata=0

In GIMP I cleaned the image with the eraser tool, then I imported it again:

Clean image in GIMP

GRASS 6.4.1 (Landsat_128047):~ > r.in.gdal input=~/NamNgum.bw.png output=NamNgum.clean

The next step was to vectorize the raster data:

GRASS 6.4.1 (Landsat_128047):~ > r.to.vect input=NamNgum.clean output=NamNgum_raw feature=area

As a last processing step I wanted to smooth the newly created vector geometry. The v.generalize module that generalizes and smooths geometries was the only processing step I did in GRASS in my previous workflow.

GRASS 6.4.1 (Landsat_128047):~ > v.generalize input=NamNgum_raw output=NamNgum_smooth method=snakes alpha=0.5 beta=0.5

To convert the data to the osm format I needed a Shapefile. I exported the layer with:

GRASS 6.4.1 (Landsat_128047):~ > v.out.ogr -c input=NamNgum_smooth dsn=~/namngum.shp type=area

The conversion I did like last time with the ogr2osm python script.

python ogr2osm.py -e 32648 -a -o ~/namngum.osm ~/namgnum.shp

Again it was necessary to finally tag the geometries in JOSM and correct the multipolygon relation. A task that needed to be done manually.

Newly created Nam Ngum in the background

Conclusion
The final result is basically the same like in the previous presented workflow, but this one is probably even easier:
With GRASS it is not necessary to clip the images, it's enough to set the computational region. Another needless step is to delete polygons that represents land area because I separated the raster pixels to value 1 or null (i.e. no data) and GRASS ignores the null data pixels during the vectorization.
A disadvantage in GRASS are the additional import and export steps.