ProvisionTNTscript.sml

  Download

More scripts: TNTscript

Syntax Highlighing:

comments, key words, predefined symbols, class members & methods, functions & classes
            
# ProvisionTNTscript.sml
# Image provisioning script to be run in TNTscript using extraction area (rectangular extents or polygon) provided by web client.
# Version 21 July 2015
#		- output formats: RVC (composite or separates), GeoTIFF, GeoJP2, PNG, JPEG, KMZ
#		- output image sized to the selected limit: 1, 2, or 4 megapixels
# uncomment the following preprocessor command to run in test mode
# with hard-coded  parameters
#$define TEST
##############   Variables to be provided by the web interface   #################
string taskID$;
string userEmail$; 
string rvcPath$;		# string with filepath to source RVC file
string rvcObject$;	# string with name of RVC raster (not the object path)
string outFormat$;	# output file format
										# possible values:  "RVC", "RVC_RGB", "GeoTIFF", "GeoJP2", "JPEG",  "KMZ, or "PNG"
string CRS$;				# EPSG code for the CRS selected by user for the output
numeric lossyComp;		# compression ratio for JP2 lossy compression, 0 for lossless
string limit$;			# size limit: possible values "1", "2" or "4"
# extraction area parameters with coordinates in Spherical Web Mercator:
string area$;			# coordinate pairs for extraction polygon vertices (empty string for rectangle)
numeric minx, miny, maxx, maxy;		# min and max extents for rectangular extraction area
##############################################################
################  Parameters for test mode  ###################
$ifdef TEST
		taskID$ = "ProvisionTNTscriptArea";
		#rvcPath$ =  "G:/2008/OK2008 NC UTM14 1m 15to1/OK2008 NC UTM14.rvc";
		#rvcObject$ = "OK 2008 NC UTM14 1n 15to1";
		rvcPath$ = "E:/SMLtest/provisioning/New_SML_Autotests/provision/ProvisionNE.rvc";
		rvcObject$ = "NE 2006 NC UTM14 1m 15to1";
		outFormat$ = "RVC";
		userEmail$ = "[email protected]";
		CRS$ = "EPSG:26915";
		limit$ = "2";
		minx = -11453574;	miny = 5269632;  maxx = -11445486;  maxy = 5275985;
#		area$ = "";
		area$ = "-11452753 5270359,-11451980 5274860,-11447455 5275282,-11446611 5271766,-11448932 5270359,-11452753 5270359";
		lossyComp = 15;
$endif
##############################################################
#############   Other variables   #############
class RVC_RASTER RastIn;		# the source image raster
class REGION2D exRegion;		# region to use for extraction
class POLYLINE poly;				# polyline outlining extraction area
class POINT2D vertex;			# point for next vertex to be added to polyline
class STRING pair$;					# next substring from area$ containing x & y coordinates
class DATETIME currentDT;		# the current local date/time
class EMAIL_MESSAGE email;
class STRING logfileName$;
class STRING readmeFilename$;
class FILE logfile;
class FILE readme;
class STRING outpath$;			# path string for output www jobs directory
class STRING outDirPath$;		# path string for output directory created in www jobs directory
class FILEPATH outFilepath;	# filepath for output directory
class STRING outputFileExtension$;
numeric err;		# error value from defined functions
string error$;		# string with error text
numeric i;				# loop counter
class ZIPFILE zip;
class FILEPATH zipFilepath;
string start$,  finish$, errpath$;
######  Variables for resampling to new cell size to meet image size limit #############
numeric numLins, numCols;		# number of lines and columns in extraction area
numeric totalCells;						# number of cells in extracted image
numeric limitCells;							# maximum allowed number of cells in output
class POINT2D rastScale;			# source raster x and y scales computed from georeference
numeric rescale;					# scale factor for rescaling line and column scales for output
string inCRS$;						#  EPSG code for CRS of source image
class SR_COORDREFSYS inCRS;		# CRS of source image
class SR_COORDREFSYS outCRS;	# output CRS
# error notification procedure
proc notifyError(string err$) {
	string php$ = "errpath$ " +  taskID$ + " > error_out"+taskID$;
	fprint(logfile, err$);
	fclose(logfile);
	fprint(readme, err$);
	fclose(readme); 
	email.SetSubject(sprintf("Error with provisioning job %s", taskID$) );
	email.SetBody(err$);
	email.SetRecipient("[email protected]");
	email.SendMail();
	email.SetRecipient("[email protected]");
	email.SendMail();
	ExecuteProcess(php$);
	Exit();
	}
# error checking procedure
func string ReportError(numeric linenum, numeric err) {
	string err$;
	err$ = sprintf("FAILED -line: %d, error: %d\n", linenum - 1, err);
	return err$; 
	}
