CropTiffsWithShapefiles.sml

  Download

More scripts: Raster

Syntax Highlighing:

comments, key words, predefined symbols, class members & methods, functions & classes
            
## CropTiffsWithShapefiles.sml
#		Script to resample georeferenced airphoto TIFF files to
#		NAD83 / New Jersey State Plane Feet and crop to
#		a rectangular area in that coordinate system.  The
#		rectangular area for each  photo TIFF is defined by a 
#		rectangle in an accompanying shapefile.
#		The output for each input TIFF is a resampled and cropped
#		TIFF file (uncompressed, with no GeoTIFF tags) with
#		Arc World (.tfw) and .prj files.  A log file (text) is
#		created in the same directory as the output files to
#		record the processing parameters and result for each
#		TIFF file.
#		The script creates and opens a dialog window to allow the user
#		to create input and output directories and to provide status
#		information on the processing.
#		Assumptions and Conditions:
#		1) The user selects a single input directory containing both
#				TIFF and matching shape files
#		2) Input TIFF filenames begin with a unique 4-digit photo number.
#		3) The shapefile that matches a particular input TIFF has a filename that 
#				begins with the same 4-digit identifier as that TIFF file.
#		4) TIFF files that do not have a matching shape file in the input directory
#				are skipped during processing and reported in the script dialog and log file.
#		5) A TIFF file whose rectangle in the matching shape file is not fully contained
#				within the map extents of the TIFF file is skipped during processing and
#				reported in the log file.
#		5) Extra nonmatching shape files in the input directory are ignored.
#		6) All output TIFF files are placed in a single directory selected by the user.
#		7) Georeferenced input TIFFs need not be georeferenced to the same
#				coordinate reference system.  (i.e, TIFFs georeferenced to New Jersey
#				State Plane feet or meters will both work).
#		 Version 26 March 2009
#		 Requires TNT 2008:74 dated 30 December 2008 or later 
###############   Global Variables   #################
class STRING inDir$, outDir$;
numeric err;																		# error flag
class STRINGLIST inTiffList, inShapeList;			# lists of input TIFF and shape files
class STRING logfilename$;										# name of log file
class FILE logfile;
class GUI_DLG dlgwin;
class GUI_CTRL_PUSHBUTTON outDirBtn;
class GUI_CTRL_EDIT_STRING indirtext, outdirtext;
class GUI_CTRL_EDIT_NUMBER numTiff, numShape, numProc;
class GUI_CTRL_LISTBOX TIFFlistbox, Shapelistbox, Proclistbox;
class GUI_CTRL_EDIT_STRING statustext;
class GUI_CTRL_PUSHBUTTON runBtn, cancelBtn, closeBtn;
class STATUSDIALOG statusD;
class STATUSCONTEXT statusC;
###############   User-defined Procedures   #############
##########################################
# error checking procedure
proc ReportError(numeric linenum, numeric err) {
	printf("FAILED -line: %d, error: %d\n", linenum - 1, err);
	PopupError(err); 
	}
###########################################
### function to find matching shapefile and return its path string
func string findShapefile(class STRINGLIST list, class STRING id$)
	{
	local numeric j;
	local numeric shapefound = 0;
	local string foundShapename$;
	local class FILEPATH shapeName;			# filepath for the current shapefile
	local class STRING shapename$;			# shape filename
	local class STRING test$;							# first nine characters of the shape filename
	for j = 1 to list.GetNumItems()
		{
		shapeName = list[j-1];
		shapename$ = shapeName.GetName();
		test$ = shapename$.substr(0, 9);
		if (test$ == id$)
			{
			shapefound = 1;
			foundShapename$ = shapename$;
			}
		}
		if (shapefound)
			return foundShapename$;
		else
			return "failed";
	}	 # end findShapefile
