Armchair mapping with Landsat imagery

In this post I want to show step by step how I extracted hydro power reservoirs and lakes in Laos from up-to-date Landsat imagery for further use in OpenStreetMap. The extraction is done straightforward with common GIS tools but without sophisticated and/or proprietary remote-sensing software.
As an example I take the Nam Leuk hydro power reservoir in Laos. This reservoir needs to be remapped, since the origin user declined the upcoming license change. Currently the reservoir is rather roughly generalized.

Obtain the images
There are different ways how to get Landsat images, my preferred way is the USGS Global Visualization Viewer (short: GLOVIS), a web application that allows comfortable browsing through the Landsat images. Use the GLOVIS preview to make sure that the desired region is not covered by clouds.
The Landsat Product Level 1 contains six visible bands, a thermal infrared and a panchromatic band. For my purpose, to extract water bodies, band 4 seems to be appropriate and sufficient.

Remove the stripes
Due to technical problems with the Landsat sensor there are undesirable stripes in all images. It is possible to remove these stripes using overlapping images taken on two different dates.
First of all I clipped the GeoTIFF images to the required extent with GDAL to accelerate the whole image processing.

gdal_translate -projwin 272500 2050000 287000 2038000 L71128047_04720111126_B40.TIF L71128047_04720111126_B40_NamLeuk.tif

Nam Leuk reservoir with stripes

I merged both images using gdal_merge.py, here it is important to specify the no-data value in the images, in this case 0.

gdal_merge.py -n 0 -o B40_NamLeuk.tif L71128047_04720110211_B40_NamLeuk.tif L71128047_04720111126_B40_NamLeuk.tif

Convert to black-white
Next step is to convert the gray value image to a black-white 2-bit image to separate water from land. A threshold has to be set, all pixel values below the threshold are set to 0 (black), all pixels above to 1 (white). I used ImageMagick, but of course there are different ways how to perform that. The threshold value had to be evaluated by trial and error, the display command is very helpful to preview the image.

convert B40_NamLeuk.tif -threshold 11% png:- | display png:-

After I had found a reasonable threshold I wrote it to a new PNG file

convert B40_NamLeuk.tif -threshold 11% png:namleuk_bw_dirty.png

Still noisy black-white image

Now there were still unwanted, noisy black pixels that needed to be set to white. It was necessary to edit the image manually in GIMP using the pencil tool. The next image shows the clean black and white image.

Clean black-white image

ImageMagick dropped the georeference, that's why it was necessary to reassign it with GDAL and converting the image back to GeoTIFF.

gdal_translate -a_srs EPSG:32648 -a_ullr 272475 2050005 286965 2038005 -of PNG namleuk_bw_clean.png namleuk_bw_clean.tif

Vectorization and smoothing
Having a clean black and white referenced image I used gdal_polygonize to vectorize the image and create a new Shapefile.

gdal_polygonize.py namleuk_bw_clean.tif -f "ESRI Shapefile" namleuk.shp

A new polygon was created for every contiguous region with the same value, i.e. there were polygons with attribute value 0 representing water and value 1 representing land.

Vectorized Nam Leuk shape

Since I was only interested in water polygons I dropped all land polyons

ogr2ogr -where "DN=0" namleuk_clean.shp namleuk.shp

The polygons still showed the original pixels from the Landsat images, that's why I wanted to smooth the geometries.

Before and after smoothing

The following steps I did in GRASS GIS. I created a new GRASS location named namleuk using the georeferenced image and imported the Shapefile.

GRASS 6.4.1 (namleuk):~ > v.in.ogr dsn=~/Desktop/namleuk_clean.shp output=namleuk_clean

To perform the smoothing I used module v.generalize:

GRASS 6.4.1 (namleuk):~ > v.generalize input=namleuk_clean output=namleuk_smooth method=snakes alpha=0.5 beta=0.5

After completing the smoothing I exported the layer again. It was necessary to specify the type as area.

GRASS 6.4.1 (namleuk):~/Data/GISBASE > v.out.ogr -s input=namleuk_smooth type=area dsn=~/Desktop/namleuk_smooth.shp

Uploading to OSM
For the final processing steps it was necessary to convert the Shapefile to an OSM file, clean up the OSM file and upload it.
To convert the layer I used ogr2osm, a very powerful and helpful Python script, that gives full control of the attribute mapping. The script also reproject the input layer if it is not yet in geographic coordinates.

python ogr2osm.py -e 32648 -o ~/Desktop/namleuk.osm ~/Desktop/namleuk_smooth.shp

After opening in JOSM I had to clean up manually the tags and the multipolygon relation, replace the existing reservoir with the newly extracted one and finally, finally upload it to the OpenStreetMap database.

Comparison hand-digitized (blue) and newly extracted reservoir (gray)