#############################################################
#########################  Program  ###########################
clear();
### set up to run php scripts to write info to MySQL database
start$ = "C:/wamp/bin/php/php5.5.12/php.exe c:/wamp/www/provisioning/start.php " +  taskID$;
finish$ = "C:/wamp/bin/php/php5.5.12/php.exe c:/wamp/www/provisioning/finish.php " +  taskID$;
errpath$ = "C:/wamp/bin/php/php5.5.12/php.exe c:/wamp/www/provisioning/error.php";
ExecuteProcess(start$);
 
########################  Set up Log File  ######################
# open log file to append status information
logfileName$ = _context.ScriptDir + "/Provisioning.log";
logfile = fopen(logfileName$, "a");
print("Logfile opened. logfileName$");
currentDT.SetCurrent();
fprintf(logfile, "\nJob initiated %s Central Standard Time\n", currentDT);
 
$ifdef TEST
		# print job info to console
		printf("TaskID: %s\n", taskID$);
		printf("User Email: %s\n", userEmail$);
		printf("RVCpath: %s\n", rvcPath$);
		printf("RVCobject: %s\n", rvcObject$);
		printf("Selected output format: %s\n", outFormat$);
		printf("Selected CRS: %s\n", CRS$);
		printf("Limit = %s\n", limit$);
		printf("LossyComp value = %d\n", lossyComp);
$endif
# print job info to log file
fprintf(logfile, "TaskID: %s\n", taskID$);
fprintf(logfile, "User Email: %s\n", userEmail$);
fprintf(logfile, "RVCpath: %s\n", rvcPath$);
fprintf(logfile, "RVCobject: %s\n", rvcObject$);
fprintf(logfile, "Selected output format: %s\n", outFormat$);
fprintf(logfile, "Selected CRS: %s\n", CRS$);
fprintf(logfile, "Limit = %s\n", limit$);
fprintf(logfile, "LossyComp value = %d\n", lossyComp);
#####################  Other  File Handling  ######################
# set up output directory
$ ifdef TEST
			outpath$ = _context.ScriptDir + "/";
$ else
			outpath$ = "C:/wamp/www/jobs/"; 
$ endif
outDirPath$ = outpath$ + taskID$;		# new directory to be created for this job
CreateDir(outDirPath$);						# create new directory
outFilepath = outDirPath$ + "/" + taskID$;	# filepath for new directory
# set filepath for zipfile in output directory
if (outFormat$ == "KMZ") then
	zipFilepath = outDirPath$ + "/" + taskID$ + ".kmz";
else
	zipFilepath = outDirPath$ + "/" + taskID$ + ".zip";
$ifdef TEST											# delete previously existing zip file with same path/name
	if ( fexists(zipFilepath) ) then
		DeleteFile(zipFilepath);
$endif
readmeFilename$ = outDirPath$ + "/readme" + taskID$ + ".txt";
# open readme file
readme = fopen(readmeFilename$, "a");
fprintf(logfile, "New directory path: %s\n", outDirPath$);
fprintf(logfile, "outFilepath = %s\n", outFilepath.GetPath() );
fprintf(logfile, "Zip file path = %s\n", zipFilepath.GetPath() );
fprintf(logfile, "readmeFilename$ = %s\n", readmeFilename$);
$ifdef TEST
		printf("New directory path: %s\n", outDirPath$);
		printf("Zip file path = %s\n", zipFilepath.GetPath() + "/" +  zipFilepath.GetName() );
		printf("readmeFilename$ = %s\n", readmeFilename$);
$endif
# set up e-mail source
email.SetSourceAddress("[email protected]");
email.SetSourceName("GeoProvisioning at MicroImages.com");
#######################  Source Image  #################################
# open the input raster to get an RVC_OBJITEM to use to set up a pipeline source
OpenRaster(RastIn, rvcPath$, rvcObject$);
class RVC_OBJITEM sourceObjItem;
sourceObjItem = RastIn.GetObjItem();
if (!sourceObjItem.IsExisting() ) {
	error$ = "Source raster cannot be found.";
	notifyError(error$);
	}
##########################################################
##########  set up pipeline source for source image  #############
class IMAGE_PIPELINE_SOURCE_RVC sourceRVC(sourceObjItem);
err = sourceRVC.Initialize();
if (err < 0)  {
	error$ = ReportError(_context.CurrentLineNum, err);
	notifyError(error$);
	}