############################
### procedure to get the output directory
proc GetOutDirectory()
	{
	outDir$ = GetDirectory("", "Select output directory:");
	class FILEPATH outDirPath(outDir$);
	outdirtext.SetValueStr(outDir$);		# write output directory name to dialog
	class DATETIME currentDT;
	currentDT.SetCurrent();
	# open log file for this run
	logfilename$ = sprintf("%s/%s.log", outDir$, outDirPath.GetNameOnly() );
	print(logfilename$);
	logfile = fopen(logfilename$, "a");
	fprint(logfile, "Log for TIFF file extraction and resampling.");
	fprintf(logfile, "\nProcessing initiated %s Local Time\n", currentDT);
	fprintf(logfile, "Input directory = %s\n", inDir$);
	fprintf(logfile, "Output directory = %s\n", outDir$);
	fprintf(logfile, "Number of tiff files found = %d\n", inTiffList.GetNumItems() );
	fprintf(logfile, "Number of shape files found = %d\n", inShapeList.GetNumItems() );
	# enable Run button
	runBtn.SetEnabled(1);
	}
############################
### procedure to get the input directory
proc GetInDirectory()
	{
	local numeric i, j;
	local numeric shapecount;
	inDir$ = GetDirectory("", "Select directory with TIFF and SHAPE files:");
	class FILEPATH inDirPath(inDir$);
	indirtext.SetValueStr(inDir$);				# write input directory name to dialog
	# clear dialog listboxes if they are still populated from a previous run
	if (TIFFlistbox.GetCount() > 0) then
		TIFFlistbox.DeleteAllItems();
	if (Shapelistbox.GetCount() > 0) then
			Shapelistbox.DeleteAllItems();
	if (Proclistbox.GetCount() > 0)
			{ 
			Proclistbox.DeleteAllItems();
			numProc.SetValueNum(0);
			}
	outDir$ = "";
	outdirtext.SetValueStr("");		# clear output directory text field 
	statustext.SetValueStr("");		# clear status message text field
	local class STRING tiffName$;			# filename of current TIFF file
	local class STRING id$;							# id$ = first nine characters of TIFF filename
	local class STRING shapefileName$;
	inTiffList = inDirPath.GetFileList("*.tif");
	inShapeList = inDirPath.GetFileList("*.shp");
	for i = 1 to 	inTiffList.GetNumItems()
		{
		TIFFlistbox.AddItem( inTiffList[i-1] );		# add TIFF filename to listbox in dialog
		# get first nine characters of TIFF filename as id to check for matching shapefile
		tiffName$ = inTiffList[i-1];
		id$ = tiffName$.substr(0, 9);
		shapefileName$ = findShapefile(inShapeList, id$);
		if (shapefileName$ <> "failed")		# matching shape file found
			{
			Shapelistbox.AddItem(shapefileName$);
			++shapecount;
			}
		else
			{
			statustext.SetValueStr( sprintf("No matching shapefile found for %s\n",   tiffName$) );
			}
		}
	# write number of files found to the dialog
	numTiff.SetValueNum(inTiffList.GetNumItems() );
	numShape.SetValueNum(shapecount);
	# enable Output Directory and Cancel buttons
	outDirBtn.SetEnabled(1);
	cancelBtn.SetEnabled(1);
	}
############################
### procedure called by the Close button
proc onClose()
	{
	fclose(logfile);
	dlgwin.Close(0);
	Exit();
	}
############################
### procedure called by the Close button
proc onCancel()
	{
	fclose(logfile);
	dlgwin.Close(0);
	Exit();
	}
