# ProvisionPoly_test.sml
# Image provisioning script for client providing polygon vertex coordinates.
# Version 5 June 2009 incorporating:
# - zoom and reprojection using the resampling filter
# - make zip file containing image, georeference, and readme files
# - log start and end processing times
# - RVC file output (currently not working)
# - TNTlite or Maximum Allowed size limits
# - Test mode
# - KMZ format with PNG
# - Mask out area outside the selection polygon
# - PNG format
# - PDF and RVC formats
# - fixes computation of output image dimensions to deal with change between projected / nonprojected CRS
# - stops making kml file for output formats other than KMZ
# - provides RVC Composite or RVC RGB Separate rasters for composite sources
# - provides anaglyph stereo rendering for orthoimage composites using NED10 as terrain
# - taskID appended to readme filename
# - use PDF.Finish() method to close PDF file so it can be zipped
# - error notification procedure added to e-mail staff on specific errors 5 June 09
# uncomment the following preprocessor command to run in test mode
# with hard-coded job parameters outside of the Job Control environment
$define TEST
############## Variables to be provided by the job file #################
#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: "GeoTIFF", "GeoJP2", "KMZ, or "RVC"
string CRS$; # EPSG code for the CRS selected by user for the output
string limit$; # size limit: possible values "tntlite" or "max"
string area$; # coordinate pairs for extraction polygon vertices
numeric lossyComp; # compression ratio for JP2 lossy compression, 0 for lossless
numeric anaglyph; # flag = 1 for rendering to anaglyph stereo
string anaglyphMode$; # color mode for anaglyph viewing;
##############################################################
################ Job parameters for test mode ###################
$ifdef TEST
# taskID$ = "50014";
# rvcPath$ = _context.ScriptDir + "/provision/ProvisionNE.rvc";
# rvcObject$ = "NE 2006 NC UTM14 1m 15to1";
# outFormat$ = "PNG";
# userEmail$ = "[email protected]";
CRS$ = "EPSG:26915";
limit$ = "max";
area$ = "182669 4737565,181360 4739802,183447 4742994,187566 4740539,186147 4736829,182669 4737565";
lossyComp = 15;
# anaglyph = 0;
anaglyphMode$ = "RedBlueMI";
$endif
##############################################################
############# VARIABLES FOR AUTOMATED TESTING ####################
# read from control file
string correctfilepath$; # path to file with correct output
string correctobjectpath$; # object path to correct object
string terrainPath$; # DEM variables for anaglyph mode
string terrainObj$;
#--------------------------------------
string outfilepath$ = CreateTempFileName(); # temporary output file for non-rvc target
string rvcfilepath$ = CreateTempFileName(); # temporary project file for rvc target
class FILEPATH outFilepath = outfilepath$;
class FILEPATH correctfile = correctfilepath$; # path to file with correct output
class RVC_OBJITEM outobjitem;
CreateProjectFile(rvcfilepath$, "Provision Poly Tests");
outobjitem.SetFilePath(rvcfilepath$);
outobjitem.SetObjectPath(correctobjectpath$);
class RVC_DESCRIPTOR desc;
desc.SetName(correctobjectpath$);
desc.SetDescription("Provision Poly Tests");
outobjitem.CreateNew(desc, "RASTER");
class RVC_RASTER RDest;
############# 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$;
###### 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$ = "php c:/wamp/www/error.php " + 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
proc ReportError(numeric linenum, numeric err) {
print("FAILED -line: %d, error: %d\n", linenum - 1, err);
Exit();
}
#############################################################
######################### Program ###########################
clear();
### set up to run php scripts to write info to MySQL database
#start$ = "php c:/wamp/www/start.php " + taskID$;
#finish$ = "php c:/wamp/www/finish.php " + taskID$;
#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.");
#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);
# printf("Anaglyph value = %d\n", anaglyph);
#$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);
#fprintf(logfile, "Anaglyph value = %d\n", anaglyph);
##################### 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.";
print("FAILED - source raster cannot be found.");
Exit();
# notifyError(error$);
}
##########################################################
########## set up pipeline source for source image #############
class IMAGE_PIPELINE_SOURCE_RVC sourceRVC(sourceObjItem);
err = sourceRVC.Initialize();
if (err < 0) {
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 #########################
# 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);
# set CRS for extraction region
exRegion.CoordRefSys = 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:");
#fprintf(logfile, "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);
#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("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);
#$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) {
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("FAILED - Extraction areas is outside the extents of the source image. No extracted image was created.");
Exit();
$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$ == "tntlite")
{
limitCells = 314368;
maxAspect = 3.333; # maximum 1024 cells in either dimension limits smaller dimension to 307 cells = 3.333:1
size$ = "TNTlite";
}
else if (limit$ == "max")
{
limitCells = 4000000;
maxAspect = 4; # maximum 4000 cells in either dimension limits smaller dimension to 1000 cells = 4:1
size$ = "Maximum Allowed";
}
#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");
#######################################################################
### if anaglyph, compute stereo parameters and set up stereo rendering filter, otherwise set up resample filter
class IMAGE_PIPELINE_FILTER rosFilter; # generic filter stage class for resample or stereo filter
if (anaglyph)
{
# string terrainPath$ = _context.ScriptDir + "/provision/ProvisionNE.rvc";
# string terrainObj$ = "NED10m";
class RVC_RASTER terrain;
OpenRaster(terrain, terrainPath$, terrainObj$);
# get region from the extents of the extraction area in the output CRS
class REGION2D exRegion2;
exRegion2.CoordRefSys = outCRS;
exRegion2.SetFromRect(exRect);
# $ifdef TEST
# print("Extents of area to crop terrain to in output CRS: ");
# printf("min x: %.6f, min y: %.6f, max x: %.6f, max y: %.6f\n\n", exRegion2.Extents.x1, exRegion2.Extents.y1, exRegion2.Extents.x2, exRegion2.Extents.y2);
# $endif
# PIPELINE SOURCE FOR TERRAIN RASTER
class IMAGE_PIPELINE_SOURCE_RVC sourceTerrain(terrain.GetObjItem() );
err = sourceTerrain.Initialize();
if ( err < 0) {
ReportError(_context.CurrentLineNum, err);
#notifyError(error$);
}
else
{
# fprint(logfile, "Set up anaglyph stereo rendering.");
# fprint(logfile, "Terrain raster initialized.");
# $ifdef TEST
# print("Set up anaglyph stereo rendering.");
# print("Terrain raster initialized.");
# $endif
}
# get georeference and CRS from the terrain raster
class IMAGE_PIPELINE_GEOREFERENCE terrainGeoref;
terrainGeoref = sourceTerrain.GetGeoreference();
class SR_COORDREFSYS terrainCRS;
terrainCRS = terrainGeoref.GetCRS();
string terrainCRS$ = "EPSG:" + terrainCRS.GetID("EPSG");
# fprintf(logfile, "terrainCRS = %s, %s\n", terrainCRS$, terrainCRS.Name);
# $ifdef TEST
# printf("terrainCRS = %s, %s\n", terrainCRS$, terrainCRS.Name);
# $endif
# if terrain CRS is different from the selected output CRS, convert the extraction
# region to the output CRS
if (terrainCRS$ <> CRS$)
{
exRegion2.ConvertTo(terrainCRS);
# $ifdef TEST
# printf("terrain extraction region converted to: EPSG:%s, %s\n", exRegion2.CoordRefSys.GetID("EPSG"), exRegion2.CoordRefSys.Name );
# $endif
}
class RECT exRECT2; # rectangular extents of the extraction region in the terrain CRS
exRECT2 = exRegion2.Extents;
# $ifdef TEST
# print("Extents of area to crop terrain to in terrain CRS: ");
# printf("min x: %.6f, min y: %.6f, max x: %.6f, max y: %.6f\n\n", exRECT2.x1, exRECT2.y1, exRECT2.x2, exRECT2.y2);
# $endif
# get coordinate transformation from terrain image coordinates to its defined CRS
class TRANS2D_MAPGEN transTerrain;
transTerrain = terrainGeoref.GetTransGen();
# transform extraction extents rect to terrain object coordinates
exRECT2 = transTerrain.ConvertRectInv(exRECT2);
# $ifdef TEST
# print("Extraction area in terrain object coordinates:");
# printf("min x = %.1f, max x = %.1f\n", exRECT2.x1, exRECT2.x2);
# printf("miny = %.1f, max y = %.1f\n", exRECT2.y1, exRECT2.y2);
# $endif
# fprint(logfile, "Extraction area in terrain object coordinates:");
# fprintf(logfile, "min x = %.1f, max x = %.1f\n", exRECT2.x1, exRECT2.x2);
# fprintf(logfile, "miny = %.1f, max y = %.1f\n", exRECT2.y1, exRECT2.y2);
# PIPELINE FILTER TO CROP TERRAIN RASTER TO TEMPORARY RASTER TO GET
# LOCAL ELEVATION PARAMETERS FOR STEREO RENDERING
class IMAGE_PIPELINE_FILTER_CROP filterCrop(sourceTerrain, exRECT2);
err = filterCrop.Initialize();
if ( err < 0) {
ReportError(_context.CurrentLineNum, err);
#notifyError(error$);
}
# PIPELINE TARGET FOR CROPPED TERRAIN
class RVC_OBJECT tempFileObj;
tempFileObj.MakeTempFile(1);
class RVC_OBJITEM terrainCropObjItem;
class RVC_DESCRIPTOR terrainCropDescriptor;
terrainCropDescriptor.SetName("tempTerrain");
terrainCropObjItem.CreateNew(tempFileObj.GetObjItem(), "RASTER", terrainCropDescriptor);
class IMAGE_PIPELINE_TARGET_RVC targetTerrainCrop(filterCrop, terrainCropObjItem );
err = targetTerrainCrop.Initialize();
if ( err < 0) {
ReportError(_context.CurrentLineNum, err);
#notifyError(error$);
}
# EXECUTE PIPELINE FOR TEMPORARY CROPPED TERRAIN
err = targetTerrainCrop.Process();
if (err < 0) {
ReportError(_context.CurrentLineNum, err);
#notifyError(error$);
}
else
{
# fprint(logfile, "Cropping terrain to temporary raster...");
# $ifdef TEST
# print("Cropping terrain to temporary raster...");
# $endif
}
# Get elevation values from the cropped terrain
class RVC_RASTER terrainCrop;
terrainCrop.Open(targetTerrainCrop.GetObjItem(), "Read");
numeric elevMin; # minimum elevation
numeric elevRange; # elevation range
elevMin = GlobalMin(terrainCrop);
elevRange = GlobalMax(terrainCrop) - elevMin;
# compute stereo rendering parameters from elevation values
numeric relativeDepthScale = 30;
numeric ratioVE = 20;
numeric maxReliefExaggeration = 5;
numeric displacementZeroValue = elevMin + 0.8 * elevRange;
numeric stereoDisplacementScale = 1.5 * relativeDepthScale / elevRange;
numeric cellsize = sqrt( abs(rastScale.x) * abs(rastScale.y) );
numeric reliefExaggeration;
reliefExaggeration = cellsize * relativeDepthScale * ratioVE / elevRange;
# fprintf(logfile, "Cropped terrain min elev = %.1f, elev range = %.1f\n", elevMin, elevRange);
# fprintf(logfile, "Computed stereoDisplacementScale = %.2f\n", stereoDisplacementScale);
# fprintf(logfile, "computed displacementZeroValue = %.2f\n", displacementZeroValue);
# fprintf(logfile, "computed reliefExaggeration = %.1f\n", reliefExaggeration);
# $ifdef TEST
# printf("Cropped terrain min elev = %.1f, elev range = %.1f\n", elevMin, elevRange);
# printf("Computed stereoDisplacementScale = %.2f\n", stereoDisplacementScale);
# printf("computed displacementZeroValue = %.2f\n", displacementZeroValue);
# printf("computed reliefExaggeration = %.1f\n", reliefExaggeration);
# $endif
if (reliefExaggeration > maxReliefExaggeration)
{
stereoDisplacementScale = maxReliefExaggeration / reliefExaggeration;
# fprintf(logfile, "adjusted stereoDisplacementScale = %.1f\n", stereoDisplacementScale);
# $ifdef TEST
# printf("adjusted stereoDisplacementScale = %.1f\n", stereoDisplacementScale);
# $endif
}
# close the cropped terrain raster.
terrainCrop.Close();
# set mode for stereo rendering
class STEREOPARM stereoparm;
stereoparm.SetMethod("Anaglyph");
stereoparm.SetAnaglyphMode(anaglyphMode$);
# PIPELINE STEREO RENDERING FILTER
class IMAGE_PIPELINE_FILTER_STEREO stereoFilt(sourceRVC, sourceTerrain,
stereoDisplacementScale, displacementZeroValue, dimensions, targetGeoref, "Nearest", 0);
stereoFilt.SetAccuracy(1);
stereoFilt.SetStereoParm(stereoparm);
rosFilter = stereoFilt;
} # end if (anaglyph)
else
{
###########################################################
##### set up pipeline filter to resample the image #######
class IMAGE_PIPELINE_FILTER_RESAMPLE filterResample(sourceRVC, dimensions, targetGeoref, "Nearest");
rosFilter = filterResample;
}
### INITIALIZE THE RESAMPLE/STEREO RENDERING FILTER
err = rosFilter.Initialize();
if ( err < 0) {
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, rosFilter);
err = sourceReg.Initialize();
if (err < 0) {
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(rosFilter, sourceReg);
err = filterMask.Initialize();
if (err < 0) {
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) {
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 (lossyComp == 0)
# {
outputFileExtension$ = ".png";
outFilepath.SetExtension(outputFileExtension$);
class IMAGE_PIPELINE_FILTER_NULLTOALPHA filterNullToAlpha(filterMask);
err = filterNullToAlpha.Initialize();
if (err < 0) {
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);
class IMAGE_PIPELINE_TARGET_RVC targetRVC(filterMask, outobjitem);
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;
case "PDF":
# create temporary Project File for temp raster
# class RVC_OBJECT outRVC;
# outFilepath.SetExtension(".rvc");
# printf("Filepath = %s\n", outFilepath);
# outRVC.MakeTempFile(1); # delete on close
# 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);
# 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) {
ReportError(_context.CurrentLineNum, err);
#notifyError(error$);
}
err = target.Process();
if (err < 0) {
ReportError(_context.CurrentLineNum, err);
#notifyError(error$);
}
################# Make PDF File from RVC Target #########################
#if (outFormat$ == "PDF")
# {
# set up a hardcopy layout with group containing the output raster and render to PDF
# fprint(logfile, "Setting up layout to render to PDF.");
# class PRINTPARMS printparms; # page size and scale settings for hardcopy layou t
# numeric inTOm = GetUnitConvDist("in", "m"); # scale factor from inches to meters
# numeric imgSize, printSize;
# set up 8.5 x 11 in page (printable area 7.5 x 10 in, aspect ratio 1.333)
# if (numLins >= numCols) # set page to portrait orientation
# {
# printparms.FullWidth = 8.5; # width in inches including margins
# printparms.FullHeight = 11; # height in inches including margins
# printparms.PrintableWidth = 7.5;
# printparms.PrintableHeight = 10;
# if (numLins > 1.333 * numCols) # set layout scale from image/page height
# {
# imgSize = numLins * abs(rastScale.y); # image ground height in meters
# printSize = 10 * inTOm; # allowed printed size in meters
# }
# else # set layout scale from image/page width
# {
# imgSize = numCols * rastScale.x; # image ground width in meters
# printSize = 7.5 * inTOm;
# }
# }
# else # set page to landscape orientation
# {
# printparms.FullWidth = 11;
# printparms.FullHeight = 8.5;
# printparms.PrintableWidth = 10;
# printparms.PrintableHeight = 7.5;
# if (numCols > 1.333 * numLins) # set layout scale from image/page width
# {
# imgSize = numCols * rastScale.x;
# printSize = 10 * inTOm;
# }
# else # set layout scale from image/page height
# {
# imgSize = numLins * abs(rastScale.y);
# printSize = 7.5 * inTOm;
# }
# }
# printparms.MapScale = imgSize / printSize;
# class GRE_LAYOUT pdfLayout; # set up the hardcopy layout
# pdfLayout.Create("pdfLayout", 1);
# pdfLayout.Hardcopy = printparms;
# fprintf(logfile, "Layout %s width = %d in, height = %d in\n", pdfLayout.Name, printparms.FullWidth, printparms.FullHeight);
# fprintf(logfile, "Image ground size = %d m, printSize = %.4f m\n", imgSize, printSize);
# fprintf(logfile, "Printparms map scale = %d\n", printparms.MapScale);
# fprintf(logfile, "Layout map scale = %d\n", pdfLayout.MapScale);
# fclose(logfile);
# logfile = fopen(logfileName$, "a");
# $ifdef TEST
# printf("Layout %s width = %d, height = %d\n", pdfLayout.Name, printparms.FullWidth, printparms.FullHeight);
# printf("Image ground size = %d, printSize = %.4f\n", imgSize, printSize);
# printf("Printparms map scale = %d\n", printparms.MapScale);
# printf("Layout map scale = %d\n", pdfLayout.MapScale);
# $endif
# class GRE_GROUP imgGroup; # set up group to contain the image
# imgGroup = GroupCreate("ImageGroup", pdfLayout);
# printf("Name of group to render image: %s\n", imgGroup.Name);
# printf("Name of layout the group belongs to: %s\n", imgGroup.Layout.Name);
# class RVC_RASTER outRast; # handle for the output RVC raster created by the pipeline
# outRastObjItem = targetRVC.GetObjItem();
# err = outRast.Open(outRastObjItem, "Read");
# printf("outRast.Open() returned: %d\n", err);
# GroupQuickAddRasterVar(imgGroup, outRast); # add output raster to group
# fprintf(logfile, "Group name = %s\n", imgGroup.Name);
# $ifdef TEST
# printf("Group name = %s\n", imgGroup.Name);
# $endif
# class PDF pdf; # set up PDF file
# pdf.SetLayout(pdfLayout);
# outFilepath.SetExtension(".pdf");
# pdf.SetPath(outFilepath);
# pdf.Write();
# pdf.Finish();
# } # end if PDF
######################################################################
if (outFormat$ == "RVC")
{
# outRVC.Close(); # close the output Project File so it can be zipped.
# fprint(logfile, "Closing output RVC file.");
class RVC_OBJITEM targetobjitem;
targetobjitem.SetFilePath(correctfilepath$);
targetobjitem.SetObjectPath(correctobjectpath$);
if (!targetobjitem.IsExisting()) {#no known good exists
print("FAILED - No correct image exists for ", outFormat$);
CopyObject(rvcfilepath$, 1, correctfilepath$);
}
OpenRaster(RDest, rvcfilepath$, correctobjectpath$);
}
else {
if (!correctfile.IsFile()) { #no known good exists
print("FAILED - No correct image exists for ", outFormat$);
CopyFile(outfilepath$, correctfilepath$);
}
RDest.OpenByName(outFilepath, outFilepath.GetNameOnly());
}
######################################################################
#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")
# {
# 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;
# case "KMZ":
# if (lossyComp == 0) then
# outFilepath.SetExtension(".pgw");
# else
# outFilepath.SetExtension(".jgw"); # reset path to point to TFW file
# }
# 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);
# }
# 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);
#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();
#ExecuteProcess(finish$);
outFilepath.Delete();
print("Process Completed");
Download