# get pipeline georeference from the source image
class IMAGE_PIPELINE_GEOREFERENCE georefSource;		# georeference from source
georefSource = sourceRVC.GetGeoreference();
inCRS = georefSource.GetCRS();
inCRS$ = "EPSG:" + inCRS.GetID("EPSG");
# get coordinate transformation from source image to its defined CRS
class TRANS2D_MAPGEN transImgToCRS = georefSource.GetTransGen();
#############################################################
##################  get extraction area  #########################
if (area$ <> "")		# area$ has vertex coordinates from polygon tool
	{
	# parse area$ to set up polyline for creating region for area of interest
	fprint(logfile, "Vertex coordinates read from log file:");
	for i = 1 to NumberTokens(area$, ",")
		{
		pair$ = GetToken(area$, ",", i);
		vertex.x = StrToNum(GetToken(pair$, " ", 1) );
		vertex.y = StrToNum(GetToken(pair$, " ", 2) );
		poly.AppendVertex(vertex);
		fprintf(logfile, "vertex %d: x = %.6f, y = %.6f\n", i, vertex.x, vertex.y);
		}
	# create region from polyline by union of empty region with polyline
	exRegion.UnionPolyLine(poly);
	}
else 		# area$ is empty, rectangle tool was used, get rectangle extents values
	{
	class RECT exRect;
	exRect.x1 = minx;		exRect.y1 = miny;		exRect.x2 = maxx;			exRect.y2 = maxy;
	# create region from the rectangle extents
	exRegion.SetFromRect(exRect);
	}
##########################
# set CRS for extraction region and convert to CRS of source image
class SR_COORDREFSYS toolCRS("EPSG:3857");		# coordinate reference system of extraction area in web interface
exRegion.CoordRefSys = toolCRS;
exRegion.ConvertTo(inCRS);
#################################
class POINT2D centerPt, centerPtImg;
# find center point of extraction area in map coordinates
centerPt.x = (exRegion.Extents.x2 - exRegion.Extents.x1) * 0.5;
centerPt.y = (exRegion.Extents.y2 - exRegion.Extents.y1) * 0.5;
# compute center point of extraction area in image coordinates
centerPtImg = transImgToCRS.ConvertPoint2DInv(centerPt);
# compute cell sizes (scales) at center of extraction area of source image from its georeference
georefSource.ComputeScale(centerPtImg, rastScale, 0);
fprint(logfile, "\nSource raster info:");
fprintf(logfile, "Raster cols: %d, lines: %d\n", sourceRVC.GetTotalColumns(), sourceRVC.GetTotalRows() );
fprintf(logfile, "Source CRS: %s, %s\n", inCRS$, inCRS.Name);
fprintf(logfile, "Native cell size of source image: x = %.3f, y = %.3f\n", rastScale.x, rastScale.y);
fprint(logfile, "Extents of area to extract in image CRS:");
fprintf(logfile, "minx: %.6f, miny: %.6f, maxx: %.6f, max y: %.6f\n", exRegion.Extents.x1, exRegion.Extents.y1, exRegion.Extents.x2, exRegion.Extents.y2);
fprintf(readme, "Details of your provisioning order %s:", taskID$);
fprintf(readme, "\nCoordinate Reference System of the source image: %s\n", inCRS.Name);
fprintf(readme, "The native cell size of the source image in meters is: x = %.3f, y = %.3f\n", rastScale.x, rastScale.y);
fprint(readme, "The extents of the area you extracted (referenced to the above CRS):");
fprintf(readme, "min x: %.6f, min y: %.6f, max x: %.6f, max y: %.6f\n", exRegion.Extents.x1, exRegion.Extents.y1, exRegion.Extents.x2, exRegion.Extents.y2);
$ifdef TEST
		print("Source raster info:");
		printf("Raster cols: %d, lines: %d\n", sourceRVC.GetTotalColumns(), sourceRVC.GetTotalRows() );
		printf("Source CRS: %s, %s\n", inCRS$, inCRS.Name);
		printf("Native cell size of source image: x = %.3f, y = %.3f\n", rastScale.x, rastScale.y);
		print("Extents of area to extract:");
		printf("minx: %.6f, miny: %.6f, max x: %.6f, max y: %.6f\n", exRegion.Extents.x1, exRegion.Extents.y1, exRegion.Extents.x2, exRegion.Extents.y2);
$endif
RastIn.Close();
##############################################################
###############  check extraction area  ###########################
# get extents of the image source as a region
class REGION2D sourceRegion;
err = sourceRVC.ComputeGeoreferenceRegion(sourceRegion);
if (err < 0)  {
	error$ = ReportError(_context.CurrentLineNum, err);
	notifyError(error$);
	}
fprint(logfile, "Source raster extents:");
fprintf(logfile, "Min x: %.4f, Min y: %.4f, Max x: %.4f, Max y: %.4f\n",
	sourceRegion.Extents.x1, sourceRegion.Extents.y1, sourceRegion.Extents.x2, sourceRegion.Extents.y2);
