PipelineResampleToUTM.sml

# PipelineResampleToUTM.sml
# Sample standalone script to demonstrate a simple image pipeline application:
# resampling/reprojecting a source raster to a different coordinate reference
# system.  In this example the image is reprojected to the UTM zone appropriate 
# for its location while maintaining the same datum. 
# 21 February 2008
# Randy Smith, MicroImages, Inc.
# Requires Version 2007:73 or later of the TNT products.
numeric err;										# error code returned
# error checking procedure
proc ReportError(numeric linenum, numeric err) {
	printf("FAILED -line: %d, error: %d\n", linenum, err);
	PopupError(err); 
	}
#####################################  Main program  ##############################
clear();
# set context so script does not exit on error, allowing manual
# handling of errors.
_context.AbortOnError = 0;
# variables for setting up UTM coordinate reference system for the output raster
class SR_COORDREFSYS utmCRS;			# UTM coordinate reference system
class SR_COORDSYS utmCoordSys;		# coordinate system for UTM zone
class SR_COORDOPDEF projectDef;		# projection parameters for the UTM zone
class SR_COORDOPMETHOD method;		# projection method
utmCoordSys.Assign("Projected2D_EN_m");	# assign projected coordinate system
method.Assign("1909");					# assign using numeric ID for Transverse Mercator projection;
												# code can be found on Details panel of Coordinate Reference System window
projectDef.Method = method;			# assign method to projection
# assign projection parameters that are the same for all UTM zones
projectDef.SetParmValueNum("LatitudeOfNaturalOrigin", 0);
projectDef.SetParmValueNum("ScaleFactorAtNaturalOrigin", 0.9996);
projectDef.SetParmValueNum("FalseEasting", 500000);
# CHOOSE INPUT RASTER to be resampled.
class RVC_OBJITEM riObjItem;					# objItem for input raster
DlgGetObject("Select raster to resample:", "Raster", riObjItem, "ExistingOnly");
# PIPELINE SOURCE
# set input raster as source
class IMAGE_PIPELINE_SOURCE_RVC source_In( riObjItem );
err = source_In.Initialize();
if (err < 0) 
	ReportError(_context.CurrentLineNum, err);
else
	print("Pipeline source intialized.");
printf("Source image has %d lines and %d columns.\n",  source_In.GetTotalRows(), source_In.GetTotalColumns() );
# CHECK THAT SOURCE HAS VALID COORDINATE REFERENCE SYSTEM
class IMAGE_PIPELINE_GEOREFERENCE sourceGeoref;
sourceGeoref = source_In.GetGeoreference();
# get coordinate reference system from the source georeference
class SR_COORDREFSYS inCRS;						# input raster's coordinate reference system
inCRS = sourceGeoref.GetCRS();
if (inCRS.IsDefined() == 0 or inCRS.IsLocal() ) {
	PopupMessage("Source coordinate reference system is undefined or local; exiting script.");
	Exit();
	}
else 
	printf("Input image coordinate reference system: %s\n", inCRS.Name );
# FIND GEOGRAPHIC (LAT/LON) COORDINATES OF IMAGE CENTER TO DETERMINE UTM ZONE
class POINT2D center, centerMap;	# point with coordinates of image center
center.x = source_In.GetTotalColumns() * 0.5;	# center in line/column coordinates
center.y = source_In.GetTotalRows() * 0.5;
# get coordinate transformation from image coordinates to map coordinates in its CRS
# in order to find map coordinates of the image center
class TRANS2D_MAPGEN transObjToMap;		# the coordinate transformation
transObjToMap = sourceGeoref.GetTransGen();
# transform image coordinates of image center to map coordinates 
centerMap = transObjToMap.ConvertPoint2DFwd(center);
# if input image's coordinate reference system is projected, convert its center coordinates 
# to Geographic (lat/lon) coordinate system to determine the appropriate UTM zone for the output image
if (inCRS.IsProjected() == 1) 	 
	{
	class SR_COORDSYS geogCoordSys("Ellipsoidal2D_Deg");	# geographic coordinate system
	class SR_COORDREFSYS geogCRS;
	geogCRS.Create(geogCoordSys, inCRS.Datum);
	# set up coordinate transformation from input CRS to geographic coordinates
	class TRANS2D_MAPGEN transToGeog;
	transToGeog.InputCoordRefSys = inCRS;
	transToGeog.OutputCoordRefSys = geogCRS;
	centerMap = transToGeog.ConvertPoint2DFwd(centerMap);	# image center in geographic coordinates
	}
