ExportMultiSectColladaKMZtrack.sml

  Download

More scripts: Advanced

Syntax Highlighing:

comments, key words, predefined symbols, class members & methods, functions & classes
            
# ExportMultiSectColladaKMZtrack.sml
# write multiple geologic cross-section images to Collada files packaged in a
# KMZ file for viewing in Google Earth. The contained KML file animates extruding the
# sections from their true below-ground positions (where they are not visible in Google
# Earth). The elevation range for the animation is automatically set to raise 
# all sections above the terrain while keeping their correct relative elevations.
# Cross-sections must have 3D Piecewise Affine (manifold) georeference model
# with 3D control points.  Each section's coordinate reference system can be any
# projected CRS or a nonprojected geographic (latitude/longitude) CRS.
# Animation uses <gx:Track> KML element supported in Google Earth 5.2 and later.
# Provides smoother and simpler animation than the multiple-placemark method.
# This script version also sets initial Google Earth camera parameters to 
# establish a view direction perpendicular to the section trend.  
# Randy Smith, MicroImages, Inc.
# 29 July 2010
class RVC_OBJITEM tempObjItemList[];	# hash of ObjItems for the selected cross-section rasters 
class RVC_OBJITEM xsecObjItemList[];	# hash of ObjItems for valid cross-section rasters
class RVC_RASTER XSEC;
class RVC_GEOREFERENCE xsecgeoref, tempgeoref;
numeric i, j, k; # counters
numeric err; # for error checking
numeric numSections;
numeric minZm, maxZm;	# min and max elevations among all section control points (always stored in meters)
minZm = 30000; # starting value higher than any terrestrial elevation in meters;
maxZm = -1000;	# starting value lower than any terrestrial elevation in meters
numeric numPts;
# number of steps in extrusion animation 
# (= number of cross-section placemarks to write to KML file) 
numeric numExtrudeSteps = 10;
clear();
##############   FUNCTIONS AND PROCEDURES   #########################
# error checking procedure
proc ReportError(numeric linenum, numeric err) {
	printf("FAILED -line: %d, error: %d\n", linenum - 1, err);
	PopupError(err); 
	}
# function to determine if three points are in counter-clockwise order
# using image coordinates (origin in upper left, positive downward);
# returns value < 0 if counter-clockwise, value > 0 if clockwise.
func ccwImageCoords(class POINT2D p1, class POINT2d p2, class POINT2d p3)
	{
	return (p2.x - p1.x)*(p3.y - p1.y) - (p2.y - p1.y)*(p3.x - p1.x);
	}
###############################################################################
##############################   MAIN PROGRAM   ###############################
############### Get input cross-sections
class STRING prompt$ = "Choose cross-section rasters to export to COLLADA:";
DlgGetObjects(prompt$, "Raster", tempObjItemList, "ExistingOnly", 1);
numSections = tempObjItemList.GetNumItems();
printf("%d rasters selected.\n", numSections );
array numeric statusArray[numSections]; # array to record whether section rasters are valid
# check for manifold georeference and delete objitems of invalid objects from xsecObjItemList
for i = 1 to numSections
	{
	XSEC.Open(tempObjItemList[i], "Read");
	# get georeference subobject and check that it is present and manifold type with valid CRS
	if (XSEC.GetDefaultGeoref(tempgeoref) )	# true if object has a georeference subobject
		{
		if (tempgeoref.GetGeorefType() == "Manifold" )
			{
			printf("Section object %d has manifold georeference.\n", i);
			if (tempgeoref.GetCoordRefSys().IsLocal() )
				{
				print("Section raster %d has arbitrary local georeference. Removing from input list.\n", i);
				}
			else
				{
				statusArray[i] = 1;
				xsecgeoref.OpenLastUsed(XSEC);
				numPts = xsecgeoref.GetNumCtrlPoints();
				class CTRLPOINT TctrlPtList[numPts];
				xsecgeoref.GetPointList(TctrlPtList);
				for j = 1 to numPts
					{ 
					if (TctrlPtList[j].dest.z < minZm) then minZm = TctrlPtList[j].dest.z;
					if (TctrlPtList[j].dest.z > maxZm) then maxZm = TctrlPtList[j].dest.z;
					}
				xsecgeoref.Close();
				}
			}
		else
			{
			printf("Section object %d does not have manifold georeference. Removing from input list.\n", i);
			statusArray[i] = 0;
			}
		}
	else
		{
		printf("Section object %d is not georeferenced. Removing from input list.\n", i);
		statusArray[i] = 0;
		}
	XSEC.Close();
	}
# loop through statusArray to copy objItems for sections with valid manifold georeference
# to the final objItemList
numeric validcount;		# count of valid sections to be moved to the new objItemList
for i = 1 to numSections
	{
	if (statusArray[i] == 1 ) 
		{
		++validcount;
		xsecObjItemList[validcount] = tempObjItemList[i];
		}
	}
if (xsecObjItemList.GetNumItems() == 0 )
	{
	print("No selected cross-section raster has valid manifold georeference. Exiting now.");
	Exit();
	}