$ifdef TEST 
		print("Source raster extents:");
		printf("Min x: %.4f, Min y: %.4f, Max x: %.4f, Max y: %.4f\n",
			sourceRegion.Extents.x1, sourceRegion.Extents.y1, sourceRegion.Extents.x2, sourceRegion.Extents.y2);
$endif
# check that extraction region is containined within the source image extents region
if (sourceRegion.TestRegion(exRegion, "FullInside") )
	{
	fprint(logfile, "Source image contains extraction extents.");
	fprint(readme, "The source image contains the provided extraction extents.");
	$ifdef TEST
			print("Source image contains extraction extents.");
	$endif
	}
else
	{
	error$ = "Extraction areas is outside the extents of the source image. No extracted image was created.";
	fprintf(logfile, error$); 
	fprint(readme, "The desired extraction areas is outside the extents of the source image.");
	fprint(readme, "No extracted image could be created.");
	$ifdef TEST
			printf("Extraction areas is outside the extents of the source image. No extracted image was created.");
	$endif
	notifyError(error$);
	return;
	}
fclose(logfile);
logfile = fopen(logfileName$, "a");
########################################################################
###############  Determine resampling parameters  ###########################  
class RECT exRect;			# extents of the area to extract in output map coordinates
### check if input and output CRS are the same or different
if (CRS$ <> inCRS$)			# need to determine extraction extents in output CRS 
	{
	class STRING inCoordUnit$, outCoordUnit$;
	numeric unitScaleX = 1;					# scale factor from input to output raster scale units
	numeric unitScaleY = 1;
	# set output CRS from the EPSG code string provided in the job file
	outCRS.Assign(CRS$);
	# reproject the region of interest to the output CRS and get its extents
	exRegion.ConvertTo(outCRS);
	exRect = exRegion.Extents;
	fprintf(logfile, "The extraction area has been reprojected to %s\n", outCRS.Name);
	fprintf(readme, "The extraction area has been reprojected to %s\n", outCRS.Name);
	fprint(logfile, "Extents of the extraction area in the output coordinate reference system:");
	fprint(readme, "Extents of the extraction area in the output coordinate reference system:");
	fprintf(logfile, "Min x: %.4f, Min y: %.4f, Max x: %.4f, Max y: %.4f\n", exRect.x1, exRect.y1, exRect.x2, exRect.y2);
	fprintf(readme, "Min x: %.4f, Min y: %.4f, Max x: %.4f, Max y: %.4f\n", exRect.x1, exRect.y1, exRect.x2, exRect.y2);
	$ifdef TEST
			printf("The extraction area has been reprojected to %s\n", outCRS.Name);
			print("Extents of the extraction area in the output coordinate reference system:"); 
			printf("Min x: %.4f, Min y: %.4f, Max x: %.4f, Max y: %.4f\n", exRect.x1, exRect.y1, exRect.x2, exRect.y2);
	$endif
	# check input & output CRS projections to set cell sizes and scaling factors for computing output image dimensions
	if (inCRS.IsProjected() <> outCRS.IsProjected() )			# one CRS projected, one not; must adjust cell sizes
		{
		inCoordUnit$ = inCRS.CoordSys.GetAxis(1).Unit.GetSymbol();
		outCoordUnit$ = outCRS.CoordSys.GetAxis(1).Unit.GetSymbol();
		printf("Input CRS coordinate unit = %s, output CRS coordinate unit = %s\n", inCoordUnit$, outCoordUnit$);
		if (inCRS.IsProjected() == 0)			# input CRS is geographic,  output CRS is projected
			{
			print("Input CRS is Geographic, output CRS is projected");
			# recompute cell sizes (scales) in meters at center of extraction area of source image from its georeference
			georefSource.ComputeScale(centerPtImg, rastScale, 1);		# cell sizes in meters
			if (outCoordUnit$ <> "m")			# get distance unit conversion scale factor from input to output units
				{
				unitScaleX = GetUnitConvDist(inCoordUnit$, outCoordUnit$);
				unitScaleY = unitScaleX;
				}
			}
		else		# input CRS is projected, output is Geographic
			{
			print("Input CRS is projected, output is Geographic");
			class TRANS2D_MAPGEN maptrans;		# coordinate transformation from input to output CRS
			maptrans.InputCoordRefSys = inCRS;
			maptrans.OutputCoordRefSys = outCRS;
			# get extraction area center point in output geographic coordinates
			numeric projectedDistM;								# projected distance in meters
			centerPt = maptrans.ConvertPoint2DFwd(centerPt);			# reprojected center point of extraction area to geographic coords
			# compute unit scale factor from meters to degrees in east-west direction at center point by getting
			# projected distance in meters from center point to a point 0.01 degree to the east 
			projectedDistM = ProjDistanceToMeters(outCRS, centerPt.x, centerPt.y, centerPt.x + 0.01, centerPt.y);
			unitScaleX = 0.01 / projectedDistM;
			printf("x projected distance for .01 deg = %.2f m, unitScaleX = %.8f\n", projectedDistM, unitScaleX);
			# compute unit scale factor from meters to degrees in north-south direction at center point by getting
			# projected distance in meters from center point to a point 0.01 degree to the south
			projectedDistM = ProjDistanceToMeters(outCRS, centerPt.x, centerPt.y, centerPt.x, centerPt.y - 0.01);
			unitScaleY = 0.01 / projectedDistM;
			printf("y projected distance for .01 deg = %.2f m, unitScaleY = %.8f\n", projectedDistM, unitScaleY);
			if (inCoordUnit$ <> "m")			# apply unit conversion from input projected units to meters
				{
				unitScaleX = unitScaleX * GetUnitConvDist(inCoordUnit$, "m");
				unitScaleY = unitScaleY * GetUnitConvDist(inCoordUnit$, "m");
				}
			}
		}
	else			# input and output CRS are either both Geographic or  both projected
		{
		if (inCRS.IsProjected() )		# both CRSs are projected but could have different units
			{
			if (inCoordUnit$ <> outCoordUnit$)		# get distance unit conversion scale factor from input to output units
				{
				unitScaleX = GetUnitConvDist(inCoordUnit$, outCoordUnit$);
				unitScaleY = unitScaleX;
				}
			}
		}
	# Compute output image dimensions at native cell size (in output coordinate system)
	rastScale.y = rastScale.y * unitScaleY;
	rastScale.x = rastScale.x * unitScaleX;
	numLins = abs(round( (exRect.y2 - exRect.y1) / rastScale.y) );
	numCols = round( (exRect.x2 - exRect.x1) / rastScale.x);
	}	# end if (CRS$ <> inCRS$)
