Sunday, February 19, 2012

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):~ > input=Laos/Landsat7/LE71280472011042PFS00/L71128047_04720110211_B40.TIF output=20110211_B40
GRASS 6.4.1 (Landsat_128047):~ > input=Laos/Landsat7/LE71280472011330EDC00/L71128047_04720111126_B40.TIF output=20111126_B40
GRASS 6.4.1 (Landsat_128047):~ > 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 ' = 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
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 format=PNG output=~/ 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):~ > input=~/ output=NamNgum.clean
The next step was to vectorize the raster data:
GRASS 6.4.1 (Landsat_128047):~ > 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 -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

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.

No comments:

Post a Comment