printf("Number of valid cross-sections to process = %d\n", xsecObjItemList.GetNumItems() );
printf("Minimum elevation in meters among all section control points = %.1f\n", minZm);
printf("Maximum elevation in meters among all section control points = %.1f\n", maxZm); 
# prompt user for the name of the output KMZ file
class STRING kmzfilename$;
prompt$ = "Choose and output KMZ file to package the cross-section data \n (spaces will be replaced with underscores):"; 
kmzfilename$ = GetOutputFileName(_context.ScriptDir, prompt$, ".kmz");
for i = 0 to kmzfilename$.length - 1
	{
	if (kmzfilename$.charAt(i) == " ") then
		kmzfilename$.replace(i, i, "_");
	}
class FILEPATH kmzFilepath(kmzfilename$);
class FILEPATH dirFilepath(kmzfilename$);	# filepath to parent directory of KMZ file
dirFilepath.StripToFolder();	# path to directory containing the KMZ
class STRING dir$;
dir$ = dirFilepath;
class STRING dataName$;
dataName$ = kmzFilepath.GetNameOnly();
# set up a WGS84 / Geographic CRS for lat/lon coordinates to position sections
class SR_COORDREFSYS geoCRS;
geoCRS.Assign("Geographic2D_WGS84_Deg");
printf("Geographic CRS = %s\n", geoCRS.Name);
####################################################################
# string with COLLADA file template for a cross-section 
class STRING dae$ = '<?xml version="1.0" encoding="utf-8"?>
   <COLLADA xmlns="http://www.collada.org/2005/11/COLLADASchema" version="1.4.1">
   <asset>
      <contributor>
         <authoring_tool>TNTmips SML Script</authoring_tool>
      </contributor>
      <created></created>
      <modified></modified>
      <unit name="" meter="1.0"/>
      <up_axis>Z_UP</up_axis>
   </asset>
   <library_images>
      <image id="Xsect-image" name="Xsect-image">
         <init_from></init_from>
      </image>
   </library_images>
   <library_materials>
      <material id="XsectID" name="Xsect">
         <instance_effect url="#Xsect-effect"/>
      </material>
   </library_materials>
   <library_effects>
      <effect id="Xsect-effect" name="Xsect-effect">
         <profile_COMMON>
            <newparam sid="Xsect-image-surface">
               <surface type="2D">
                  <init_from>Xsect-image</init_from>
               </surface>
            </newparam>
            <newparam sid="Xsect-image-sampler">
               <sampler2D>
                  <source>Xsect-image-surface</source>
               </sampler2D>
            </newparam>
            <technique sid="COMMON">
               <phong>
                  <emission>
                     <color>0.000000 0.000000 0.000000 1</color>
                  </emission>
                  <ambient>
                     <color>0.000000 0.000000 0.000000 1</color>
                  </ambient>
                  <diffuse>
                     <texture texture="Xsect-image-sampler" texcoord="UVSET0"/>
                  </diffuse>
                  <specular>
                     <color>0.330000 0.330000 0.330000 1</color>
                  </specular>
                  <shininess>
                     <float>20.000000</float>
                  </shininess>
                  <reflectivity>
                     <float>0.100000</float>
                  </reflectivity>
                  <transparent>
                     <color>1 1 1 1</color>
                  </transparent>
                  <transparency>
                     <float>0.000000</float>
                  </transparency>
               </phong>
            </technique>
            <extra>
               <technique profile="GOOGLEEARTH">
                  <double_sided>1</double_sided>
               </technique>
            </extra>
         </profile_COMMON>
      </effect>
   </library_effects>
   <library_geometries>
      <geometry id="mesh1-geometry" name="mesh1-geometry">
         <mesh>
            <source id="mesh1-geometry-position">
               <float_array id="mesh1-geometry-position-array" count="">
               		 </float_array>
               <technique_common>
                  <accessor source="#mesh1-geometry-position-array" count="" stride="3">
                     <param name="X" type="float"/>
                     <param name="Y" type="float"/>
                     <param name="Z" type="float"/>
                  </accessor>
               </technique_common>
            </source>
            <source id="mesh1-geometry-normal">
               <float_array id="mesh1-geometry-normal-array" count="3">-0.000000 -1.000000 -0.000000 </float_array>
               <technique_common>
                  <accessor source="#mesh1-geometry-normal-array" count="1" stride="3">
                     <param name="X" type="float"/>
                     <param name="Y" type="float"/>
                     <param name="Z" type="float"/>
                  </accessor>
               </technique_common>
            </source>
            <source id="mesh1-geometry-uv">
               <float_array id="mesh1-geometry-uv-array" count="">
               		 </float_array>
               <technique_common>
                  <accessor source="#mesh1-geometry-uv-array" count="" stride="2">
                     <param name="S" type="float"/>
                     <param name="T" type="float"/>
                  </accessor>
               </technique_common>
            </source>
            <vertices id="mesh1-geometry-vertex">
               <input semantic="POSITION" source="#mesh1-geometry-position"/>
            </vertices>
            <triangles material="Xsect" count="">
               <input semantic="VERTEX" source="#mesh1-geometry-vertex" offset="0"/>
               <input semantic="NORMAL" source="#mesh1-geometry-normal" offset="1"/>
               <input semantic="TEXCOORD" source="#mesh1-geometry-uv" offset="2" set="0"/>
               <p> </p>
            </triangles>
         </mesh>
      </geometry>
   </library_geometries>
   <library_visual_scenes>
      <visual_scene id="CrossSectionScene" name="CrossSectionScene">
         <node id="Model" name="Model">
            <node id="mesh1" name="mesh1">
               <instance_geometry url="#mesh1-geometry">
                  <bind_material>
                     <technique_common>
                        <instance_material symbol="Xsect" target="#XsectID">
                           <bind_vertex_input semantic="UVSET0" input_semantic="TEXCOORD" input_set="0"/>
                        </instance_material>
                     </technique_common>
                  </bind_material>
               </instance_geometry>
            </node>
         </node>
      </visual_scene>
   </library_visual_scenes>
   <scene>
      <instance_visual_scene url="#CrossSectionScene"/>
   </scene>