else			# CRS$ == inCRS$, no reprojection required
	{
	outCRS = inCRS;
	exRect = exRegion.Extents;
	# get source image coordinates for extraction area and compute raw extracted image dimensions
	class RECT imgRect;				# extents rectangle in image coordinates
	imgRect = transImgToCRS.ConvertRectInv(exRect);
	numLins = abs( round(imgRect.y2 - imgRect.y1) );
	numCols = round(imgRect.x2 - imgRect.x1);
	}
fprintf(logfile, "Output image dimensions at native cell size: %d lines, %d columns\n", numLins, numCols); 
$ifdef TEST
		printf("Output image dimensions at native cell size: %d lines, %d columns\n", numLins, numCols);
$endif
###################################################################
#### adjust output image dimensions and cell size if necessary to match output size limits
numeric maxAspect;		# maximum allowed aspect ratio at maximum size
class STRING size$;
if (limit$ == "1")
	{
	limitCells = 1000000;
#	maxAspect = 3.333;		# maximum 1024 cells in either dimension limits smaller dimension to 307 cells = 3.333:1
	size$ = "1 megapixel";
	}
else if (limit$ == "2")
	{
	limitCells = 2000000;
#	maxAspect = 2;					# maximum 2000 cells in either dimension limits smaller dimension to 1000 cells = 4:1
	size$ = "2 megapixels";
	}
else if (limit$ == "4")
	{
	limitCells = 4000000;
#	maxAspect = 4;					# maximum 4000 cells in either dimension limits smaller dimension to 1000 cells = 4:1
	size$ = "4 megapixels";
	}
fprintf(logfile, "Size limit: %s, %d\n", limit$, limitCells);
fprintf(readme, "Your selected size limit is %s, limiting the output image to %d cells\n", size$, limitCells); 
$ifdef TEST
		printf("Selected size limit: %s\n", limit$);
$endif
# reduce larger dimension of extraction area if maximum aspect ratio is exceeded
# maximum aspect ratio of extraction area = 4:1; reduce larger dimension if exceeded
#if (numLins > maxAspect * numCols) then
#		numLins = maxAspect * numCols;
#if (numCols > maxAspect * numLins) then
#	numCols = maxAspect * numLins;
# rescale image dimensions and cell size if total number of pixels exceeds limit
totalCells = numLins * numCols;
if (totalCells > limitCells)
	{
	rescale = sqrt(totalCells / limitCells);
	numLins = numLins / rescale;
	numCols = numCols / rescale;
 	rastScale.y = rastScale.y * rescale;
	rastScale.x = rastScale.x * rescale;
	fprintf(logfile, "Image resampled to coarser cell size: x = %.3f, y = %.3f\n", rastScale.x, rastScale.y); 
	fprintf(readme, "Your extracted image was resampled to a coarser cell size: x = %.3f, y = %.3f\n", rastScale.x, rastScale.y); 
	$ifdef TEST
			printf("Image resampled to coarser cell size: x = %.3f, y = %.3f\n", rastScale.x, rastScale.y);
	$endif
	}
