Image processing workshop:

Part 1: Optical image processing

Study zone: Jockey's Ridge, East Coast of North Carolina, USA
http://skagit.meas.ncsu.edu/~helena/measwork/jockeys/Page1.html

#######################################################
PART 1: Download, reproject and import
        (not done during workshop due to file size)

1. Download from GLCF (GeoTIFF data are in UTM projection, size around 650MB)
   ftp://ftp.glcf.umiacs.umd.edu/glcf/Landsat/WRS2/p014/r035/p014r035_7x19990923.ETM-EarthSat-Orthorectified/

2. Reproject ("warp") to projection of current location (here NC State Plane)
   and cut out to boundaries of current region. The script header must be edited.
   The resulting channels are written in Erdas/IMG format:
   
   sh landsat7_warp2region.sh

3. Bulk import of all channels (now data projection matches that of location):
   for i in *.img ; do 
      NAME=`echo $i | cut -d'_' -f4`
      r.in.gdal $i out=$NAME
   done

#######################################################
PART 2: Import of data subset in UTM
        (done during workshop)

0. All sample data in original formats are here:
   $HOME/workshop_material/GRASS-imagery/
   on the liveCD.

1. Expand data tarball: landsat7_NC_p014r035_1999_subset.tar.gz
   cd /tmp/
   tar xvfz landsat7_NC_p014r035_1999_subset.tar.gz

2. Start GRASS
   grass62
   
3. We create a location from the data set (startup screen):
   -> Button "Georeferenced file"
      -> Name of new location: nc_utm18
      -> Path to new location: <nothing to change>
      -> path to georeferenced file: /tmp/p014r035_7p19990923_z18_nn80_cut.img
         => Define location
	 
	 The location does not appear automatically, so we go into
	 the "GIS Data Directory" field and press <Enter> there to re-read
	 the list of locations.

    Now we select the new location "nc_utm18" and mapset "PERMANENT" and
    enter GRASS. Note that the region settings are still at 1 pixel.

4. Bulk import of all channels:
   cd /tmp
   for i in *.img ; do 
      NAME=`echo $i | cut -d'_' -f4`
      r.in.gdal $i out=$NAME
   done
  
#######################################################
PART 3: Import and visualization in GRASS

1. RGB display. To start, we first set the region to the blue channel (res. is 28.5m):
   g.region rast=nn10 -p
   d.mon x0
   d.rgb b=nn10 g=nn20 r=nn30

2. Color enhancement:
   i.landsat.rgb b=nn10 g=nn20 r=nn30
   d.rgb b=nn10 g=nn20 r=nn30

   (to revert, just use 'r.colors map col=grey'

3. Color composite:
   r.composite b=nn10 g=nn20 r=nn30 out=landsat.rgb

   This color composite can be viewed in QGIS (enable GRASS Plugin,
   load GRASS raster map) or with:
   d.rast landsat.rgb

4. Panchromatic channel:
   g.region rast=nn80 -p
   d.erase
   d.rast nn80

5. Roads vector map overlay:

   # get projection info:
   ogrinfo -so 91479899.shp 91479899
   
   # see location projection info
   g.proj -we
   
   # reproject to UTM18N (param. order warning 'ogr2ogr .. OUT IN'):
   ogr2ogr -t_srs "`g.proj -wef`" bts_roads_UTM18.shp 91479899.shp
   
   # import: 
   v.in.ogr bts_roads_UTM18.shp out=roads
   d.vect -c roads

4. Aerial infrared photo import and overlay:

   # 71037545.tif is available from USGS (http://nationalmap.gov/)
   #  or NC1 (http://gisdata.usgs.net/website/NC_OneMap/viewer.asp)
   gdalinfo 71037545.tif
   
   #reproject to current location projection (UTM 18N), import:
   # -t_srs is target spatial ref. system; -tr is target resolution
   gdalwarp -t_srs "`g.proj -wef`" -tr 1 1 71037545.tif aerial1998_UTM18.tif
   r.in.gdal aerial1998_UTM18.tif out=airph98

   # set current region/resolution to one of the input maps:
   g.region rast=airph98.blue -p
   r.composite b=airph98.blue g=airph98.green r=airph98.red out=airph98.rgb
   d.erase
   d.rast airph98.rgb
   
   g.region n=n+1000 s=s-1000 w=w-1000 e=e+1000 -p
   d.erase
   d.his i=nn10 h=airph98.rgb

#######################################################
PART 4: Image fusion:

   # Brovey fusion:
   i.fusion.brovey -l ms1=nn10 ms2=nn20 ms3=nn30 pan=nn80 out=lsat.fusion
   g.region -p rast=lsat.fusion.red
   d.erase
   
   # for viz. purposes we swap assignment to blue anc green color gun of d.rgb:
   d.rgb r=lsat.fusion.red b=lsat.fusion.green g=lsat.fusion.blue
   d.vect roads col=yellow

   r.composite r=lsat.fusion.red b=lsat.fusion.green g=lsat.fusion.blue out=lsat.fusion.rgb

   # Perspective view ('ned' is the DEM):
   nviz el=ned col=lsat.fusion.rgb vect=roads

#######################################################
PART 5: Image classification:


   g.region rast=airph98.rgb -p
   d.erase
   d.rast airph98.rgb

   g.list group

   # we add a subgroup from all three RGB channels:
   i.group group=airph98 subgroup=airph98 in=airph98.blue,airph98.green,airph98.red

   # digitizing of training areas (not too big areas...) Use 'x' to save-quit:
   r.digit airph98.train
   d.rast.leg airph98.train

   # generating signatures:
   i.gensigset airph98.train group=airph98 subgroup=airph98 signature=airph98
   
   # assign pixels to classes using combined segmentation/radiometric approach:
   i.smap airph98.train group=airph98 subgroup=airph98 signature=airph98  output=airph98.smap
   d.rast.leg airph98.smap
   
   # filter tiny areas away (5x5 is 25m^2, take dominant class in moving window):
   r.neighbors airph98.smap out=airph98.smap.filt method=mode size=5
   d.rast.leg airph98.smap.filt
   
   # vectorize:
   r.to.vect airph98.smap.filt out=airph98_smap feature=area
   d.rast airph98.rgb
   d.vect airph98_smap
   # ... better use gis.m or QGIS for vector transparency