</COLLADA>
';
# variables for modifying the COLLADA template for each section
class DATETIME dt;
class STRING dt$, format$;
class XMLDOC daedoc;
class XMLNODE colladaRootNode;
class XMLNODE node, node1, node2, node3;
numeric arraycount;
###############################################################
# string with template for KML file
class STRING kmldoc$;
class STRING kmlfile$ = dir$ + "/" + "doc.kml";
class FILEPATH kmlFilepath(kmlfile$);
kmldoc$ = '<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">
<Folder>
	<name></name>
	<open>1</open>
</Folder>
</kml>';
# variables for modifying the KML file template
class XMLDOC kmldoc;
class XMLNODE kmlRootNode;
class XMLNODE folder, document, placemark, track, model;
numeric startDate;
numeric date;
numeric altitudeIncrement;
##############################
# parse KML file template to reference and position the COLLADA models
err = kmldoc.Parse(kmldoc$);
if (err < 0) {	PopupError(err); Exit(); }
kmlRootNode = kmldoc.GetRootNode();	# get root node 
folder = kmlRootNode.FindChild("Folder");		# get folder node that will hold a Document for each cross-section
node = folder.FindChild("name");
node.SetText(dataName$);		# set name of data collection
##################################################################
################ Create KMZ file (ZIP file with .kmz file extension
class ZIPFILE kmz;
###########################################################################
# variables for processing the cross-section rasters and their control points
class RVC_DESCRIPTOR xsecDescriptor;
class STRING xsecName$;
class SR_COORDREFSYS xsecCRS;
class STRING unit$;
numeric geogCRSflag;
class SR_COORDAXIS axisX;
numeric unitConvToMeters;		# scale factor to convert CRS map unit to meters
numeric metersConvToUnit;		# scale factor to convert meters to CRS map unit
class BOUNDARYLIST bndryList;		# container for boundary point arrays
array numeric bndryArray[1];	# array indices for the current boundary
array numeric xyz_array[1];	# holds 3D map coordinates for all points: X Y Z X Y Z etc.
array numeric uv_array[1];	# holds normalized texture coordinates: S T S T etc.
# arrays to hold map coordinates (destination) used to determine coordinate
# ranges and compute local 3D position coordinates for the COLLADA model 
array numeric xarray[1];
array numeric yarray[1];
array numeric zarray[1];
numeric minX, minY, minZ;		# minimum values for each of x, y, and z coordinates for the section control points in CRS units
numeric maxX, maxY;				# maximum value for X and Y from control points
numeric maxZ;						# maximum z value from control points in CRS units
class POINT2D globalMinPt, globalMaxPt;	# global minimum and maximum extents points for group of sections
globalMinPt.x = 360;	globalMinPt.y = 90;		# starting values greater than any anticipated values
globalMaxPt.x = -360; globalMaxPt.y = -90;
numeric sectionLength, sectionAzimuth, sectionElev;
array numeric sectionAzimuths[numSections];		# array to hold azimuths for all sections
# variables for processing control point coordinates
class STRING xyzString, uvString;		# string arrays of coordinates to write to COLLADA
numeric col, lin, s, t;
class POINT3D destPt;
class POINT3D originPt;		# origin point for COLLADA model to place it correctly in Google Earth
numeric angleToNorth;
altitudeIncrement = (maxZm - minZm) / numExtrudeSteps;####################
printf("Altitude increment for %d extrusion steps = %.1f\n", numExtrudeSteps, altitudeIncrement);
# variables for creating and processing triangulation from the control points
class TRIANGULATOR triangulator; # Triangulator class to create a triangulation in image space for the control points 
class TRIANGLENODES triNodes;
class CTRLEDGE hardedge;
numeric nodeNum, nodeIndex, diffNode;
numeric numTriangles;
class POINT2D pt, nodePt;
array numeric triNodeIndexArray[3];
array numeric tempArray[3];
numeric testCCW;
class STRING pString;	# string array of indices into coordinate arrays to represent their triangulations in COLLADA file
#########################################################################
#########################################################################
#### loop through objItemList of valid sections to process them
startDate = 2000; 
for k = 1 to validcount
	{
	printf("\nProcessing section %d of %d\n", k, validcount);
	XSEC.Open(xsecObjItemList[k], "Read");
	xsecgeoref.OpenLastUsed(XSEC);
	xsecDescriptor = XSEC.GetDescriptor();
	xsecName$ = xsecDescriptor.GetShortName();
	printf("Name of cross-section raster %d is: %s\n", k, xsecName$);
	# get Coordinate Reference System from section
	xsecCRS = xsecgeoref.GetCoordRefSys();
	printf("Section %d CRS = %s\n", k, xsecCRS.Name);
	if (xsecCRS.IsProjected() )	# projected CRS with cartesian map coordinates
		{
		printf("Section %d has projected CRS.\n", k);
		axisX = xsecCRS.CoordSys.GetAxis();
		unit$ = axisX.Unit.GetName();
		# get scale factors for converting between units in the projected CRS and meters.
		# Scale factor from unit to meters is required as a parameter in the COLLADA file.
		# Elevations in manifold georeference are always stored as meters, so need to convert
		# elevation values to units in projected CRS as well.
		unitConvToMeters = GetUnitConvDist(unit$, "meter"); 
		printf("Map unit = %s, conversion factor to meters = %.4f\n", unit$, unitConvToMeters);
		metersConvToUnit = GetUnitConvDist("meter", unit$);
		}
	else 
		{					# nonprojected = geographic (latitude longitude)
		printf("Section %d has geographic CRS.\n", k);
		geogCRSflag = 1;	# set flag indicating CRS is geographic 
		unit$ = "meters";
		metersConvToUnit = 1;
		unitConvToMeters = 1;
		# find center point of section in lat/lon coordinates to determine a UTM
		# zone to use for a temporary projected coordinate system for position coordinates
		# in the COLLADA file
		class TRANS2D_MAPGEN transObjToMap;
		xsecgeoref.GetTransParm(transObjToMap, 0, xsecgeoref.GetCalibModel() );
		class RECT3D extents;
		XSEC.GetExtents(extents);
		class POINT2D center;
		center = extents.center;
		printf("center object coordinates: x = %.1f, y = %.1f\n", center.x, center.y);
		center = transObjToMap.ConvertPoint2DFwd(center);
		printf("center map coordinates: x = %.6f, y = %.6f\n", center.x, center.y);
		# variables for setting up the UTM CRS
		class SR_COORDREFSYS utmCRS;
		class SR_COORDSYS utmCoordSys;		# coordinate system
		class SR_COORDOPDEF projectDef;		# projection parameters and method
		class SR_COORDOPMETHOD method;
		class SR_DATUM datum;
		datum = xsecCRS.Datum;
		utmCoordSys.Assign("Projected2D_EN_m");	# assign projected coordinate system
		method.Assign("1909");					# assign using numeric ID for Transverse Mercator projection
		projectDef.Method = method;		# assign method to projection
		# assign projection parameters that are the same for all UTM zones
		projectDef.SetParmValueNum("LatitudeOfNaturalOrigin", 0);
		projectDef.SetParmValueNum("ScaleFactorAtNaturalOrigin", 0.9996);
		projectDef.SetParmValueNum("FalseEasting", 500000);
		# assign false northing based on whether north or south of equator
		if ( center.y < 0) then
			projectDef.SetParmValueNum("FalseNorthing", 10000);
		else
			projectDef.SetParmValueNum("FalseNorthing", 0);
		# assign longitude of natural origin (central meridian of UTM zone) based on longitude
		numeric zoneNum, centMerid;		# UTM zone and longitude of central meridian
		zoneNum = ceil( (180 + center.x) / 6 );
		centMerid = (zoneNum * 6) - 183;
		projectDef.SetParmValueNum("LongitudeOfNaturalOrigin", centMerid);
		# create defined UTM coordinate reference system
		utmCRS.Create(utmCoordSys, datum, projectDef);
		# set up coordinate transformation from section's geographic CRS to the defined UTM CRS
		class TRANS2D_MAPGEN transGeogUTM;
		transGeogUTM.InputCoordRefSys = xsecCRS;
		transGeogUTM.OutputCoordRefSys = utmCRS;
		}
	# set up coordinate transformation from section CRS to geoCRS
	class TRANS2D_MAPGEN transToGeo;
	transToGeo.InputCoordRefSys = xsecCRS;
	transToGeo.OutputCoordRefSys = geoCRS;
	# report dimensions of the cross-section raster
	printf("Cross section %d dimensions: %d lines, %d columns\n", k, XSEC.$Info.NumLins, XSEC.$Info.NumCols);
	###############################################################
	#################  make a PNG file from the cross-section raster 
	# PIPELINE SOURCE
	class IMAGE_PIPELINE_SOURCE_RVC sourceXS(xsecObjItemList[k]);
	err = sourceXS.Initialize();
	if ( err < 0) ReportError(_context.CurrentLineNum, err);
	# PIPELINE FILTER NULL TO ALPHA
	# converts null mask to alpha channel for output PNG
	class IMAGE_PIPELINE_FILTER_NULLTOALPHA alpha(sourceXS);
	err = alpha.Initialize();
	if ( err < 0) ReportError(_context.CurrentLineNum, err);
	# PIPELINE TARGET
	class STRING png$ = dir$ + "/" + xsecName$ + ".png";
	class FILEPATH pngFilepath(png$);
	printf("PNG filepath = %s\n", png$);
	class IMAGE_PIPELINE_TARGET_PNG targetPNG(alpha, pngFilepath, "none");
	err = targetPNG.Initialize();
	if ( err < 0) ReportError(_context.CurrentLineNum, err); 
	# Execute the pipeline
	err = targetPNG.Process();
	if ( err < 0) ReportError(_context.CurrentLineNum, err);
	#####################################################################
	############### process manifold georeference control points
	# get the list of control points as an array of class CTRLPOINT
	printf("Getting list of georeference control points for section %d.\n", k);
	numeric numPts = xsecgeoref.GetNumCtrlPoints(); 
	class CTRLPOINT ctrlPtList[numPts];
	xsecgeoref.GetPointList(ctrlPtList);
	numeric numBndry = xsecgeoref.GetNumBoundaries();
	numeric numEdges = xsecgeoref.GetNumCtrlEdges();
	printf("Number of control points = %d\n", numPts);
	printf("Number of control edges = %d\n",  numEdges);
	printf("Number of triangulation boundaries = %d\n\n", numBndry );
	###### get boundary of triangulated control points if used (may be more than one boundary)
	print("Get control point boundary if present and list indices:");
	numeric numBndryPts;		# number of points in current boundary
	if (numBndry > 0)
		{
		if (bndryList.GetNumBoundaries() > 0) then
			bndryList.ClearBoundaries();
		for j = 1 to numBndry
			{
			numBndryPts = xsecgeoref.ReadBoundary(j-1, bndryArray);	# indexing of boundaries is 0-based
			printf("Number of points in boundary %d = %d\n", j, numBndryPts);
			for i = 1 to numBndryPts {
				printf("Boundary index %d = %d\n", i, bndryArray[i] + 1);
				}
			printf("\nAdd control point boundary to boundary list.\n");
			bndryList.AddBoundary(bndryArray, numBndryPts);
			ResizeArrayClear(bndryArray, 1);
			}
		}
	printf("Check boundary list for section %d:\n", k);
	printf("Number of boundaries in boundary list = %d\n", bndryList.GetNumBoundaries() );
	# resize and clear array holding all map coordinates (from destination points) and
	# texture coordinates (source points).
	ResizeArrayClear(xyz_array, numPts * 3);
	ResizeArrayClear(uv_array, numPts * 2);
	# resize and clear arrays for x and y coordinates used to determine coordinate ranges and
	# compute local 3D position coordinates for the COLLADA model
	ResizeArrayClear(xarray, numPts);
	ResizeArrayClear(yarray, numPts);
	ResizeArrayClear(zarray, numPts);
	minX = 100000000;	minY = 100000000;		# starting values greater than any anticipated real values
	maxX = 0;  maxY = 0;
	# loop through points to read the destination (map) coordinates to arrays and determine
	# minimum values for x and y to use to create offset coordinates for each axis of
	# the COLLADA model; 
	for i = 1 to numPts
		{
		printf("Control Pt %d:\n", i);
		# get map point from current control point
		destPt = ctrlPtList[i].dest;
		if (geogCRSflag)  	# if input CRS is geographic, convert to our virtual UTM coordinates
			{
			destPt = transGeogUTM.ConvertPoint3DFwd(destPt);
			destPt.z = ctrlPtList[i].dest.z;
			}
		xarray[i] = destPt.x;
		yarray[i] = destPt.y;
		zarray[i] = destPt.z * metersConvToUnit;  # apply scaling from meters to coordinate system unit
		printf("x = %.6f, y = %.6f\n", xarray[i], yarray[i]);
		if (xarray[i] < minX) then minX = xarray[i];
		if (yarray[i] < minY) then minY = yarray[i];
		if (xarray[i] > maxX) then maxX = xarray[i];
		if (yarray[i] > maxY) then maxY = yarray[i];
		}
	printf("\nMinimum coordinate values for offset for this section: x = %.4f, y = %.4f\n", minX, minY);
	printf("Maximum coordinate values for this section: x = %.4f, y = %.4f\n", maxX, maxY);
	# convert min and max coordinates to geographic and update global min and max extents of the set of sections
	class POINT2D ptmin, ptmax;
	ptmin.x = minX; ptmin.y = minY;
	ptmax.x = maxX; ptmax.y = maxY;
	if (geogCRSflag) {	# if input geographic CRS, convert UTM origin point back to original CRS
		ptmin = transGeogUTM.ConvertPoint2DInv(ptmin);
		ptmax = transGeogUTM.ConvertPoint2DInv(ptmax);
		}
	else {
		ptmin = transToGeo.ConvertPoint2DFwd(ptmin);
		ptmax = transToGeo.ConvertPoint2DFwd(ptmax); 
		}
	if (ptmin.x < globalMinPt.x) then globalMinPt.x = ptmin.x;
	if (ptmin.y < globalMinPt.y) then globalMinPt.y = ptmin.y;
	if (ptmax.x > globalMaxPt.x) then globalMaxPt.x = ptmax.x;
	if (ptmax.y > globalMaxPt.y) then globalMaxPt.y = ptmax.y;
	printf("Current global extents: min lon = %.8f deg, min lat = %.8f deg, max lon = %.8f, max lat = %.8f\n", globalMinPt.x, globalMinPt.y, globalMaxPt.x, globalMaxPt.y);
	# find azimuth of section
	sectionAzimuth = 0;
	Displacement3Dd(minX, minY, 0, maxX, maxY, 0, sectionLength, sectionAzimuth, sectionElev);
	sectionAzimuths[k] = sectionAzimuth;
	printf("Section azimuth = %.1f\n", sectionAzimuths[k]);	 
	# set min x, min y as origin point for the model; find its coordinates in WGS84 Geographic CRS and
	# find angle to north at this location in original CRS to set position and orientation of model in Google Earth
	originPt.x = minX;  originPt.y = minY;		# coordinates in CRS
	# if input geographic CRS, convert UTM origin point back to original CRS
	if (geogCRSflag) then originPt = transGeogUTM.ConvertPoint3DInv(originPt);
	angleToNorth = xsecCRS.ComputeAngleToNorth(originPt);
	originPt = transToGeo.ConvertPoint3DFwd(originPt);		# reproject to WGS 84 / Geographic 
	# set Z value for origin point in meters to provide correct vertical position for model in Google Earth
	originPt.z = minZm;	# minZ in meters for all sections
	printf("origin point: angle to north = %.1f, lon = %.6f, lat = %.6f z = %.1f m\n\n", angleToNorth, originPt.x, originPt.y, originPt.z);
	# loop through points again to write offset XYZ values to xyzString, read image coordinates and
	# write to uvString, and add original points to triangulator to determine mapping of vertices to
	# triangles for position mesh in COLLADA file
	xyzString.Clear();
	uvString.Clear();
	pString.Clear();
	triangulator.ClearAll();
	minZ = minZm * metersConvToUnit;		# origin elevation in CRS units as offset for z
	print("Create strings with vertex coordinates for position and texture meshes in COLLADA file:"); 
	for i = 1 to numPts
		{
		# add offset 3D map coordinates to the xyzString
		xyzString += sprintf("%.6f %.6f %.6f ", xarray[i] - minX, yarray[i] - minY, zarray[i] - minZ); 
		printf("xyzString = %s\n", xyzString);
		# read image coordinates
		col = ctrlPtList[i].source.x;		lin = ctrlPtList[i].source.y;
		# convert to texture coordinates (normalized to 0 to 1 range, origin in lower left) 
		s = col / XSEC.$Info.NumCols;
		t = 1 - (lin / XSEC.$Info.NumLins);		# inverted due to lower left origin
		# add to uvString
		printf("Image coordinates: col = %.6f, lin = %.6f\n", col, lin);
		printf("Texture coordinates: s = %.6f, t = %.6f\n", s, t);
		uvString += sprintf("%.6f %.6f ", s, t);
		printf("uvString = %s\n", uvString);
		triangulator.AddPoint(ctrlPtList[i].source);
		}
	# add hard edges from georeference (if any) to the triangulator
	printf("\nAdd hard edges to triangulator:\n");
	if (numEdges > 0)
		{
		for i = 1 to xsecgeoref.GetNumCtrlEdges()
			{ 
			xsecgeoref.ReadCtrlEdge(i-1, hardedge);	# internal indexing of hard edges is 0-based, so use i-1
			# internal indexing of hard edge points is 0-based, so add 1 to correspond to 1-based control point numbering
			printf("Hardedge %d: start = %d, end = %d\n", i-1, hardedge.start + 1, hardedge.end + 1);
			triangulator.AddHardEdgeAsIndices(hardedge.start + 1, hardedge.end + 1);
			}
		}
	# triangulate the control points
	err = triangulator.Triangulate(0, "ShortestEdge");
	printf("\nTriangulating, triangulation returned %d\n", err); 
	printf("Number of nodes = %d\n", triangulator.GetNumNodes() );
	printf("Number of edges = %d\n", triangulator.GetNumEdges() );
	printf("Number of triangles = %d\n", triangulator.GetNumTriangles() );
	## add control point boundary read from the georeference object to the triangulator.
	## This must be done AFTER the points have been triangulated; points and triangles outside
	## the boundary are then discarded.
	if (numBndry > 0) 
		{
		printf("\nAdd control point boundary to triangulator.\n");
		triangulator.SetBoundaries(bndryList);
		# check number of elements again
		printf("Number of nodes = %d\n", triangulator.GetNumNodes() );
		printf("Number of edges = %d\n", triangulator.GetNumEdges() );
		printf("Number of triangles = %d\n\n", triangulator.GetNumTriangles() );
		}
	# The triangulator creates a new set of Node Numbers but provides a GetIndex() method
	# to find the node index (1-based) which is set by the order in which the points were added
	# to the triangulator.  Node indices therefore map to the georeference control point number.  
	# Hard edges added to the set result in additional nodes (each at the
	# same location as one of the original points). 
	print("Loop through nodes to report indices and coordinates.");
	for i = 1 to triangulator.GetNumNodes()
		{
		triangulator.GetPoint(i, pt);
		printf("Node %d, index = %d, x = %.6f, y = %.6f\n", i, triangulator.GetIndex(i), pt.x, pt.y);
		}
	numTriangles = triangulator.GetNumTriangles();
	# loop through the triangles to get the nodes for each and find the corresponding georeference control points
	for i = 1 to numTriangles
		{
		err = triangulator.GetTriangleNodes(i, triNodes);
		printf("\nTriangle %d, GetTriangleNodes returned %d\n", i, err);
		for j = 1 to 3 	# loop through nodes of triangle to get the control point number (node index) for each node
			{
			nodeNum = triNodes.GetNode(j);
			nodeIndex = triangulator.GetIndex(nodeNum);
			if (nodeIndex > numPts)	# check if node was added with a hard edge and if so find the
													# corresponding node from the original control point set
				{
				triangulator.GetPoint(nodeNum, nodePt);
				nodeNum = triangulator.FindDifferentNode(nodePt, nodeNum);
				nodeIndex = triangulator.GetIndex(nodeNum);
				}
			printf("Node %d = %d, Point Index = %d\n", j, nodeNum, nodeIndex);
			tempArray[j] = nodeIndex;	# add node index to the array of triangle node indices
			}
		printf("Temp array indices: %d, %d, %d\n", tempArray[1], tempArray[2], tempArray[3]);
		testCCW = ccwImageCoords(ctrlPtList[ tempArray[1] ].source,	ctrlPtList[ tempArray[2] ].source, ctrlPtList[ tempArray[3] ].source);
		printf("counter-clockwise test returned %.4f\n", testCCW);
		triNodeIndexArray[1] = tempArray[1];
		if (testCCW > 0) 	# if nodes not in counterclockwise order, reorder nodes 2 and 3
			{
			triNodeIndexArray[2] = tempArray[3];
			triNodeIndexArray[3] = tempArray[2];
			}
		else
			{
			triNodeIndexArray[2] = tempArray[2];
			triNodeIndexArray[3] = tempArray[3];
			}
		printf("Triangle node indices in counter-clockwise order: %d, %d, %d\n", triNodeIndexArray[1], triNodeIndexArray[2], triNodeIndexArray[3]);
		# loop through the nodes for the current triangle to write index positions to the "primitive" string (pString)
		# for the COLLADA file.  For each triangle need three values, in order: index for vertex, index for normals,
		# index for texture. All indices are 0-based; index for vertex and texture are equal to each other and equal
		# node index - 1.
		for j = 1 to 3
			{
			pString += sprintf("%d %d %d ", triNodeIndexArray[j] - 1, 0, triNodeIndexArray[j] - 1);
			}
		printf("pString = %s\n\n", pString);
		}
	# close the georeference object and its parent raster object
	xsecgeoref.Close();
	XSEC.Close();
	# get current date and time for inclusion in the COLLADA file
	dt.SetCurrent();
	format$ = "%Y-%m-%dT%H:%M:%SZ";
	dt$ = dt.GetString(format$);
	printf("formatted date and time: %s\n", dt$);
	#######################################################################
	# parse the onesection.dae COLLADA file template and get the root node
	print("Making COLLADA file.");
	err = daedoc.Parse(dae$);
	if (err < 0) {
		PopupError(err);
		break; 
		}
	colladaRootNode = daedoc.GetRootNode();
	# set created and modified dates, units, and scale factor to meters
	node = colladaRootNode.FindChild("asset");
	node1 = node.FindChild("created");
	node1.SetText(dt$);
	node1 = node.FindChild("modified");
	node1.SetText(dt$);
	node1 = node.FindChild("unit");
	node1.SetAttribute("name", unit$);
	node1.SetAttribute("meter", sprintf("%.4f", unitConvToMeters) );
	### set the output PNG file as the library image
	node = daedoc.GetElementByID("Xsect-image");
	node1 = node.FindChild("init_from");
	node1.SetText(pngFilepath.GetName() );
	### set the position mesh count and array
	arraycount = numPts * 3;
	node = daedoc.GetElementByID("mesh1-geometry-position-array");
	node.SetAttribute("count", NumToStr(arraycount) );
	node.SetText(xyzString);
	### set the count for the position array accessor
	node1 = node.GetParent();
	node2 = node1.FindChild("technique_common");
	node3 = node2.FindChild("accessor");
	node3.SetAttribute("count", NumToStr(numPts) );
	### set the texture mesh count and array
	arraycount = numPts * 2;
	node = daedoc.GetElementByID("mesh1-geometry-uv-array");
	node.SetAttribute("count", NumToStr(arraycount) );
	node.SetText(uvString);
	### set the count for the uv array accessor
	node1 = node.GetParent();
	node2 = node1.FindChild("technique_common");
	node3 = node2.FindChild("accessor");
	node3.SetAttribute("count", NumToStr(numPts) );
	### set the material and primitive for the triangles element
	node = daedoc.GetElementByID("mesh1-geometry");
	node1 = node.FindChild("mesh");
	node2 = node1.FindChild("triangles");
	node2.SetAttribute("count", NumToStr(numTriangles) );
	node3 = node2.FindChild("p");
	node3.SetText(pString);
	############## write out the COLLADA file
	class STRING daefile$ = dir$ + "/" + xsecName$ + ".dae";
	printf("\nDAE filepath = %s\n", daefile$);
	daedoc.Write(daefile$);
	class FILEPATH daeFilepath(daefile$);
	##################################################################
	#############
	############# add the section PNG file and COLLADA file to the KMZ
	kmz.AddFile(kmzFilepath, pngFilepath, "files/" + xsecName$ + ".png");
	kmz.AddFile(kmzFilepath, daeFilepath, "files/" + xsecName$ + ".dae");
	# delete the source files
	DeleteFile(png$);
	DeleteFile(daefile$);
	##################################################################
	############
	############  add info to the KML structure for this section
	# create <Document> element for this section in <Folder>
	document = folder.NewChild("Document");
	node = document.NewChild("name");
	node.SetText(xsecName$);
	# create <Placemark> for the cross-section with name
	placemark = folder.NewChild("Placemark");
	placemark.NewChild("name", xsecName$);
	# create <gx:Track> element for <Placemark> and set <altitudeMode>
	track = placemark.NewChild("gx:Track");
	track.NewChild("altitudeMode", "absolute");
	# create <Model> element in <gx:Track>
	model = track.NewChild("Model");
	model.SetAttribute("id", sprintf("model%d", k) ); 
	# create <Orientation> element for <Model>
	node = model.NewChild("Orientation");
	node.NewChild("heading", NumToStr(angleToNorth) );
	node.NewChild("tilt", "0");
	node.NewChild("roll", "0");
	# create <Scale> element for <Model>
	node = model.NewChild("Scale");
	node.NewChild("x", "1.0");
	node.NewChild("y", "1.0");
	node.NewChild("z", "1.0");
	# create <Link> element for <Model>
	node = model.NewChild("Link");
	node.NewChild("href", "files/" + xsecName$ + ".dae");
	# create <ResourceMap> element for <Model>
	node = model.NewChild("ResourceMap");
	node1 = node.NewChild("Alias");
	node1.NewChild("targetHref", "files/" + xsecName$ +".png");
	node1.NewChild("sourceHref", "files/" + xsecName$ +".png");
	# create <when>, <gx:coord>, and <gx:angles> elements for each extrusion step
	for i = 1 to numExtrudeSteps + 1
		{
		date = startDate + ( (i-1) * 100);	# increment date by 100 years for each step
		# create <when> element for step with this date
		node = track.NewChild("when", NumToStr(date) );
		# create <gx:coord> element for step, with X, Y, and Z values, separated by spaces, as element text
		node = track.NewChild("gx:coord", sprintf("%.6f %.6f %.1f", originPt.x, originPt.y, originPt.z + (i-1) * altitudeIncrement ) );
	
		# create <gx:angles> element with heading, tilt, and roll angles (set all to 0)
		# even though we set angles to 0, without this element for each step the section is not correctly oriented
		node = track.NewChild("gx:angles", "0 0 0");
		printf("Extrude step %d: altitude = %.1f m\n", i - 1, originPt.z + (i-1) * altitudeIncrement);
		}
	}	# end:  for k = 1 to validcount