############################
### procedure called by the Run button
proc onRun()
	{
	local numeric success;
	statusD.Create();
	statusD.SetSelfAsDefault();
	statusC = statusD.CreateContext();
	### set output Coordinate Reference System 
	### NAD83 / New Jersey State Plain ftUS
	class SR_DATUM datum;
	datum.Assign("2347");
	class SR_COORDSYS coordsys;
	coordsys.Assign("1202");
	class SR_COORDOPDEF projection;
	projection.Assign("25196");
	class SR_COORDREFSYS crsNJSPFT;
	crsNJSPFT.Create(coordsys, datum, projection);
	fprintf(logfile, "Output coordinate reference system = %s\n", crsNJSPFT.Name);
	###########################################
	### loop through the list of input TIFF files to perform extraction
	local numeric i;   # counter
	statusC.BarInit(inShapeList.GetNumItems(), 0);
	for i = 1 to inTiffList.GetNumItems()
		{
		fprintf(logfile, "\nTiff file %d: %s\n", i, inTiffList[i-1] );
		string shapepath$;			# full path string for the matching shapefile
		class STRING tiffName$ = inTiffList[i-1];	# filename of current TIFF file
		class STRING id$ = tiffName$.substr(0, 9);				# id$ = first four characters of TIFF filename
		# call custom function to find shapefile name with first 4 characters matching those of TIFF file
		shapepath$ = findShapefile(inShapeList, id$);
		if (shapepath$ == "failed")		# no matching shape file found
			{
			fprintf(logfile, "No matching shapefile found for %s\n",   tiffName$);
			statusC.Message = sprintf("No matching shapefile found for %s\n",   tiffName$);
			statustext.SetValueStr( sprintf("No matching shapefile found for %s\n",   tiffName$) );
			}
		else			# a shape file matches the current TIFF file
			{
			statusC.Message = sprintf("Processing %s...", tiffName$);
			statusC.BarIncrement(1, 0);
			statustext.SetValueStr(sprintf("Processing %s...", tiffName$) );
			++ success;
			# construct full filepath for matching shapefile 
			class FILEPATH shapeFilepath;
			shapeFilepath = inDirPath;
			shapeFilepath.Append(shapepath$);
			fprintf(logfile, "matching shapefile = %s\n", shapeFilepath);
			# construct full path for current TIFF file
			class FILEPATH tiffFilepath;
			tiffFilepath = inDirPath;
			tiffFilepath.Append(tiffName$);
			# Set up pipeline source for TIFF file
			class IMAGE_PIPELINE_SOURCE_TIFF sourceTIFF(tiffFilepath);
			err = sourceTIFF.Initialize();
			# get pipeline georeference from the source image
			class IMAGE_PIPELINE_GEOREFERENCE georefSource;		# georeference from source
			georefSource = sourceTIFF.GetGeoreference();
			# get coordinate reference system from source TIFF image
			class SR_COORDREFSYS crsSource;
			crsSource = georefSource.GetCRS(); 
			fprintf(logfile, "source TIFF coordinate reference system = %s\n", crsSource.Name);
			# get extents of the source image as a region
			class REGION2D sourceRegion;
			sourceTIFF.ComputeGeoreferenceRegion(sourceRegion);
			# get centroid of the extents region to use to compute source cell size
			class POINT2D centerPt;
			centerPt = sourceRegion.GetCentroid();
			fprintf(logfile, "Center point of TIFF: x = %.2f, y = %.2f\n", centerPt.x, centerPt.y);
			# get coordinate transformation from source image to its defined CRS
			class TRANS2D_MAPGEN transImgToCRS = georefSource.GetTransGen();
			# compute center point of source TIFF image in image coordinates
			class POINT2D centerPtImg;
			centerPtImg = transImgToCRS.ConvertPoint2DInv(centerPt);
			fprintf(logfile, "Center point of TIFF in image coordinates: x = %.1f, y = %.1f\n", centerPtImg.x, centerPtImg.y); 
			# compute cell sizes (scales) in meters at center of TIFF image from its georeference
			class POINT2D rastScale;
			georefSource.ComputeScale(centerPtImg, rastScale, 1);
			fprintf(logfile, "TIFF cell size (m): x = %.1f, y = %.1f\n", rastScale.x, rastScale.y);
			# convert cell sizes to US survey feet to use to find dimensions of output image from
			# the extents of the rectangle in the shape file
			numeric metersToFeet = GetUnitConvDist("m", "ft");
			rastScale.x = rastScale.x * metersToFeet;
			rastScale.y = rastScale.y * metersToFeet;
			fprintf(logfile, "TIFF cell size (ft): x = %.1f, y = %.1f\n", rastScale.x, rastScale.y);
			# convert source region to New Jersey State Plane feet if necessary
			if (sourceRegion.CoordRefSys.Name <> crsNJSPFT.Name)
				{
				fprintf(logfile, "Converting TIFF extents region %s to %s to test containment.\n", tiffName$, crsNJSPFT.Name);
				sourceRegion.ConvertTo(crsNJSPFT);
				}
			# open the shapefile
			class RVC_SHAPE shape;
			err = shape.OpenByName(shapeFilepath, shapeFilepath.GetName(), "Read");
			# get the georeference from the shapefile
			class RVC_GEOREFERENCE shapeGeoref;
			shape.GetDefaultGeoref(shapeGeoref);
			fprintf(logfile, "Shapefile coordinate reference system = %s\n", shapeGeoref.GetCoordRefSys().Name);
			# get the rectangle from the shapefile as a region and set its CRS
			class REGION2D clipRegion;
			clipRegion = shape.ReadPolygons(0);
			clipRegion.CoordRefSys =shapeGeoref.GetCoordRefSys();
			# convert extraction region to New Jersey State Plane feet if necessary
			if (clipRegion.CoordRefSys <> crsNJSPFT.Name)
				{
				fprintf(logfile, "Converting clip region for %s to %s\n", tiffName$, crsNJSPFT.Name);
				clipRegion.ConvertTo(crsNJSPFT);
				}
			### check extraction area
			if (sourceRegion.TestRegion(clipRegion, "FullInside") )
				{
				fprintf(logfile, "Source image %s contains the extraction area.\n", tiffName$);
				fprintf(logfile, "clipRegion minX = %.4f, maxX = %.4f\n", clipRegion.Extents.x1, clipRegion.Extents.x2);
				fprintf(logfile, "clipRegion minY = %.4f, maxY = %.4f\n", clipRegion.Extents.y1, clipRegion.Extents.y2);
				# compute image extract dimensions in raster lines and columns
				numeric numLins, numCols;
				numLins = abs(round( (clipRegion.Extents.y2 - clipRegion.Extents.y1) / rastScale.y));
				numCols = round( (clipRegion.Extents.x2 - clipRegion.Extents.x1) / rastScale.x);
				fprintf(logfile, "Image dimensions of extraction area: %d lines, %d columns\n", numLins, numCols);
				class IMAGE_PIPELINE_DIMENSIONS dimensions;
				dimensions.SetTotalDimensions(numCols, numLins);
				# set up affine georeference for target image
				class TRANS2D_AFFINE transAffine;
				transAffine.ApplyScale(rastScale.x, rastScale.y);
				transAffine.ApplyOffset(clipRegion.Extents.x1, clipRegion.Extents.y2);
				class IMAGE_PIPELINE_GEOREFERENCE targetGeoref;
				targetGeoref.SetTransAffine(crsNJSPFT, transAffine);
				# set up pipeline filter to resample & extract the image
				class IMAGE_PIPELINE_FILTER_RESAMPLE filterResample(sourceTIFF, dimensions, targetGeoref, "Nearest");
				err = filterResample.Initialize();
				if (err < 0) then ReportError(_context.CurrentLineNum, err);
				# construct filepath for output TIFF file
				class FILEPATH outFilepath(outDir$);
				outFilepath.Append(id$ + "crp.tif");
				fprintf(logfile, "Filepath for the output TIFF file = %s\n", outFilepath);
				# set up pipeline filter target and process
				class IMAGE_PIPELINE_TARGET_TIFF tiffTarget(filterResample, outFilepath, "ArcWorld");
				class IMAGE_PIPELINE_TARGET_TIFF_SETTINGS settings;
				settings.SetGeoTag("None");			# no GeoTiff tags
				tiffTarget.SetParms(settings);
				err = tiffTarget.Initialize();
				if (err < 0) then ReportError(_context.CurrentLineNum, err);
				err = tiffTarget.Process();
				if (err < 0) then ReportError(_context.CurrentLineNum, err);
				fprintf(logfile, "File %s successfully extracted.\n", tiffName$);
				numProc.SetValueNum(success);
				Proclistbox.AddItem(	outFilepath.GetName() );
				} # end: if (sourceRegion.TestRegion(clipRegion, "FullInside") )
			else
				{
				fprintf(logfile, "Extraction area for %s is outside the extents of the source image.  No extracted image created.\n", tiffName$); 
				}
			} # end: else a shapefile matches current tiff file
		} # end: 	for i = 1 to inTiffList.GetNumItems()
	statusD.Destroy();
	runBtn.SetEnabled(0);
	cancelBtn.SetEnabled(0);
	statustext.SetValueStr("Processing complete.");
	} # end of onRun()
