Image classification: Difference between revisions
(moved into own page) |
(new tutorial) |
||
Line 31: | Line 31: | ||
* {{cmd|i.smap}} - Performs contextual (image segmentation) image classification using sequential maximum a posteriori (SMAP) estimation. | * {{cmd|i.smap}} - Performs contextual (image segmentation) image classification using sequential maximum a posteriori (SMAP) estimation. | ||
== Small classification tutorial == | |||
Required: GRASS >= 6.4.0 | |||
We are using Landsat satellite data prepared for the [http://www.grassbook.org/data_menu3rd.php North Carolina] sampling site. It is a multispectral dataset where each channel is stored in a separate map. | |||
We are using the North Carolina location here. | |||
=== Warming up with Landsat TM === | |||
# look what maps are there | |||
g.list rast | |||
# set to one of the color channels and print current region settings | |||
g.region rast=lsat7_2002_10 -p | |||
# look at some Landsat map (1: blue, 2: green, 3: red, 4: near infrared etc.) | |||
d.mon x0 | |||
d.rast lsat7_2002_10 | |||
d.rast lsat7_2002_20 | |||
d.rast lsat7_2002_30 | |||
d.rast lsat7_2002_40 | |||
We can easily query the spectrum for individual points: | |||
d.what.rast lsat7_2002_10,lsat7_2002_20,lsat7_2002_30,lsat7_2002_40 | |||
Histogram of a channel: | |||
d.histogram lsat7_2002_10 | |||
Create RGB view on the fly, true and false color: | |||
d.rgb b=lsat7_2002_10 g=lsat7_2002_20 r=lsat7_2002_30 | |||
d.rgb b=lsat7_2002_70 g=lsat7_2002_20 r=lsat7_2002_30 | |||
d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_30 | |||
d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_10 | |||
Check thermal channel (encoded to fit into 0.255 data range): | |||
# http://landsathandbook.gsfc.nasa.gov/handbook/handbook_htmls/chapter11/chapter11.html | |||
# Look up LMIN/LMAX: "Table 11.2 ETM+ Spectral Radiance Range watts/(meter squared * ster * µm)" | |||
# we won't go into details here, r.mapcalc would be used to convert to Kelvoin/degree Celsius. | |||
# look at encoded thermal channel | |||
d.rast lsat7_2002_61 | |||
NDVI as example for map algebra: | |||
r.mapcalc "ndvi = 1.0 * (lsat7_2002_40 - lsat7_2002_30)/(lsat7_2002_40 + lsat7_2002_30)" | |||
d.rast.leg ndvi | |||
r.colors ndvi color=ndvi | |||
d.rast.leg ndvi | |||
Show selected NDVI "classes": | |||
d.rast ndvi val=0.3-1.0 | |||
d.rast ndvi val=-1.0--0.1 | |||
d.rast ndvi val=0.0-0.2 | |||
d.rast ndvi val=0.0-0.3 | |||
Look at panchromatic channel, compare resolution: | |||
g.region rast=lsat7_2002_80 -p | |||
d.erase -f | |||
d.rast lsat7_2002_80 | |||
d.zoom | |||
d.rast lsat7_2002_30 | |||
d.rast lsat7_2002_80 | |||
=== Image Classification === | |||
Set region setting to red image: | |||
g.region rast=lsat7_2002_30 -p | |||
Create a group | |||
i.group group=lsat7_2002 subgroup=lsat7_2002 \ | |||
input=lsat7_2002_10,lsat7_2002_20,lsat7_2002_30,lsat7_2002_40,lsat7_2002_50,lsat7_2002_70 | |||
==== Unsupervised classification==== | |||
Generate unsupervised statistics | |||
i.cluster group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig classes=10 reportfile=lsat7_2002.txt | |||
# look at report file | |||
gedit lsat7_2002.txt | |||
Assign pixels to classes, check quality of assignment: | |||
i.maxlik group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig class=lsat7_2002_class reject=lsat7_2002_reject | |||
d.rast.leg lsat7_2002_class | |||
d.mon x1 | |||
d.rast.leg lsat7_2002_reject | |||
d.mon x2 | |||
d.histogram lsat7_2002_reject | |||
Play with different cluster separations: | |||
i.cluster group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig classes=10 reportfile=lsat7_2002.txt separation=1.5 | |||
i.maxlik group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig class=lsat7_2002_class2 reject=lsat7_2002_reject2 | |||
d.mon x4 | |||
d.rast.leg lsat7_2002_class2 | |||
d.rast.leg lsat7_2002_reject2 | |||
d.rast.leg lsat7_2002_reject | |||
d.rast.leg lsat7_2002_reject2 | |||
d.rast.leg lsat7_2002_reject | |||
d.rast.leg lsat7_2002_reject2 | |||
d.rast.leg lsat7_2002_class2 | |||
Compare to RGB: | |||
d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_10 | |||
d.rast lsat7_2002_class2 cat=1 -o | |||
d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_10 | |||
d.rast lsat7_2002_class cat=1 -o | |||
d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_10 | |||
d.rast lsat7_2002_class cat=1 -o | |||
We do a 3rd attempt with more classes | |||
i.cluster group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig classes=15 reportfile=lsat7_2002.txt | |||
i.maxlik group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig class=lsat7_2002_class3 reject=lsat7_2002_reject3 | |||
d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_10 | |||
d.rast lsat7_2002_class3 cat=1 -o | |||
====Supervised classification==== | |||
Rlaunch GUI for digitizer: | |||
g.gui wxpython | |||
# ... digitize training areas of water, asphalt, forest, uncovered soil; polygons/centroids with attribute | |||
# advantage: we know what the classes are | |||
Assign also numerical ID for each class in new column | |||
v.db.select lsat7_training | |||
v.db.addcol lsat7_training col="id integer" | |||
v.db.select lsat7_training | |||
v.db.update lsat7_training col=id where="name = 'water'" val=1 | |||
v.db.update lsat7_training col=id where="name = 'forest'" val=2 | |||
v.db.update lsat7_training col=id where="name = 'asphalt'" val=3 | |||
v.db.update lsat7_training col=id where="name = 'soil'" val=4 | |||
v.db.select lsat7_training | |||
Transform vector training map to raster model: | |||
v.to.rast in=lsat7_training out=lsat7_training use=attr col=id labelcol=name --o | |||
d.mon x0 | |||
d.rast.leg lsat7_training | |||
Generate statistics from training areas | |||
i.gensigset group=lsat7_2002 subgroup=lsat7_2002 sig=lsat7_2002_smap training=lsat7_training | |||
Perform supervised classification | |||
i.smap group=lsat7_2002 subgroup=lsat7_2002 sig=lsat7_2002_smap out=smap | |||
d.rast.leg smap | |||
d.rgb b=lsat7_2002_10 g=lsat7_2002_20 r=lsat7_2002_30 | |||
Vectorize result: | |||
r.to.vect smap out=smap feat=area | |||
d.rgb b=lsat7_2002_10 g=lsat7_2002_20 r=lsat7_2002_30 | |||
d.vect smap type=boundary col=red | |||
That's it :) | |||
== Further reading classification with GRASS == | == Further reading classification with GRASS == |
Revision as of 16:29, 28 October 2009
Image classification
Classification methods in GRASS
radiometric unsupervised |
radiometric supervised 1 |
radiometric supervised 2 |
radiometric & geometric supervised | |
Preprocessing | i.cluster | i.class (monitor digitizing) | i.gensig (using training maps) | i.gensigset (using training maps) |
Computation | i.maxlik | i.maxlik | i.maxlik | i.smap |
Interactive setup
- i.class - Generates spectral signatures for an image by allowing the user to outline regions of interest.
- The resulting signature file can be used as input for i.maxlik or as a seed signature file for i.cluster.
Processing
- i.cluster - Generates spectral signatures for land cover types in an image using a clustering algorithm.
- The resulting signature file is used as input for i.maxlik, to generate an unsupervised image classification.
- i.gensig - Generates statistics for i.maxlik from raster map layer.
- i.gensigset - Generate statistics for i.smap from raster map layer.
Unsupervised classification
- i.maxlik - Classifies the cell spectral reflectances in imagery data.
- Classification is based on the spectral signature information generated by either i.cluster, i.class, or i.gensig.
Supervised classification
- i.smap - Performs contextual (image segmentation) image classification using sequential maximum a posteriori (SMAP) estimation.
Small classification tutorial
Required: GRASS >= 6.4.0
We are using Landsat satellite data prepared for the North Carolina sampling site. It is a multispectral dataset where each channel is stored in a separate map.
We are using the North Carolina location here.
Warming up with Landsat TM
# look what maps are there g.list rast # set to one of the color channels and print current region settings g.region rast=lsat7_2002_10 -p
# look at some Landsat map (1: blue, 2: green, 3: red, 4: near infrared etc.) d.mon x0 d.rast lsat7_2002_10 d.rast lsat7_2002_20 d.rast lsat7_2002_30 d.rast lsat7_2002_40
We can easily query the spectrum for individual points:
d.what.rast lsat7_2002_10,lsat7_2002_20,lsat7_2002_30,lsat7_2002_40
Histogram of a channel:
d.histogram lsat7_2002_10
Create RGB view on the fly, true and false color:
d.rgb b=lsat7_2002_10 g=lsat7_2002_20 r=lsat7_2002_30 d.rgb b=lsat7_2002_70 g=lsat7_2002_20 r=lsat7_2002_30 d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_30 d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_10
Check thermal channel (encoded to fit into 0.255 data range):
- http://landsathandbook.gsfc.nasa.gov/handbook/handbook_htmls/chapter11/chapter11.html
- Look up LMIN/LMAX: "Table 11.2 ETM+ Spectral Radiance Range watts/(meter squared * ster * µm)"
- we won't go into details here, r.mapcalc would be used to convert to Kelvoin/degree Celsius.
# look at encoded thermal channel d.rast lsat7_2002_61
NDVI as example for map algebra:
r.mapcalc "ndvi = 1.0 * (lsat7_2002_40 - lsat7_2002_30)/(lsat7_2002_40 + lsat7_2002_30)" d.rast.leg ndvi r.colors ndvi color=ndvi d.rast.leg ndvi
Show selected NDVI "classes":
d.rast ndvi val=0.3-1.0 d.rast ndvi val=-1.0--0.1 d.rast ndvi val=0.0-0.2 d.rast ndvi val=0.0-0.3
Look at panchromatic channel, compare resolution:
g.region rast=lsat7_2002_80 -p d.erase -f d.rast lsat7_2002_80 d.zoom d.rast lsat7_2002_30 d.rast lsat7_2002_80
Image Classification
Set region setting to red image:
g.region rast=lsat7_2002_30 -p
Create a group
i.group group=lsat7_2002 subgroup=lsat7_2002 \ input=lsat7_2002_10,lsat7_2002_20,lsat7_2002_30,lsat7_2002_40,lsat7_2002_50,lsat7_2002_70
Unsupervised classification
Generate unsupervised statistics
i.cluster group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig classes=10 reportfile=lsat7_2002.txt
# look at report file gedit lsat7_2002.txt
Assign pixels to classes, check quality of assignment:
i.maxlik group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig class=lsat7_2002_class reject=lsat7_2002_reject d.rast.leg lsat7_2002_class d.mon x1 d.rast.leg lsat7_2002_reject d.mon x2 d.histogram lsat7_2002_reject
Play with different cluster separations:
i.cluster group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig classes=10 reportfile=lsat7_2002.txt separation=1.5 i.maxlik group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig class=lsat7_2002_class2 reject=lsat7_2002_reject2 d.mon x4 d.rast.leg lsat7_2002_class2 d.rast.leg lsat7_2002_reject2 d.rast.leg lsat7_2002_reject d.rast.leg lsat7_2002_reject2 d.rast.leg lsat7_2002_reject d.rast.leg lsat7_2002_reject2 d.rast.leg lsat7_2002_class2
Compare to RGB:
d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_10 d.rast lsat7_2002_class2 cat=1 -o d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_10 d.rast lsat7_2002_class cat=1 -o d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_10 d.rast lsat7_2002_class cat=1 -o
We do a 3rd attempt with more classes
i.cluster group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig classes=15 reportfile=lsat7_2002.txt i.maxlik group=lsat7_2002 subgroup=lsat7_2002 sigfile=lsat7_2002_sig class=lsat7_2002_class3 reject=lsat7_2002_reject3 d.rgb b=lsat7_2002_70 g=lsat7_2002_50 r=lsat7_2002_10 d.rast lsat7_2002_class3 cat=1 -o
Supervised classification
Rlaunch GUI for digitizer:
g.gui wxpython # ... digitize training areas of water, asphalt, forest, uncovered soil; polygons/centroids with attribute # advantage: we know what the classes are
Assign also numerical ID for each class in new column
v.db.select lsat7_training v.db.addcol lsat7_training col="id integer" v.db.select lsat7_training v.db.update lsat7_training col=id where="name = 'water'" val=1 v.db.update lsat7_training col=id where="name = 'forest'" val=2 v.db.update lsat7_training col=id where="name = 'asphalt'" val=3 v.db.update lsat7_training col=id where="name = 'soil'" val=4 v.db.select lsat7_training
Transform vector training map to raster model:
v.to.rast in=lsat7_training out=lsat7_training use=attr col=id labelcol=name --o d.mon x0 d.rast.leg lsat7_training
Generate statistics from training areas
i.gensigset group=lsat7_2002 subgroup=lsat7_2002 sig=lsat7_2002_smap training=lsat7_training
Perform supervised classification
i.smap group=lsat7_2002 subgroup=lsat7_2002 sig=lsat7_2002_smap out=smap
d.rast.leg smap d.rgb b=lsat7_2002_10 g=lsat7_2002_20 r=lsat7_2002_30
Vectorize result:
r.to.vect smap out=smap feat=area d.rgb b=lsat7_2002_10 g=lsat7_2002_20 r=lsat7_2002_30 d.vect smap type=boundary col=red
That's it :)
Further reading classification with GRASS
- Micha Silver: Analyzing acacia tree health in the Arava with GRASS GIS
- Perrygeo: Impervious surface deliniation with GRASS
- Dylan Beaudette: Working with Landsat Data
- Dylan Beaudette: Canopy Quantification via Image Classification