##########  Get map extents of sections (already in geographic coordinates) and use to set viewpoint parameters 
class RECT sectionsExtents;
sectionsExtents.pt1 = globalMinPt;
sectionsExtents.pt2 = globalMaxPt;
# find diagonal dimension of the extents in meters
numeric ew, ns, range;
ew = sectionsExtents.GetWidth() * 111319.5;
ns = sectionsExtents.GetHeight() * 111319.5;
range = sqrt( sqr(ew) + sqr(ns) );
printf("Long dimension of sections extents in degrees = %.8f\n", range);
# find average section azimuth and compute heading perpendicular to it
numeric heading;
sectionAzimuth = 0;
for i = 1 to numSections {
	sectionAzimuth += sectionAzimuths[i];
	printf("Cumulative section azimuth = %.1f\n", sectionAzimuth);
	}
heading = (sectionAzimuth / numSections) - 90;
if (heading < 0) then heading += 360;
printf("Computed LookAt heading = %.1f\n\n", heading);
 
# create <LookAt> element in KML to set viewer settings
class XMLNODE lookat;
lookat = folder.NewChild("LookAt");
lookat.NewChild("longitude", sprintf("%.8f", sectionsExtents.center.x) ); 
lookat.NewChild("latitude", sprintf("%.8f", sectionsExtents.center.y) );
lookat.NewChild("altitude", sprintf("%.1f", maxZm) );
lookat.NewChild("altitudeMode", "absolute");
lookat.NewChild("tilt", "70");
lookat.NewChild("heading", sprintf("%.1f", heading) );
lookat.NewChild("range", sprintf("%.1f", range) );
############## Write the KML structure to a file and add to KMZ
kmldoc.Write(kmlfile$);
kmz.AddFile(kmzFilepath, kmlFilepath, "doc.kml");
DeleteFile(kmlfile$);
print("End.");