More scripts: Raster
Syntax Highlighing:
comments, key words, predefined symbols, class members & methods, functions & classes
# PipelineNDVIfromTIFF.sml
# Sample standalone script to demonstrate an image pipeline application.
# Computes a Normalized Difference Vegetation Index (NDVI) raster from Red and
# Near-infrared bands directly from an input GeoTIFF file containing a 4-band multispectral
# QuickBird or Ikonos satellite image. The GeoTIFF file does not need to have been previously linked
# or imported to TNTmips.
# NDVI is calculated for each image cell using the formula:
# 1000 * (Near Infrared - Red) / (Near Infrared + Red)
# Values in the resulting NDVI raster can range from -1000 to + 1000 (using Signed 16-bit raster type).
# If the input image is projected, there may be 0-value cells around the edges. A ValidityNear filter
# is used to identify cells that have zero value in both bands; the filter creates a null mask
# for each band with such cells marked as null.
# Path radiance (atmospheric scattering / haze) is an additive effect that varies between
# bands, and so should be removed before calculating NDVI. This script uses the minimum (non-null)
# raster value in each band as an approximation to the Path Radiance value for that band. This value
# is subtracted using the SCALEOFFSET filter (scale = 1.00, offset value negative).
# The LINEAR filter (linear combination of band values) is used to perform the addition and
# subtraction of the band values.
# The NDVI image is produced as a single 16-bit signed integer RVC raster in a Project File.
# 8 January 2008
# Randy Smith, MicroImages, Inc.
# Requires Version 2008:74 or later of the TNT products.
numeric err;
# error checking procedure
proc ReportError(numeric linenum, numeric err) {
printf("FAILED -line: %d, error: %d\n", linenum - 1, err);
PopupError(err);
}
################################################################################
################################# Main program ###############################
clear();
# set context so script does not exit on error, allowing manual
# handling of errors.
_context.AbortOnError = 0;
# CHOOSE GEOTIFF FILE containing QuickBird or IKONOS 4-band multispectral image
class FILEPATH inFilePath(GetInputFileName("","Choose GeoTIFF file with source image:", ".tif"));
# CHOOSE OUTPUT RASTER
# DlgGetObject creates only an ObjItem for the output; the new raster will then be created
# by the TARGET_RVC with the desired compression. The ObjectPath of the ObjItem is not
# valid until the output raster is actually created by the pipeline.
class RVC_OBJITEM ndviOutObjItem;
DlgGetObject("Choose new raster for NDVI result:", "Raster", ndviOutObjItem, "NewOnly");
# PIPELINE SOURCE
# set GeoTIFF file containing a four-band QuickBird or Ikonos image as source
class IMAGE_PIPELINE_SOURCE_TIFF source_tif( inFilePath );
err = source_tif.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Source TIFF file initialized");
}
if (source_tif.GetNumSamples() == 4)
printf("Number of bands in TIFF file = %d\n", source_tif.GetNumSamples() );
else {
print("Selected TIFF file does not contain 4 bands; exiting now.");
Exit();
}
# Get georeference from the GeoTIFF file to print what CRS is used.
class IMAGE_PIPELINE_GEOREFERENCE georef;
georef = source_tif.GetGeoreference();
class SR_COORDREFSYS crs = georef.GetCRS();
printf("Coordinate reference system = %s\n", crs.Name );
# PIPELINE FILTER_SELECT to select red and near-infrared bands for processing (don't need to carry the rest
# forward through the pipeline).
# The source TIFF file consists of 4 bands (samples) per pixel, which have numeric indices 0 = blue,
# 1 = green, 2 = red, 3 = near-infrared. These samples are selected by index number in this filter.
numeric numSamples = 2; # number of bands from 4-band TIFF file to be used
array numeric srcSamples[2]; # array to contain the numeric indices of the bands to select
srcSamples[1] = 2; # index for red band
srcSamples[2] = 3; # index for near-infrared band
class IMAGE_PIPELINE_FILTER_SELECT filter_select(source_tif, srcSamples, numSamples);
err = filter_select.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized band selection filter");
}
# the result of this FILTER_SELECT is an "image" with two samples per pixel, with the red value as index 0
# and the near-infrared value as index 1.
class STRING pixelType = filter_select.GetPixelType();
printf("Pixel type in image = %s\n", pixelType);
class STRING datatypeRed$, datatypeNIR$;
class IMAGE_PIPELINE_PIXEL pixel;
pixel = filter_select.GetPixelValueMax();
datatypeRed$ = pixel.GetSample(0).GetDataType();
datatypeNIR$ = pixel.GetSample(1).GetDataType();
printf("Datatype of red band = %s\n", datatypeRed$ );
printf("Datatype of near-infrared band = %s\n", datatypeNIR$ );
# PIPELINE FILTER_VALIDITYNEAR
# Filter to set null for cells based on nearness (distance) to a specified value
# In this case test for zero value in both bands (assumed to be non-image cells) and set
# such cells to null (0 value in null mask created by the filter)
# validity flag
numeric resultValid = 0; # set 1 if pixels near test value should be valid, 0 if invalid;
# we want 0 values to be invalid, so set this parameter to false (0)
numeric distance = 0.5; # distance to test must be greater than 0; use a decimal value less
# than 1 to test for exact match to an integer value
class IMAGE_PIPELINE_FILTER_VALIDITYNEAR validity(filter_select, "Orthogonal", resultValid);
# define two 0-value samples to provide the set of values
class IMAGE_PIPELINE_SAMPLE redNullSamp("Gray", "UINT16", 0);
class IMAGE_PIPELINE_SAMPLE nirNullSamp("Gray", "UINT16", 0);
# define a text pixel with the two 0-value samples to pass to the validity filter
class IMAGE_PIPELINE_PIXEL testPixel(redNullSamp, nirNullSamp);
# add the pixel value to test and the test distance to the validity filter;
validity.AddValue(testPixel, distance);
array numeric samplesToTest[2]; # array to set which image samples to test (in this case, both)
samplesToTest[1] = 0; # sample index for red
samplesToTest[2] = 1; # sample index for near-infrared
validity.SetSamplesTest(samplesToTest, 2);
err = validity.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized ValidityNear filter.");
}
# GET MINIMUM AND MAXIMUM VALUES of the red and near-infrared bands.
# First call method on stage to compute the range of values in the bands; using the "Exact"
# value for the RangeType parameter in this method to compute the actual range of values using the
# rasters without sampling (i.e. computing using the pyramids, if any). The "Default" value returns
# the range of the raster datatype.
validity.ComputePixelRange("Exact");
class IMAGE_PIPELINE_PIXEL minVals, maxVals; # container to hold the minimum and maximum value for each band
minVals = validity.GetPixelValueMin();
maxVals = validity.GetPixelValueMax();
printf("Minimum value of red band = %d\n", minVals.GetSample(0).GetValue() );
printf("Maximum value of red band = %d\n", maxVals.GetSample(0).GetValue() );
printf("Minimum value of near-infrared band = %d\n", minVals.GetSample(1).GetValue() );
printf("Maximum value of near-infrared band = %d\n", maxVals.GetSample(1).GetValue() );
# PIPELINE FILTER_SCALEOFFSET applies scale (multiply by constant) and offset (add constant value) to each cell
# Use with scale = 1 and negative offset to subtract the minimum valid value from each of the bands
# as an approximate correction for path-radiance.
class IMAGE_PIPELINE_FILTER_SCALEOFFSET offsetMin(validity);
# set scale = 1 and negative offset equal to minimum value for red band (0) and near-infrared band (1)
offsetMin.SetScaleOffset(0, 1, minVals.GetSample(0).GetValue() * -1);
offsetMin.SetScaleOffset(1, 1, minVals.GetSample(1).GetValue() * -1);
err = offsetMin.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized Scale-offset filter for path-radiance correction.");
}
# PIPELINE FILTER_LINEAR produces a linear combination of the band values using a coefficient for each band
# Cell Value = a * Band1value + b * Band2value + ...
# Use here to subtract band values to provide dividend for the NDVI expression: NIR - Red;
# Set number of samples in target to 1; coefficient for NIR = 1 and for Red = -1
numeric numTargetSamples = 1;
numeric targetSampleNum = 0; # 0th target sample
numeric sourceSampleNum = 0; # 0th source sample (red)
numeric coefValue; # value of the coefficient
class IMAGE_PIPELINE_FILTER_LINEAR linearSubtract(offsetMin, numTargetSamples);
coefValue = -1; # set coefficient value = -1 for red
linearSubtract.SetCoefficient(targetSampleNum, sourceSampleNum, coefValue);
sourceSampleNum = 1; # 1st (near-infrared) source sample
coefValue = 1; # set coefficient value = 1 for near-infrared
linearSubtract.SetCoefficient(targetSampleNum, sourceSampleNum, coefValue);
err = linearSubtract.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized linear filter for subtracting red from near-infrared.");
}
# PIPELINE FILTER_LINEAR
# Use here to add band values to provide divisor for the NDVI expression: NIR + Red;
# Set number of samples in target to 1; coefficient for both NIR and Red = 1;
class IMAGE_PIPELINE_FILTER_LINEAR linearAdd(offsetMin, numTargetSamples);
sourceSampleNum = 0; # 0th source sample (red)
linearAdd.SetCoefficient(targetSampleNum, sourceSampleNum, coefValue);
sourceSampleNum = 1; # 1st source sample (near-infrared)
linearAdd.SetCoefficient(targetSampleNum, sourceSampleNum, coefValue);
err = linearAdd.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized linear filter for adding red and near-infrared.");
}
# PIPELINE FILTER_DIVIDE: filter to divide two images
# Use to perform the final NDVI division using the results from the two Linear filters
# The stage takes two source stages, the dividend stage followed by the divisor stage,
# Output of this filter is a single 64-bit floating-point image with possible range -1.00 to + 1.00
class IMAGE_PIPELINE_FILTER_DIVIDE divide(linearSubtract, linearAdd);
err = divide.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized divide filter to compute NDVI.");
}
# PIPELINE FILTER_SCALEOFFSET
# Scale floating-point values to range -1000 to + 1000 using scale factor = 1000 and offset = 0
class IMAGE_PIPELINE_FILTER_SCALEOFFSET scale1000(divide);
scale1000.SetScaleOffset(0, 1000, 0);
err = scale1000.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized scale-offset filter to rescale NDVI.");
}
# PIPELINE FILTER_DATATYPE: filter to change datatype
# Change rescaled NDVI from 64-bit floating point to 16-bit signed integer
class IMAGE_PIPELINE_FILTER_DATATYPE changeToSINT16(scale1000, "SINT16", "Invalidate");
err = changeToSINT16.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized datatype filter to convert to signed 16-bit.");
}
# PIPELINE TARGET
class IMAGE_PIPELINE_TARGET_RVC target_rvc(changeToSINT16, ndviOutObjItem);
target_rvc.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized NDVI target.");
}
# EXECUTE the pipeline process
print("Processing...");
target_rvc.Process();
print("Done.");
Download