coastal.sml

# SML created 7 May 2003 by Dan Glasser
# Script takes an input elevation raster (suggest cb_elev.rvc/DEM_16bit).
#
# Based on the range values set in the script (specific to DEM_16bit),
# an output classification raster is created.
#
# This script uses 3 elevation ranges and creates a color map to separate
# the 3 classes 0-1099 = red, 1100-1249 = green, >=1250 = blue in output raster.
#
# Georeferencing subobjects are copied from the input raster.
#
# To create the appropriate legend a table is generated that stores the raster
# value and it's associated legend label.  Setting up the legend requires a 
# "Color Sample" legend with the [Label...] set to Label.Label
#
# To modify change the number of classes as desired, add entries to the array
# "range[]".  Make appropriate change to foreach loop setting class values.
#

$warnings 3  # Set warning level 3

numeric classes = 3;

class RVC_RASTER Rin, Rclassified;
class RVC_COLORMAP colormap;
numeric lin, col;
class DATABASE db;
class DBTABLEINFO tableinfo, internal;
class DBFIELDINFO valuefield, labelfield, internalvalue;

GetInputRaster(Rin);

GetOutputRaster(Rclassified, NumLins(Rin), NumCols(Rin), "8-bit unsigned");

# make separate classes here
array numeric range[classes];
range[1] = 0;
range[2] = 1100;
range[3] = 1250;

# class values are 1, 2, and 3.  0 is Null
foreach Rin[lin, col]
{
	if(Rin >= range[3])
		Rclassified[lin,col] = 3;
	else if (Rin >= range[2])
		Rclassified[lin,col] = 2;
	else if (Rin >= range[1])
		Rclassified[lin,col] = 1;
}
SetNull(Rclassified, 0);

# set color map (1=R, 2=G, 3=B)
ColorMapSetColorRGB(colormap, 1, 255, 0, 0);
ColorMapSetColorRGB(colormap, 2, 0, 255, 0);
ColorMapSetColorRGB(colormap, 3, 0, 0, 255);
ColorMapWriteToRastVar(Rclassified, colormap, "ClassifyColormap", "Color map created for classification");

# open database so that legend labels can be set up
db = OpenRasterDatabase(Rclassified);
internal = DatabaseGetTableInfo(db, "Internal");
internalvalue = FieldGetInfoByName(internal, "Value");

# create table and add fields for value and label
tableinfo = TableCreate(db, "Label", "Legend labels");
valuefield = TableAddFieldInteger(tableinfo, "Value", 4);
labelfield = TableAddFieldString(tableinfo, "Label", 20);

# setup up relation to relate table values/labels to internal values
tableinfo.RelatedOnly = 1;
internalvalue.IsIndexed = 1;
valuefield.IsIndexed = 1;
valuefield.IsKey = 1;
valuefield.Key = internalvalue;

# create records for legend labels
numeric i;
for i=1 to classes-1
{
	TableNewRecord(tableinfo, i, sprintf("%d-%d", range[i], range[i+1]-1));
}
TableNewRecord(tableinfo, i, sprintf(">=%d", range[classes]));

print ("Cols:  ", Rclassified.$Info.NumCols, " Lins:  ",  Rclassified.$Info.NumLins, " Type:  ", Rclassified.$Info.Type);

CloseRaster(Rin);
CloseRaster(Rclassified);

CopySubobjects(Rin, Rclassified, "georef");