More scripts: Raster
Syntax Highlighing:
comments, key words, predefined symbols, class members & methods, functions & classes
# PipelineNDVIfromTIFF_CropAndMaskFromRegion.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 computations are restricted to areas of the image delineated by
# polygons in an input vector object. The polygons are converted to a region which is used to
# mask the areas outside the polygons and to crop the image to the extent of the polygons.
# 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.
# 9 January 2008
# Randy Smith, MicroImages, Inc.
# updated 17 January 2014
# 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"));
# 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");
}
printf("Number of bands in TIFF file = %d\n", source_tif.GetNumSamples() );
if (source_tif.GetNumSamples() <> 4) {
print("Selected TIFF file does not contain 4 bands; exiting now.");
Exit();
}
# Get georeference and coordinate reference system from the GeoTIFF image.
class IMAGE_PIPELINE_GEOREFERENCE georefImg;
georefImg = source_tif.GetGeoreference();
class SR_COORDREFSYS imgCRS = georefImg.GetCRS();
if (!imgCRS.IsDefined() ) {
print("Image is not georeferenced; exiting now.");
Exit();
}
else if (imgCRS.IsLocal() ) {
print("Image has local engineering georeference; exiting now.");
Exit();
}
else
printf("Image coordinate reference system = %s\n", imgCRS.Name );
# Get coordinate transformation from image to its defined CRS
class TRANS2D_MAPGEN imgToCRS = georefImg.GetTransGen();
# CHOOSE VECTOR OBJECT WITH POLYGONS TO INDICATE AREAS TO PROCESS
class RVC_VECTOR polyVect;
class RVC_OBJITEM vectObjItem;
class RVC_GEOREFERENCE vectGeoref;
class SR_COORDREFSYS vectCRS;
DlgGetObject("Choose vector object with polygons outlining areas to process:", "Vector", vectObjItem, "ExistingOnly");
polyVect.Open(vectObjItem, "Read");
if (polyVect.$Info.NumPolys == 0) {
print("Vector object contains no polygons; exiting now.");
Exit();
}
err = polyVect.GetDefaultGeoref(vectGeoref);
if (err < 0 ) {
ReportError(_context.CurrentLineNum, err);
print("Polygon vector must be georeferenced and is not; exiting now.");
Exit();
}
vectCRS = vectGeoref.GetCoordRefSys();
printf("Vector coordinate reference system = %s\n", vectCRS.Name);
# 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 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$ );
# CONVERT INPUT VECTOR TO A REGION FOR USE AS A PIPELINE SOURCE
class REGION2D RegFromVect = ConvertVectToRegion(polyVect, GetLastUsedGeorefObject(polyVect) );
# if vector has different coordinate reference system from image, reproject region
if (vectCRS.Name <> imgCRS.Name) {
printf("Converting region to image CRS: %s\n", imgCRS.Name);
RegFromVect.ConvertTo(imgCRS);
}
class RECT vecRegionExtents = RegFromVect.Extents;
printf("Minimum map extents of polyVect region: %.2f, %.2f\n", vecRegionExtents.x1, vecRegionExtents.y1);
printf("Maximum map extents of polyVect region: %.2f, %.2f\n\n", vecRegionExtents.x2, vecRegionExtents.y2);
# GET EXTENTS OF THE IMAGE SOURCE AS A REGION
class REGION2D imgReg;
err = filter_select.ComputeGeoreferenceRegion(imgReg);
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
# CHECK THAT REGION FROM VECTOR IS CONTAINED WITHIN IMAGE REGION
if (imgReg.TestRegion(RegFromVect, "FullInside") ) {
print("Image region contains vector extents.");
}
else {
print("Input vector is not contained within the source image. Exiting now.");
Exit();
}
# PIPELINE SOURCE FOR REGION CREATED FROM THE VECTOR POLYGONS
# Construct using source image as a reference so the "image dimensions" of the region source match,
# enabling the region source to be used as a mask in the MASKVALIDITY filter
class IMAGE_PIPELINE_SOURCE_REGION sourceReg(RegFromVect, filter_select);
err = sourceReg.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized region source.");
printf("Size of region = %d lines, %d columns\n", sourceReg.GetTotalRows(), sourceReg.GetTotalColumns() );
}
# Get the extents of the region as a RECT in image coordinates
# to use to crop the source image. These are the extents of the region used to
# create the region source, not those of the image used as its reference.
class RECT regExtents;
regExtents = sourceReg.GetRegionExtents();
printf("Minimum values of cropping rectangle from polygons: x = %.2f, y = %.2f\n", regExtents.x1, regExtents.y1);
printf("Maximum values of cropping rectangle from polygons: x = %.2f, y = %.2f\n", regExtents.x2, regExtents.y2);
# SET UP PIPELINE FILTER TO MASK PORTIONS OF SOURCE IMAGE
# mask area outside the region created from the vector polygons
class IMAGE_PIPELINE_FILTER_MASKVALIDITY filterMask(filter_select, sourceReg);
err = filterMask.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized image mask filter.");
}
# SET UP PIPELINE FILTER TO CROP THE SOURCE IMAGE
class IMAGE_PIPELINE_FILTER_CROP filter_Crop(filterMask, regExtents);
err = filter_Crop.Initialize();
if (err < 0 )
ReportError(_context.CurrentLineNum, err);
else {
print("Initialized image crop filter.");
}
# 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_Crop, "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