##################################################
#################   Main Program   #####################
clear();
class STRING xml$='<?xml version="1.0"?>
<!DOCTYPE root SYSTEM "smlforms.dtd">
<root>
	<dialog id="dlgxtiff" Title="Crop/Resample TIFFs Using Shapefiles" Buttons="">
		<pane Orientation="Horizontal" HorizResize="Fixed">
			<pushbutton Name=" Input Directory " OnPressed="GetInDirectory()"/>
			<edittext id="indirtext" width="30" ReadOnly="true"/>
		</pane>
		<pane Orientation="Horizontal" HorizResize="Fixed">
			<pushbutton id="outDirBtn" Name=" Output Directory " Enabled="false" OnPressed="GetOutDirectory()"/>
			<edittext id="outdirtext" Width="30" ReadOnly="true"/>
		</pane>
		<groupbox ExtraBorder="3">
			<pane Orientation="Horizontal" HorizResize="Fixed">
				<pane Orientation="Vertical" HorizResize="Fixed" VertResize="Fixed">
					<listbox id="TIFFlistbox" Width="10"/>
					<pane Orientation="Horizontal" HorizResize="Fixed">
						<editnumber id="numTiff" Width="2" Precision="0" BlankZero="true" ReadOnly="true"/>
						<label>  TIFF files   </label>
					</pane>
				</pane>
				<pane Orientation="Vertical" HorizResize="Fixed" VertResize="Fixed">
					<listbox id="Shapelistbox" Width="15"/>
					<pane Orientation="Horizontal" HorizResize="Fixed">
						<editnumber id="numShape" Width="2" Precision="0" BlankZero="true" ReadOnly="true"/>
						<label>  matching shapefiles   </label>
					</pane>
				</pane>
				<pane Orientation="Vertical" HorizResize="Fixed" VertResize="Fixed">
					<listbox id="Proclistbox" Width="10"/>
					<pane Orientation="Horizontal" HorizResize="Fixed">
						<editnumber id="numProc" Width="2" Precision="0" BlankZero="true" ReadOnly="true"/>
						<label>  successfully processed</label>
					</pane>
				</pane>
			</pane>
		</groupbox>
		<pane Orientation="Horizontal" HorizAlign="left">
			<edittext id="statustext" Width="20" ReadOnly="true"/>
		</pane>
		<pane Orientation="Horizontal" HorizAlign="right">
			<pushbutton id="runBtn" Name=" Run " Enabled="false" OnPressed="onRun()"/>
			<pushbutton id="cancelBtn" Name=" Cancel " Enabled="false" OnPressed="onCancel()"/>
			<pushbutton id="closeBtn" Name=" Close " OnPressed="onClose()"/>
		</pane>
	</dialog>