printf("Geographic coordinates of image center: x = %.4f, y = %.4f\n", centerMap.x, centerMap.y);
# SET REMAINING UTM COORDINATE SYSTEM PARAMETERS
# assign false northing based on whether image center is north or south of equator
class STRING ns$;
if ( centerMap.y < 0) {
	projectDef.SetParmValueNum("FalseNorthing", 10000);
	ns$ = "S";
	}
else	{
	projectDef.SetParmValueNum("FalseNorthing", 0);
	ns$ = "N";
	}
# assign longitude of natural origin (central meridian of UTM zone) based on longitude of image center
numeric zoneNum, centMerid;		# UTM zone and longitude of its central meridian
zoneNum = ceil( (180 + centerMap.x) / 6 );
centMerid = (zoneNum * 6) - 183;
projectDef.SetParmValueNum("LongitudeOfNaturalOrigin", centMerid);
printf("Image will be reprojected to UTM zone %d%s\n", zoneNum, ns$);
# CREATE DEFINED OUTPUT UTM COORDINATE REFERENCE SYSTEM WITH SAME DATUM AS INPUT IMAGE CRS
err = utmCRS.Create(utmCoordSys, inCRS.Datum, projectDef);
if (err < 0) 
	ReportError(_context.CurrentLineNum, err);
else
	printf("Output coordinate reference system = %s\n", utmCRS.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 process.
class RVC_OBJITEM rastOutObjItem;
DlgGetObject("Choose raster for resampled output", "Raster", rastOutObjItem, "NewOnly");
# get line and column cell sizes from source's georeference
class POINT2D scaleIn;		# line and column cell sizes as x and y values of POINT2D
sourceGeoref.ComputeScale(center, scaleIn, 1);	# cell sizes in meters
scaleIn.y = abs(scaleIn.y);
printf("Source image cell sizes: line = %.2f m, col = %.2f m\n", scaleIn.y, scaleIn.x);
# PROMPT USER TO ENTER DESIRED OUTPUT LINE AND COLUMN CELL SIZES
numeric lineCellSize, colCellSize;			# specified cell size for the output raster
string prompt$ = "Enter desired line cell size for output raster:";
lineCellSize = PopupNum(prompt$, scaleIn.y, 0, 1000, 2);
prompt$ = "Enter desired column cell size for output raster:";
colCellSize = PopupNum(prompt$, scaleIn.x, 0, 1000, 2);
# SET RESAMPLING METHOD BASED ON RELATIVE CELL SIZE OF INPUT AND OUTPUT
numeric inCellArea, outCellArea;
inCellArea = scaleIn.x * scaleIn.y;
outCellArea = colCellSize * lineCellSize;
string rsmpMethod$;
if (outCellArea > 2 * inCellArea) then rsmpMethod$ = "Average";
	else rsmpMethod$ = "CubicConvolution";
printf("Resampling method: %s\n", rsmpMethod$);
# PIPELINE FILTER to resample source image
class IMAGE_PIPELINE_FILTER_RESAMPLE filter_rsmp(source_In, utmCRS, lineCellSize, colCellSize, rsmpMethod$);
err = filter_rsmp.Initialize();
if (err < 0) 
	ReportError(_context.CurrentLineNum, err);
else
	print("Resample filter initialized.");
# PIPELINE TARGET: set up the target for the pipeline
class IMAGE_PIPELINE_TARGET_RVC target_rvc(filter_rsmp, rastOutObjItem);
target_rvc.SetCompression("DPCM", 0);
err = target_rvc.Initialize();
if (err < 0) 
	ReportError(_context.CurrentLineNum, err);
else
	print("Pipeline target initialized.");
# EXECUTE pipeline process
print("Processing...");
target_rvc.Process();
print("Done.");