else 
	{
	fprint(logfile, "Output image has cell size of input.");
	fprint(readme, "Your extracted image has the same cell size as the source image."); 
	$ifdef TEST
			print("Output image has cell size of input.");
	$endif
	}
# set output image dimensions for resampling filter to use
class IMAGE_PIPELINE_DIMENSIONS dimensions;
dimensions.SetTotalDimensions(numCols, numLins);
fprintf(logfile, "Size of target image: %d lines, %d columns\n", dimensions.GetTotalRows(), dimensions.GetTotalColumns() );
$ifdef TEST
		printf("Size of target image: %d lines, %d columns\n", dimensions.GetTotalRows(), dimensions.GetTotalColumns() );
$endif
fclose(logfile);
logfile = fopen(logfileName$, "a");
#################################################################
##############   Set up affine georeference for target image  ############# 
class TRANS2D_AFFINE transAffine;
transAffine.ApplyScale(rastScale.x,rastScale.y);
transAffine.ApplyOffset(exRect.x1,exRect.y2);
fprint(logfile, "Affine transformation parameters:");
fprintf(logfile, "Forward: x scale = %.6f, y scale = %.6f\n", transAffine.GetForwardScaleX(), transAffine.GetForwardScaleY() );
fprintf(logfile, "Forward shear = %.6f\n", transAffine.GetForwardShear() );
$ifdef TEST
		print("Affine transformation parameters:");
		printf("Forward: x scale = %.6f, y scale = %.6f\n", transAffine.GetForwardScaleX(), transAffine.GetForwardScaleY() );
		printf("Forward shear = %.6f\n", transAffine.GetForwardShear() );
$endif
class IMAGE_PIPELINE_GEOREFERENCE targetGeoref;
targetGeoref.SetTransAffine(outCRS, transAffine);
fclose(logfile);
logfile = fopen(logfileName$, "a");
###########################################################
#####   set up pipeline filter to resample the image    #######
class IMAGE_PIPELINE_FILTER_RESAMPLE filterResample(sourceRVC, dimensions, targetGeoref, "Nearest");
### INITIALIZE THE RESAMPLE FILTER
err = filterResample.Initialize();
	if ( err < 0)  {
		error$ = ReportError(_context.CurrentLineNum, err);
		notifyError(error$);
		}
fprintf(logfile, "resample or stereo filter initialized, returned %d\n", err);
fclose(logfile);
logfile = fopen(logfileName$, "a");
###########################################################
####  set up region source to use to mask area outside polygon  ######
class IMAGE_PIPELINE_SOURCE_REGION sourceReg(exRegion, filterResample);
err = sourceReg.Initialize();
if (err < 0)  {
	error$ = ReportError(_context.CurrentLineNum, err);
	notifyError(error$);
	}
fprintf(logfile, "sourceReg initialized, returned %d\n", err);
fclose(logfile);
logfile = fopen(logfileName$, "a");
#############################################################
####   set up filter to mask the image using the extraction region   #######
class IMAGE_PIPELINE_FILTER_MASKVALIDITY filterMask(filterResample, sourceReg);
err = filterMask.Initialize();
if (err < 0)  { 
	error$ = ReportError(_context.CurrentLineNum, err);
	notifyError(error$);
	}