</root>';
### parse XML text for the dialog into memory
class XMLDOC dlgdoc;
err = dlgdoc.Parse(xml$); 
### get the dialog element from the parsed XML document
class XMLNODE dlgnode;
dlgnode = dlgdoc.GetElementByID("dlgxtiff");
### set the dialog XML element as the source for the GUI_DLG class instance
dlgwin.SetXMLNode(dlgnode);
dlgwin.CreateModeless();
### get handles for dialog controls
outDirBtn = dlgwin.GetCtrlByID("outDirBtn");
indirtext = dlgwin.GetCtrlByID("indirtext");
outdirtext = dlgwin.GetCtrlByID("outdirtext");
numTiff = dlgwin.GetCtrlByID("numTiff");
numShape = dlgwin.GetCtrlByID("numShape");
numProc = dlgwin.GetCtrlByID("numProc");
TIFFlistbox = dlgwin.GetCtrlByID("TIFFlistbox");
Shapelistbox = dlgwin.GetCtrlByID("Shapelistbox");
Proclistbox = dlgwin.GetCtrlByID("Proclistbox");
statustext = dlgwin.GetCtrlByID("statustext");
runBtn = dlgwin.GetCtrlByID("runBtn");
cancelBtn = dlgwin.GetCtrlByID("cancelBtn");
closeBtn = dlgwin.GetCtrlByID("closeBtn");
dlgwin.Open();
WaitForExit();