fprintf(logfile, "filterMask initialized, returned %d\n", err);
fclose(logfile);
logfile = fopen(logfileName$, "a");
############################################################
#####  set up pipeline filter target and process  ##########
class IMAGE_PIPELINE_TARGET target;			# generic pipeline target, to be assigned specific type
switch (outFormat$)
	{
	case  "GeoJP2":
		outputFileExtension$ = ".jp2";
		outFilepath.SetExtension(outputFileExtension$);
		class IMAGE_PIPELINE_TARGET_J2K_SETTINGS j2k_settings;
		# set compression parameters
		if (lossyComp == 0) then
			j2k_settings.SetReversible(1);		# set lossless compression
		else
			j2k_settings.SetTargetRatio(lossyComp);	# set compression ratio
		class IMAGE_PIPELINE_TARGET_J2K targetJ2K(filterMask, outFilepath, j2k_settings, "ArcWorld");
		target = targetJ2K;
		break;
	case "GeoTIFF":
		outputFileExtension$ = ".tif";
		outFilepath.SetExtension(outputFileExtension$);
		class IMAGE_PIPELINE_TARGET_TIFF_SETTINGS tiff_settings;
		# set compression parameters
		if (lossyComp == 0) then
			tiff_settings.SetCompression("LZW");	# set lossless compression
		else
			tiff_settings.SetCompression("JPEG_DCT");		# set lossy compression
		class IMAGE_PIPELINE_TARGET_TIFF targetTiff(filterMask, outFilepath, "ArcWorld");
		targetTiff.SetParms(tiff_settings);
		target = targetTiff;
		break;
	case "JPEG":
		outputFileExtension$ = ".jpg";
		outFilepath.SetExtension(outputFileExtension$);
		class IMAGE_PIPELINE_TARGET_JPEG targetJPEG(filterMask, outFilepath, 80,  "ArcWorld");
		target = targetJPEG;
		break;
	case "PNG":
		outputFileExtension$ = ".png";
		outFilepath.SetExtension(outputFileExtension$);
		class IMAGE_PIPELINE_FILTER_NULLTOALPHA filterNullToAlpha(filterMask);
		err = filterNullToAlpha.Initialize();
		if (err < 0)  {
			error$ = ReportError(_context.CurrentLineNum, err);
			notifyError(error$);
			}
		fprintf(logfile, "filterNullToAlpha initialized, returned %d\n", err);
		fclose(logfile);
		logfile = fopen(logfileName$, "a");
		class IMAGE_PIPELINE_TARGET_PNG targetPNG(filterNullToAlpha, outFilepath, "ArcWorld");
		target = targetPNG;
		break;
	case "KMZ":
          if (false){
			outputFileExtension$ = ".png";
			outFilepath.SetExtension(outputFileExtension$);
			class IMAGE_PIPELINE_FILTER_NULLTOALPHA filterNullToAlpha(filterMask);
			err = filterNullToAlpha.Initialize();
			if (err < 0)  {
				 error$ = ReportError(_context.CurrentLineNum, err);
				notifyError(error$);
				}
			fprintf(logfile, "filterNullToAlpha initialized, returned %d\n", err);
			fclose(logfile);
			logfile = fopen(logfileName$, "a");
			class IMAGE_PIPELINE_TARGET_PNG targetPNG(filterNullToAlpha, outFilepath, "All");
			target = targetPNG;
			}
		else
			{
			outputFileExtension$ = ".jpg";
			outFilepath.SetExtension(outputFileExtension$);
			class IMAGE_PIPELINE_TARGET_JPEG targetJPEG(filterMask, outFilepath, 80,  "All");
			target = targetJPEG;
			}
		break;
	case "RVC":
		# create output Project File
		class RVC_OBJECT outRVC;
		outFilepath.SetExtension(".rvc");
		printf("Filepath = %s\n", outFilepath);
		outRVC.MakeFile(outFilepath, "Provisioning result"); 
		class RVC_OBJITEM outFileObjItem;			# ObjItem for the project file to serve as parent
		outFileObjItem = outRVC.GetObjItem();
		# RVC_OBJITEM for the new raster object
		class RVC_OBJITEM outRastObjItem;	# ObjItem for the output raster
		class RVC_DESCRIPTOR descriptor;
		descriptor.SetName(taskID$);				# use taskID for raster name
		descriptor.SetDescription("Provisioning result");
		outRastObjItem.CreateNew(outFileObjItem, "RASTER", descriptor);
		class IMAGE_PIPELINE_TARGET_RVC targetRVC(filterMask, 	outRastObjItem);
 		targetRVC.SetCompression("JP2", lossyComp);
		target = targetRVC;
		break;
	case "RVC_RGB":
		class RVC_OBJECT outRVC;
		outFilepath.SetExtension(".rvc");
		printf("Filepath = %s\n", outFilepath);
		outRVC.MakeFile(outFilepath, "Provisioning result"); 
		class RVC_OBJITEM outFileObjItem;			# ObjItem for the project file to serve as parent
		outFileObjItem = outRVC.GetObjItem();
		# Create RVC_OBJITEMs for the RGB output rasters
		class RVC_OBJITEM outRastObjItemList[];		# ObjItem list for the RGB output rasters
		class RVC_OBJITEM outRastObjItem;
		class RVC_DESCRIPTOR descriptor;
		for i = 1 to 3
			{
			descriptor.SetName( sprintf("%s_%d", taskID$, i) );
			descriptor.SetDescription( sprintf("Provisioning result component %d", i) );
			outRastObjItemList[i].CreateNew(outFileObjItem, "RASTER", descriptor);
			}
		class IMAGE_PIPELINE_TARGET_RVC targetRVC(filterMask, outRastObjItemList);
		targetRVC.SetCompression("JP2", lossyComp);
		target = targetRVC;
		break;
	}
fprintf(logfile, "Making %s\n",  outFormat$);
fprintf(logfile, "outFilepath = %s\n", outFilepath.GetPath() );
fclose(logfile);
logfile = fopen(logfileName$, "a");
if (outFormat$ == "RVC") then
	fprint(readme, "You selected MicroImages Project File format for your extract.");
else
	fprintf(readme, "You selected %s format for your extract\n.", outFormat$);
$ifdef TEST
	printf("Making %s\n", outFormat$);
$endif
  
err = target.Initialize();
printf("target initialization returned: %d\n", err);
if (err < 0)  {
	error$ = ReportError(_context.CurrentLineNum, err);
	notifyError(error$);
	}
err = target.Process();
if (err < 0)  {
	error$ = ReportError(_context.CurrentLineNum, err);
	notifyError(error$);
	}
######################################################################
if (outFormat$ == "RVC") 
	{
	outRVC.Close();						# close the output Project File so it can be zipped.
	fprint(logfile, "Closing output RVC file.");
	}
######################################################################
if (outFormat$ == "RVC_RGB")
	{
	class RVC_RASTER rastGray;
	class CONTRAST linearCon;
	targetRVC.GetObjItemList( outRastObjItemList );
	for i = 1 to 3
		{
		err = rastGray.Open( outRastObjItemList[i], "Write" );
		err = linearCon.ComputeStandard(rastGray, "linear", "LINEAR", "Default linear contrast table");
		rastGray.Close();
		}
	outRVC.Close();						# close the output Project File so it can be zipped.
	fprint(logfile, "Closing output RVC file.");
	}
######################################################################
fclose(readme);
##############################  Make ZIP File  ##########################
# Add output image file and delete raw file
class STRING imageName$ = outFilepath.GetName();
zip.AddFile(zipFilepath, outFilepath,  imageName$);
DeleteFile(outFilepath);
if (outFormat$ <> "PDF"  && outFormat$ <> "RVC" && outFormat$ <> "RVC_RGB" && outFormat$ <> "KMZ")
	{
	# Add PRJ file
	outFilepath.SetExtension(".prj");			# reset path to point to PRJ file
	class STRING name$ = outFilepath.GetName();
	zip.AddFile(zipFilepath, outFilepath,  name$);
	DeleteFile(outFilepath);
	# Add ArcWorld file; extension varies depending on output file format
	switch (outFormat$)
		{
		case "GeoJP2":
			outFilepath.SetExtension(".j2w");			# reset path to point to J2W file
		break;
		case "GeoTIFF":
			outFilepath.SetExtension(".tfw");			# reset path to point to TFW file
		break;
		case "JPEG":
			outFilepath.SetExtension(".jgw");			# reset path to point to J2W file
		break;
		case "PNG":
			outFilepath.SetExtension(".pgw");			# reset path to point to J2W file
		break;
		}
	name$ = outFilepath.GetName();
	zip.AddFile(zipFilepath, outFilepath,  name$);
	DeleteFile(outFilepath);
	}
if (outFormat$ == "KMZ")		# add KML file to KMZ file
	{
	outFilepath.SetExtension(".kml");
	name$ = outFilepath.GetName();
	zip.AddFile(zipFilepath, outFilepath,  name$);
	DeleteFile(outFilepath); 
	}
 else {
        # Add readme file
        class Filepath readmeFilepath(readmeFilename$);
        printf("readmeFilepath = %s + %s\n", readmeFilepath.GetPath(),  readmeFilepath.GetName() );
        zip.AddFile(zipFilepath, readmeFilepath, readmeFilepath.GetName() );
        DeleteFile(readmeFilename$);
        }
######################################################################
currentDT.SetCurrent();
$ifdef TEST
		printf("Job completed %s Central Standard Time\n\n", currentDT);
$endif
fprintf(logfile, "Job completed %s Central Standard Time\n", currentDT);
ExecuteProcess(finish$); 
fprintf(logfile, "##########################################################\n\n");
fclose(logfile);
#################################################################
#########  SEND EMAIL NOTIFICATION  ############################
class STRING mailFilename$;
class STRING emailbody1$, emailbody2$, emailbody3$;
class STRING emailbody$;
email.SetSubject(sprintf("Data Provisioning for Job %s Complete", taskID$) );
mailFilename$ = _context.ScriptDir + "/" + "emailbody1.txt";
emailbody1$ = TextFileReadFull(mailFilename$);
emailbody2$ = sprintf("Please download your sample geodata file now from http://www.geoprovisioning.com/jobs/%s/%s%s\n\n", taskID$, taskID$, ".zip");
mailFilename$ = _context.ScriptDir + "/" + "emailbody3.txt";
emailbody3$ = TextFileReadFull(mailFilename$);
emailbody$ = sprintf("%s%s%s", emailbody1$, emailbody2$, emailbody3$);
email.SetBody(emailbody$);
email.SetRecipient(userEmail$);
email.SendMail();