# ExportTerrainCollada.sml
# Processes a selected DEM and overlapping IMAGE into a 3D model incorporating
# a custom terrain surface and matching image overlay that can be viewed in 3D
# in Google Earth. DEM and IMAGE extents do not have to match exactly, and they
# can have different coordinate reference systems. The DEM and IMAGE are automatically
# cropped to their common extents if necessary.
# The surface and image texture are output as a COLLADA model packaged in
# a KMZ file (compressed KML). The COLLADA model is referenced as a placemark
# in the KML file stored within the KMZ. The placemark includes the WGS84 Geographic
# coordinates of the origin (center point) of the model, its height in meters
# above sea level, and its orientation (angle to north). This placemark information
# allows Google Earth to correctly locate and orient the terrain model in 3D space.
# The COLLADA model specifies a triangular mesh defining the custom 3D terrain
# derived from the DEM and a mapping to corresponding normalized image
# coordinates to allow the image (exported as a JPEG, PNG, or TIFF file) to be overlaid as
# a texture onto the triangular mesh in Google Earth. The triangular mesh coordinates
# are read from a TIN created from the DEM. The COLLADA model has projected coordinates
# derived from the coordinate reference system of the IMAGE (unless it has a
# non-projected CRS, in which case a location-appropriate UTM CRS is created and used).
# The custom script dialog provides interactive selection of the input DEM and IMAGE
# and the output KMZ file. The elevation units for the DEM (meters or feet) must
# be correctly specified to provide proper vertical scaling for the model. A
# vertical offset of the model above the Google Earth terrain surface can be
# set to avoid intersection of the model and the Google Earth terrain (portions
# of the model below the Google Earth terrain are not visible).
# NOTE: MODELS WITH IMAGE TEXTURES LARGER THAN ABOUT 1024 X 1024 OR WITH HIGHLY
# DETAILED TERRAIN MESHES MAY LOAD VERY SLOWLY IN GOOGLE EARTH AND DEGRADE GOOGLE
# EARTH PERFORMANCE.
# Randy Smith, MicroImages, Inc.
# 28 February 2011
class RVC_RASTER DEM, IMAGE;
class RVC_OBJITEM demObjItem, imgObjItem;
class STRING prompt$, err$;
class RVC_GEOREFERENCE georefDEM, georefIMG;
class SR_COORDREFSYS crsDEM, crsIMG;
class POINT2D demCellSizesPt, imgCellSizesPt;
numeric demCellSize, imgCellSize;
numeric err;
numeric i, j, k, m; # loop counters
numeric tempfile1exists, tempfile2exists, tempfile3exists;
class STRING zUnit$; # length unit for DEM elevation values
numeric zScaleDEM;
numeric zOffset;
class STRING kmzfilename$;
class GUI_DLG dlgwin;
class GUI_CTRL_COMBOBOX elevunitsCombo, formatCombo;
class GUI_CTRL_PUSHBUTTON imgBtn, kmzBtn;
class GUI_CTRL_EDIT_STRING demText, imgText, kmzText;
class STRING outFormat$;
numeric canceled;
############# STATUS DIALOG COMPONENTS #####################
class STATUSDIALOG statusD; # status dialog to show process status
class STATUSCONTEXT statusC; # context for the status dialog
############################################################################
#######
####### USER-DEFINED FUNCTIONS AND PROCEDURES
# procedure automatically called when the dialog window is initialized
proc onInitDialog()
{
# disable the OK button
dlgwin.SetOkEnabled(0);
# add elevation units to the DEM Elevation Units combobox
elevunitsCombo = dlgwin.GetCtrlByID("elevunitsCombo");
elevunitsCombo.AddItem("meters", "meters");
elevunitsCombo.AddItem("feet", "feet");
elevunitsCombo.SetSelectedItemIndex(0); # set meters as default unit shown
# add formats to the Image Format combobox
formatCombo = dlgwin.GetCtrlByID("formatCombo");
formatCombo.AddItem("JPEG", "JPEG");
formatCombo.AddItem("PNG", "PNG");
formatCombo.AddItem("TIFF", "TIFF");
formatCombo.SetSelectedItemIndex(0); # set JPEG as the default format
# get handles for pushbuttons to enable them
imgBtn = dlgwin.GetCtrlByID("imgBtn");
kmzBtn = dlgwin.GetCtrlByID("kmzBtn");
# get handles for the edittext controls to show the names of the input rasters in the dialog
demText = dlgwin.GetCtrlByID("demText");
imgText = dlgwin.GetCtrlByID("imgText");
kmzText = dlgwin.GetCtrlByID("kmzText");
}
# procedure called when DEM button is pressed on the dialog window
# Get ObjItem for input DEM and write its display string to its edittext
proc getDEM()
{
prompt$ = "Choose input DEM raster:";
DlgGetObject(prompt$, "Raster", demObjItem, "ExistingOnly");
demText.SetValueStr(demObjItem.GetDisplayString() );
imgBtn.SetEnabled(1);
}
# procedure called when IMAGE button is pressed on the dialog window
# Get ObjItem for input IMAGE and write its display string to its edittext
proc getIMAGE()
{
prompt$ = "Choose input IMAGE raster:";
DlgGetObject(prompt$, "Raster", imgObjItem, "ExistingOnly");
imgText.SetValueStr(imgObjItem.GetDisplayString() );
kmzBtn.SetEnabled(1);
}
# procedure called when KMZ button is pressed on the dialog window
# Get filename/location for the output KMZ file.
proc getKMZ()
{
prompt$ = "Choose an output KMZ file to package the image/terrain models:";
kmzfilename$ = GetOutputFileName(_context.ScriptDir, prompt$, ".kmz");
kmzText.SetValueStr(kmzfilename$);
printf("Full path to KMZ file: %s\n", kmzfilename$);
dlgwin.SetOkEnabled(1);
}
# procedure called when OK button on the dialog window is pressed.
# Get values for maximum image tile size, DEM elevation units, and Z offset
proc onDlgOK()
{
zUnit$ = elevunitsCombo.GetValueStr();
zOffset = dlgwin.GetCtrlByID("zOffsetNum").GetValueNum();
outFormat$ = formatCombo.GetValueStr();
# get scale factor to convert DEM elevation units to meters
zScaleDEM = GetUnitConvDist(zUnit$, "meters");
printf("DEM elevation unit = %s\n", zUnit$);
printf("Scale factor from DEM elevation units to meters = %.4f\n", zScaleDEM);
printf("Offset of terrain model above true elevation = %d m\n", zOffset);
}
# procedure called when Cancel button on the dialog window is pressed.
# Exit script.
proc onDlgCancel()
{
dlgwin.Close(1);
canceled = 1;
}
############################
# error checking procedure for image pipeline routines
proc ReportError(numeric linenum, numeric err)
{
printf("FAILED -line: %d, err: %d\n", linenum, err);
PopupError(err);
}
#############################
# function to get and check the georeference and coordinate reference system of the input DEM and IMAGE
func checkGeoref(
class RVC_RASTER RAST,
class STRING name$,
var class RVC_GEOREFERENCE georef,
var class SR_COORDREFSYS crs,
var class STRING e$)
{
if (RAST.GetDefaultGeoref(georef) )
{
crs = georef.GetCoordRefSys();
if (crs.IsDefined() )
{
if (crs.IsLocal() <> 1)
{
# printf("%s Coordinate Reference System = %s\n", name$, crs.Name);
return 1;
}
else
{
e$ = sprintf("Selected %s has a local arbitrary coordinate reference system. Exiting now.\n", name$);
return 0;
}
}
else
{
e$ = sprintf("Selected %s has an undefined coordinate reference system. Exiting now.\n", name$);
return 0;
}
}
else
{
e$ = sprintf("Selected %s is not georeferenced. Exiting now.\n", name$);
return 0;
}
} # end checkGeoref()
######################################
# function to return the smaller of two values
func smaller(numeric a, numeric b)
{
if (a < b) then
return a;
else return b;
}
###############################################
# function to create a UTM CRS for raster that has lat / lon CRS
func class SR_COORDREFSYS makeUTMcrs(class RVC_RASTER RAST, class RVC_GEOREFERENCE inGeoref)
{
local class RECT3D extentsObj;
local class POINT2D center;
local class STRING demCRSunit$;
local class TRANS2D_MAPGEN transObjToMap;
inGeoref.GetTransParm(transObjToMap, 0, inGeoref.GetCalibModel() );
RAST.GetExtents(extentsObj);
center = extentsObj.center;
center = transObjToMap.ConvertPoint2DFwd(center);
printf("Image center: longitude = %.6f, latitude = %.6f\n", center.x, center.y);
# variables for setting up UTM CRS
local class SR_COORDREFSYS utmCRS; # the CRS to create
local class SR_COORDSYS utmCoordSys; # coordinate system
local class SR_COORDOPDEF projectDef; # projection parameters and method
local class SR_COORDOPMETHOD method;
local class SR_DATUM datum;
datum = inGeoref.GetCoordRefSys().Datum; # get datum from the input raster
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
local numeric zoneNum, centMerid; # UTM zone and longitude of central meridian
zoneNum = ceil( (180 + center.x) / 6 );
printf("UTM zone number = %d\n", zoneNum);
centMerid = (zoneNum * 6) - 183;
projectDef.SetParmValueNum("LongitudeOfNaturalOrigin", centMerid);
# create defined UTM coordinate reference system
utmCRS.Create(utmCoordSys, datum, projectDef);
printf("CRS for reprojection: %s\n", utmCRS.name);
printf("Central meridian = %d\n", centMerid);
return utmCRS;
} # end makeUTMcrs()
######################################################
# function to resample a raster to a specified CRS and cell size
func class RVC_OBJITEM resampleRast(
class RVC_OBJITEM objItemS, # objItem of input raster
class RVC_OBJITEM objItemT, # objItem of output raster
class SR_COORDREFSYS crs, # CRS to resample to
numeric cellsize )
{
local class IMAGE_PIPELINE_SOURCE_RVC sourceRVC(objItemS);
err = sourceRVC.Initialize();
if (err < 0) ReportError(_context.CurrentLineNum, err);
local class IMAGE_PIPELINE_FILTER_RESAMPLE resample(sourceRVC, crs, cellsize, cellsize, "CubicConvolution");
err = resample.Initialize();
if (err < 0) ReportError(_context.CurrentLineNum, err);
local class IMAGE_PIPELINE_TARGET_RVC targetRVC(resample, objItemT);
targetRVC.SetCompression("DPCM", 0);
err = targetRVC.Initialize();
if (err < 0) ReportError(_context.CurrentLineNum, err);
targetRVC.Process();
objItemT = targetRVC.GetObjItem();
return objItemT;
} # end resampleRast()
#########################################################
# function to crop a raster
func class RVC_OBJITEM cropImage(
class RVC_OBJITEM objItemS, # objItem of input raster
class RECT intersect, # common image extents RECT to use to crop input raster
class RVC_OBJITEM objItemT) # objItem of output
{
local class IMAGE_PIPELINE_SOURCE_RVC sourceRVC(objItemS);
err = sourceRVC.Initialize();
if (err < 0) ReportError(_context.CurrentLineNum, err);
local class RECT sourceRect;
sourceRect.pt1.x = 0; sourceRect.pt1.y = 0;
sourceRect.pt2.x = sourceRVC.GetTotalColumns();
sourceRect.pt2.y = sourceRVC.GetTotalRows();
if (sourceRect.ContainsRect(intersect) )
{
print("Intersect extents are fully contained within DEM extents.");
local class IMAGE_PIPELINE_FILTER_CROP crop(sourceRVC, intersect);
err = crop.Initialize();
if (err < 0) ReportError(_context.CurrentLineNum, err);
local class IMAGE_PIPELINE_TARGET_RVC targetRVC(crop, objItemT);
err = targetRVC.Initialize();
if (err < 0) ReportError(_context.CurrentLineNum, err);
targetRVC.Process();
objItemT = targetRVC.GetObjItem();
return objItemT;
}
else
{
print("Intersect extents not fully contained within DEM extents. Cannot crop DEM.");
}
} # end cropImage()
################################################################
# 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
######################
clear(); # clear the console window
class STRING dlgxml$;
dlgxml$ = '<?xml version="1.0"?>
<!DOCTYPE root SYSTEM "smlforms.dtd">
<root>
<dialog id="dlgid" Title="Export DEM and IMAGE to KML/COLLADA" OnOK="onDlgOK()" OnCancel="onDlgCancel()">
<groupbox Name=" Select Inputs " HorizResize="Expand" ExtraBorder="3">
<pane Orientation="horizontal">
<pushbutton id="demBtn" Name=" DEM... " OnPressed="getDEM()"/>
<edittext id="demText" Width="40" ReadOnly="true"/>
</pane>
<pane Orientation="horizontal">
<pushbutton id="imgBtn" Name=" IMAGE... " Enabled="false" OnPressed="getIMAGE()"/>
<edittext id="imgText" Width="40" ReadOnly="true"/>
</pane>
</groupbox>
<groupbox Name=" Select Output File " HorizResize="Expand" ExtraBorder="3">
<pane Orientation="horizontal">
<pushbutton id="kmzBtn" Name = " KMZ... " Enabled="false" OnPressed="getKMZ()"/>
<edittext id="kmzText" Width="40" ReadOnly="true"/>
</pane>
</groupbox>
<pane Orientation="horizontal">
<label>DEM elevation units: </label>
<combobox id="elevunitsCombo" HorizResize="Fixed"/>
<label>Image Format: </label>
<combobox id="formatCombo" HorizResize="Fixed"/>
</pane>
<pane Orientation="horizontal">
<label>Offset of model above true elevation: </label>
<editnumber id="zOffsetNum" Width="4" Precision="0" Default="20" MinVal="0" HorizResize="Fixed"/>
<label> meters</label>
</pane>
</dialog>
</root>';
##############################################
### parse XML text for the dialog into memory;
### return an error code (number < 0 ) if there are syntax errors
class XMLDOC dlgdoc;
err = dlgdoc.Parse(dlgxml$);
if ( err < 0 ) {
PopupError( err ); # Popup an error dialog. "Details" button shows syntax errors.
Exit();
}
################################################
# get the dialog element from the parsed XML document and
# show error message if the dialog element can't be found
class XMLNODE dlgnode;
dlgnode = dlgdoc.GetElementByID("dlgid");
if ( dlgnode == 0 ) {
PopupMessage("Could not find dialog node in XML document");
Exit();
}
################################################
# Set the XML dialog element as the source for the GUI_DLG class
# instance we are using for the dialog window.
err = dlgwin.SetXMLNode(dlgnode);
if ( err < 0 ) {
PopupError( err ); # Popup an error dialog.
Exit();
}
dlgwin.SetOnInitDialog(onInitDialog);
#############################
# open the dialog window
dlgwin.DoModal();
if (canceled) then Exit();
####################
# Open the input rasters using the selected ObjItems
DEM.Open(demObjItem, "Read");
IMAGE.Open(imgObjItem, "Read");
# check DEM and IMAGE georeference and CRS
if ( checkGeoref(DEM, "DEM", georefDEM, crsDEM, err$ ) <> 1 or checkGeoref(IMAGE, "IMAGE", georefIMG, crsIMG, err$) <> 1 )
{
PopupMessage(err$);
Exit();
}
printf("DEM Coordinate Reference System = %s\n", crsDEM.Name);
printf("IMAGE Coordinate Reference System = %s\n", crsIMG.Name);
#########################
# Set up KMZ file and ZIP file structure
class FILEPATH kmzFilepath(kmzfilename$);
class FILEPATH dirFilepath(kmzfilename$);
dirFilepath.StripToFolder();
class STRING dir$;
dir$ = dirFilepath;
printf("Directory containing KMZ file: %s\n", dir$);
# ZIP file structure to use for KMZ
class ZIPFILE kmz;
########################################
# get cell sizes for the DEM and IMAGE
DEM.ComputeScaleFromGeoref(georefDEM, demCellSizesPt, 1);
printf("Cell sizes of DEM: col = %.1f m, lin = %.1f m\n", demCellSizesPt.x, demCellSizesPt.y);
demCellSize = smaller(demCellSizesPt.x, demCellSizesPt.y);
printf("Smaller cell size of DEM: %.1f\n", demCellSize);
IMAGE.ComputeScaleFromGeoref(georefIMG, imgCellSizesPt, 1);
printf("Cell sizes of IMAGE: col = %.1f m, lin = %.1f m\n", imgCellSizesPt.x, imgCellSizesPt.y);
imgCellSize = smaller(imgCellSizesPt.x, imgCellSizesPt.y);
printf("Smaller cell size of IMAGE: %.1f\n", imgCellSize);
##########################################################
### Create status dialog and its context
##########################################################
statusD.Create();
statusC = statusD.CreateContext();
### check IMAGE for projected CRS; if not projected, it is geographic (lat / lon);
### reproject to a UTM CRS with zone appropriate for its location
class SR_COORDREFSYS crsIMGp;
if (crsIMG.IsProjected() <> 1)
{
tempfile1exists = 1;
class RVC_OBJECT tempfile1;
tempfile1.MakeTempFile(1); # set to delete on close
class RVC_OBJITEM imgpObjItem; # for image converted to projected CRS
class RVC_DESCRIPTOR imgpDescript;
imgpDescript.SetName("IMAGE");
imgpDescript.SetDescription("");
imgpObjItem.CreateNew(tempfile1.GetObjItem(), "RASTER", imgpDescript);
# compute a UTM CRS appropriate to center point and reproject inage raster to UTM
print("Reprojecting IMAGE with Geographic CRS to UTM.");
statusC.Message = "Reprojecting IMAGE with Geographic CRS to UTM.";
# call function to create a UTM CRS, computing the zone from the center of the image extents
crsIMGp = makeUTMcrs(IMAGE, georefIMG);
# call user-defined procedure to resample raster to specified CRS and cell size
imgpObjItem = resampleRast(IMAGE.GetObjItem(), imgpObjItem, crsIMGp, imgCellSize);
# close the original IMAGE raster and open the resampled raster as IMAGE
IMAGE.Close();
IMAGE.Open(imgpObjItem, "Read");
georefIMG.OpenLastUsed(IMAGE);
crsIMG = crsIMGp;
}
else
print("IMAGE has a projected CRS, no reprojection necessary.");
#############################################
### if necessary, reproject DEM to match CRS of IMAGE
if (crsDEM.GetID("MicroImages") <> crsIMG.GetID("MicroImages") )
{
tempfile2exists = 1;
class RVC_OBJECT tempfile2;
tempfile2.MakeTempFile(1); # set to delete on close
class RVC_OBJITEM dempObjItem;
class RVC_DESCRIPTOR dempDescript;
class SR_COORDREFSYS crsDEMp;
dempDescript.SetName("DEMp");
dempDescript.SetDescription("");
dempObjItem.CreateNew(tempfile2.GetObjItem(), "RASTER", dempDescript);
print("Reprojecting DEM to match CRS of IMAGE.");
statusC.Message = "Reprojecting DEM to match CRS of IMAGE.";
# call user-defined procedure to resample raster to specified CRS and cell size
dempObjItem = resampleRast(demObjItem, dempObjItem, crsIMG, demCellSize);
# close the original DEM raster and open the resampled raster as DEM
DEM.Close();
DEM.Open(dempObjItem, "Read");
georefDEM.OpenLastUsed(DEM);
}
else
print("DEM and IMAGE have the same CRS, DEM does not need to be reprojected.");
#################################################
# get extents of IMAGE and DEM in map coordinates
class RECT3D extentsIMG, extentsDEM; # get extents of each raster in object coordinates
class RECT extentsRECTimg, extentsRECTdem; # ConvertRectFwd method below returns a RECT, not RECT3D
IMAGE.GetExtents(extentsIMG);
DEM.GetExtents(extentsDEM);
# get coordinate transformations from image to map coordinates
class TRANS2D_MAPGEN transIMGobjToMap, transDEMobjToMap;
georefIMG.GetTransParm(transIMGobjToMap, 0, georefIMG.GetCalibModel() );
georefDEM.GetTransParm(transDEMobjToMap, 0, georefDEM.GetCalibModel() );
# convert extents to map coordinates
extentsRECTimg = transIMGobjToMap.ConvertRectFwd(extentsIMG);
extentsRECTdem = transDEMobjToMap.ConvertRectFwd(extentsDEM);
printf("IMAGE extents: x1 = %.1f, y1 = %.1f, x2 = %.1f, y2 = %.1f\n", extentsRECTimg.pt1.x, extentsRECTimg.pt1.y, extentsRECTimg.pt2.x, extentsRECTimg.pt2.y);
printf("DEM extents: x1 = %.1f, y1 = %.1f, x2 = %.1f, y2 = %.1f\n", extentsRECTdem.pt1.x, extentsRECTdem.pt1.y, extentsRECTdem.pt2.x, extentsRECTdem.pt2.y);
# check if extents of DEM and IMAGE are different
if (extentsRECTimg.pt1 != extentsRECTdem.pt1 or extentsRECTimg.pt2 != extentsRECTdem.pt2)
{
# check that extents overlap
if (extentsRECTimg.Overlaps(extentsRECTdem) )
{
# find intersect of the two extents RECTs in map coordinates
class RECT intersect; # new extents RECT
intersect.pt1 = extentsRECTimg.pt1; intersect.pt2 = extentsRECTimg.pt2; # set to same extents as IMAGE
intersect.Intersect(extentsRECTdem); # intersect with DEM extents
printf("common extents: x1 = %.1f, y1 = %.1f, x2 = %.1f, y2 = %.1f\n", intersect.pt1.x, intersect.pt1.y, intersect.pt2.x, intersect.pt2.y);
class RVC_OBJECT tempfile3;
tempfile3.MakeTempFile(1); # set to delete on close
tempfile3exists = 1;
numeric cropIMAGE = 1; # flags to indicate whether to crop, set to 1 by default
numeric cropDEM = 1;
# check for full containment cases where don't need to crop one of the inputs
if (extentsRECTdem.ContainsRect(extentsRECTimg) ) then
cropIMAGE = 0; # if DEM contains IMAGE, don't need to crop IMAGE
else if (extentsRECTimg.ContainsRect(extentsRECTdem) ) then
cropDEM = 0; # if IMAGE contains DEM, don't need to crop DEM
class RVC_DESCRIPTOR iDescript;
if (cropDEM)
{
class RVC_OBJITEM demiObjItem;
iDescript.SetName("DEM"); iDescript.SetDescription("");
demiObjItem.CreateNew(tempfile3.GetObjItem(), "RASTER", iDescript);
print("Cropping the DEM to the common extents.");
statusC.Message = "Cropping the DEM to the common extents.";
# convert intersect RECT to object coordinates and adjust to integer values
intersect = transDEMobjToMap.ConvertRectInv(intersect);
intersect.pt1.x = ceil(intersect.pt1.x);
intersect.pt1.y = ceil(intersect.pt1.y);
intersect.pt2.x = floor(intersect.pt2.x);
intersect.pt2.y = floor(intersect.pt2.y);
demiObjItem = cropImage(DEM.GetObjItem(), intersect, demiObjItem);
DEM.Close();
DEM.Open(demiObjItem, "Read");
georefDEM.OpenLastUsed(DEM);
georefDEM.GetTransParm(transDEMobjToMap, 0, georefDEM.GetCalibModel() );
}
if (cropIMAGE)
{
class RVC_OBJITEM imgiObjItem;
iDescript.SetName("IMAGE"); iDescript.SetDescription("");
imgiObjItem.CreateNew(tempfile3.GetObjItem(), "RASTER", iDescript);
print("Cropping the IMAGE to the common extents.");
statusC.Message = "Cropping the IMAGE to the common extents.";
# convert intersect RECT to object coordinates and adjust to integer values
intersect = transIMGobjToMap.ConvertRectInv(intersect);
intersect.pt1.x = ceil(intersect.pt1.x);
intersect.pt1.y = ceil(intersect.pt1.y);
intersect.pt2.x = floor(intersect.pt2.x);
intersect.pt2.y = floor(intersect.pt2.y);
imgiObjItem = cropImage(IMAGE.GetObjItem(), intersect, imgiObjItem);
IMAGE.Close();
IMAGE.Open(imgiObjItem, "Read");
georefIMG.OpenLastUsed(IMAGE);
georefIMG.GetTransParm(transIMGobjToMap, 0, georefIMG.GetCalibModel() );
}
if (tempfile2exists) then tempfile2.Close();
}
else
{
print("IMAGE and DEM do not overlap. Exiting now.");
SetStatusMessage("IMAGE and DEM do not overlap. Exiting now.");
Exit();
}
}
else # DEM and IMAGE have the same extents
print("DEM and IMAGE have the same extents, no cropping necessary.");
printf("Final DEM dimensions: %d columns, %d lines\n", DEM.$Info.NumCols, DEM.$Info.NumLins);
printf("Final IMAGE dimensions: %d columns, %d lines\n", IMAGE.$Info.NumCols, IMAGE.$Info.NumLins);
#################################
### get unit in DEM's projected coordinate reference system
### and set conversion factors between this unit and meters
class STRING unit$;
class SR_COORDAXIS axisX;
numeric unitConvToMeters;
numeric metersConvToUnit;
crsDEM = georefDEM.GetCoordRefSys();
axisX = crsDEM.CoordSys.GetAxis();
unit$ = axisX.Unit.GetName();
unitConvToMeters = GetUnitConvDist(unit$, "meter");
printf("Map unit = %s, conversion factor to meters = %.4f\n", unit$, unitConvToMeters);
metersConvToUnit = GetUnitConvDist("meter", unit$);
# set up a WGS84 / Geographic CRS for lat/lon coordinates to position terrain models in Google Earth
class SR_COORDREFSYS geoCRS;
geoCRS.Assign("Geographic2D_WGS84_Deg");
printf("Geographic CRS for Google Earth = %s\n", geoCRS.Name);
# set up coordinate transformation from DEM CRS to geoCRS
class TRANS2D_MAPGEN transToGeo;
transToGeo.InputCoordRefSys = crsDEM;
transToGeo.OutputCoordRefSys = geoCRS;
# Find minimum values for the DEM map coordinates and elevation to use to set coordinate offsets
# in the COLLADA models and to position them in the KMZ file
numeric minX, minY, minZ, minZm;
minX = extentsDEM.pt1.x; minY = extentsDEM.pt1.y;
minZ = GlobalMin(DEM); # minimum elevation value in DEM elevation units
minZm = minZ * zScaleDEM; # minimum elevation in meters to set origin in Google Earth
printf("Minimum DEM z value in native elevation units = %.1f\n", minZ);
printf("Minimum DEM z value in neters = %.1f\n", minZm);
######################################
# string with COLLADA file template
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="meters" meter="1.0"/>
<up_axis>Z_UP</up_axis>
</asset>
<library_images>
<image id="terrain-image" name="terrain-image">
<init_from></init_from>
</image>
</library_images>
<library_materials>
<material id="terrainID" name="terrain">
<instance_effect url="#terrain-effect"/>
</material>
</library_materials>
<library_effects>
<effect id="terrain-effect" name="terrain-effect">
<profile_COMMON>
<newparam sid="terrain-image-surface">
<surface type="2D">
<init_from>terrain-image</init_from>
</surface>
</newparam>
<newparam sid="terrain-image-sampler">
<sampler2D>
<source>terrain-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="terrain-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="terrain" 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="TerrainScene" name="TerrainScene">
<node id="Model" name="Model">
<node id="mesh1" name="mesh1">
<instance_geometry url="#mesh1-geometry">
<bind_material>
<technique_common>
<instance_material symbol="terrain" target="#terrainID">
<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="#TerrainScene"/>
</scene>
</COLLADA>';
######################################
### string with KML file template
class STRING kmldoc$;
class STRING kmlfile$ = dir$ + "/" + "doc.kml"; # filepath to 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>0</open>
</Folder>
</kml>';
# variables for modifying the KML file template
class XMLDOC kmldoc;
class XMLNODE kmlRootNode;
class XMLNODE folder, document, placemark, model;
class XMLNODE node, node1, node2, node3;
#########################################
# parse the KML file template to reference and position the COLLADA models
err = kmldoc.Parse(kmldoc$);
kmlRootNode = kmldoc.GetRootNode(); # get root node
folder = kmlRootNode.FindChild("Folder"); # get folder node that will hold a Document for each cross-section
# set the name of the folder to show in the Google Earth layer list
node = folder.FindChild("name");
node.SetText(kmzfilename$);
###############################################################
################# make a JEPG file from the drape raster
class RVC_DESCRIPTOR imgDescriptor;
imgDescriptor = IMAGE.GetDescriptor();
class STRING imgName$ = kmzFilepath.GetNameOnly();
class STRING ext$; # file extension for image file
class STRING img$ = dir$ + "/" + imgName$;
class FILEPATH imgFilepath(img$); # filepath for texture image, lacking extension
statusC.Message = "Making image file for model texture.";
# PIPELINE SOURCE
class IMAGE_PIPELINE_SOURCE_RVC sourceIMG(IMAGE.GetObjItem() );
err = sourceIMG.Initialize();
if ( err < 0) ReportError(_context.CurrentLineNum, err);
# PIPELINE TARGET
class IMAGE_PIPELINE_TARGET target; # generic pipeline target, to be assigned specific type
switch (outFormat$)
{
case "JPEG":
ext$ = ".jpg";
imgFilepath.SetExtension(ext$);
class IMAGE_PIPELINE_TARGET_JPEG targetJPG(sourceIMG, imgFilepath, 75, "none");
target = targetJPG;
break;
case "PNG":
ext$ = ".png";
imgFilepath.SetExtension(ext$);
class IMAGE_PIPELINE_TARGET_PNG targetPNG(sourceIMG, imgFilepath, "none");
target = targetPNG;
break;
case "TIFF":
ext$ = ".tif";
imgFilepath.SetExtension(ext$);
class IMAGE_PIPELINE_TARGET_TIFF_SETTINGS tiff_settings;
tiff_settings.SetCompression("JPEG_DCT");
class IMAGE_PIPELINE_TARGET_TIFF targetTIFF(sourceIMG, imgFilepath, "none");
targetTIFF.SetParms(tiff_settings);
target = targetTIFF;
break;
}
# Initialize the pipeline target
err = target.Initialize();
if ( err < 0) ReportError(_context.CurrentLineNum, err);
# Execute the pipeline
err = target.Process(1); # show status bar in status dialog
if ( err < 0) ReportError(_context.CurrentLineNum, err);
# add the JPEG file to the KMZ (ZIP) file
kmz.AddFile(kmzFilepath, imgFilepath, "files/" + imgName$ + ext$);
# delete the source JPEG file
DeleteFile(imgFilepath);
###########################################################
################ make TIN from DEM to generate COLLADA meshes
class RVC_OBJECT tempfile;
tempfile.MakeTempFile(1); # set to delete on close
class FILEPATH tempfilePath(tempfile.GetObjItem().GetFilePath());
class RVC_TIN demTIN;
class RVC_GEOREFERENCE georefTIN;
CreateTIN(demTIN, tempfilePath, "TIN");
numeric zTolerance = 1;
numeric maxNodes = 100000;
numeric zDelta = 1;
numeric shortEdge = 3;
numeric longEdge = 50000;
print("\nMaking TIN from DEM...");
statusC.Message = "Making TIN from DEM";
RasterToTINIterative2(demTIN, DEM, zTolerance, maxNodes, zDelta, shortEdge, longEdge);
demTIN.OpenAttached("Read"); # RasterToTINInterative2 closes the newly-created TIN, so must open it
demTIN.GetDefaultGeoref(georefTIN);
class TRANS2D_MAPGEN transTinObjToMap;
georefTIN.GetTransParm(transTinObjToMap, 0, georefTIN.GetCalibModel() );
# set origin point for the COLLADA model at the center of the TIN
class RECT3D extentsTIN;
demTIN.GetExtents(extentsTIN);
class POINT3D originPt;
originPt = extentsTIN.center;
originPt = transTinObjToMap.ConvertPoint3DFwd(originPt); # convert to projected map coordinates
printf("Origin point for tile in projected coordinates: x = %.1f, y = %.1f\n", originPt.x, originPt.y);
# find angle to north at the origin point to orient the COLLADA model in Google Earth
numeric angleToNorth = georefTIN.GetCoordRefSys().ComputeAngleToNorth(originPt);
# reproject origin point to WGS84 / Geographic CRS for placing model in Google Earth
originPt = transToGeo.ConvertPoint3DFwd(originPt);
originPt.z = minZm;
printf("Origin point for tile in WGS84/Geographic coordinates: x = %.8f, y = %.8f\n", originPt.x, originPt.y);
printf("Origin point z value (meters): %.1f\n", originPt.z);
printf("Origin point z value in DEM elevation units: %.1f\n", originPt.z / zScaleDEM);
#############################################################
############## process TIN to COLLADA meshes
class POINT3D pt;
class STRING xyzString, uvString;
numeric s, t;
numeric numPts = demTIN.$Info.NumNodes;
numeric numTriangles = demTIN.$Info.NumTriangles;
numeric k;
# TIN object coordinates are natively offset to a local origin at lower left corner, and may be scaled as well.
# For TIN created from DEM, node in lower left corner cell is at middle of cell and gets TIN object
# coordinates 0.5, 0.5. Origin of TIN object coordinates is at lower left corner of the lower left cell of the DEM.
# Get scales for TIN x and y coordinates (as a POINT2D) from its georeference
class POINT2D mapScalesTIN;
demTIN.ComputeScaleFromGeoref(georefTIN, mapScalesTIN, 0);
numeric zScaleTIN = demTIN.GetZScale(); # get Z scale and offset from TIN object info
numeric zOffsetTIN = demTIN.GetZOffset();
printf("TIN result: number of nodes = %d, number of triangles = %d\n", numPts, numTriangles);
printf("TIN scales: x = %.4f, %.4f, %.4f\n", mapScalesTIN.x, mapScalesTIN.y, zScaleTIN);
class POINT2D ptList[numPts]; # array of POINT2Ds to hold drape image coordinates for each node;
statusC.Message = sprintf("Processing %d TIN nodes...", numPts);
statusC.BarInit(numPts, 0);
# loop through TIN nodes to get x, y, and z values and write to xyzString, find corresponding normalized
# image coordinates and write them to uvString.
for k = 1 to numPts
{
statusC.BarIncrement(1, 0);
# get x, y, and z values (object coordinates) from Node table.
pt.x = demTIN.Node[ k ].NODE.x;
pt.y = demTIN.Node[ k ].NODE.y;
pt.z = demTIN.Node[ k ].NODE.z * zScaleTIN + zOffsetTIN; # pt elevation in DEM's elevation units
pt.z = pt.z * zScaleDEM - originPt.z; # rescale elevations to meters and subtract origin point elevation
# apply scale factors from georeference to TIN x and y object coordinates
# so they have same scale as map coordinates and add to xyzString for COLLADA model
class POINT2D meshPt;
meshPt.x = (pt.x - extentsTIN.center.x) * mapScalesTIN.x;
meshPt.y = (pt.y - extentsTIN.center.y) * mapScalesTIN.y;
xyzString += sprintf("%.2f %.2f %.1f ", meshPt.x, meshPt.y, pt.z);
# printf("xyzString = %s\n", xyzString);
# get map coordinates of the TIN node from its raw (unscaled) object coordinates using previously-obtained
# object-to-map coordinate transformation
pt = transTinObjToMap.ConvertPoint3DFwd(pt);
# find the image coordinates of the point
pt = transIMGobjToMap.ConvertPoint3DInv(pt);
# add image coordinates to point list (for determining node order for triangles in later step)
ptList[ k ] = pt;
# convert image coordinates to texture coordinates (normalized to 0 to 1 range, origin in lower left)
s = pt.x / (IMAGE.$Info.NumCols);
t = 1 - (pt.y / IMAGE.$Info.NumLins); # inverted due to lower left origin
uvString += sprintf("%.6f %.6f ", s, t);
} # end for k = 1 to numPts
# loop through TIN triangles to get the node index for each triangle vertex and reorder to
# counter-clockwise order if necessary
array numeric triNodeIndexArray[3];
array numeric tempArray[3];
numeric testCCW;
class STRING pString;
statusC.Message = sprintf("Processing %d TIN triangles...", numTriangles);
statusC.BarInit(numTriangles, 0);
for k = 1 to numTriangles
{
statusC.BarIncrement(1, 0);
tempArray[1] = demTIN.Triangle[ k ].TRIANGLE.Node1;
tempArray[2] = demTIN.Triangle[ k ].TRIANGLE.Node2;
tempArray[3] = demTIN.Triangle[ k ].TRIANGLE.Node3;
testCCW = ccwImageCoords(ptList[ tempArray[1] ], ptList[ tempArray[2] ], ptList[ tempArray[3] ]);
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];
}
# 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 are
# node index - 1.
for m = 1 to 3
{
pString += sprintf("%d %d %d ", triNodeIndexArray[ m ] - 1, 0, triNodeIndexArray[ m ] - 1);
}
} # end for k = 1 to numTriangles
# close objects and temporary files that are no longer needed
demTIN.Close();
DEM.Close();
IMAGE.Close();
tempfile.Close();
#######################################
# get current date and time for inclusion in the COLLADA file
class DATETIME dt;
dt.SetCurrent();
class STRING dt$, format$;
format$ = "%Y-%m-%dT%H:%M:%SZ";
dt$ = dt.GetString(format$);
printf("formatted date and time: %s\n", dt$);
##############################################
########### make the COLLADA file
print("Making COLLADA file.");
### parse the COLLADA template
class XMLDOC daedoc;
err = daedoc.Parse(dae$);
class XMLNODE colladaRootNode;
colladaRootNode = daedoc.GetRootNode();
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 JPEG file as the library image
node = daedoc.GetElementByID("terrain-image");
node1 = node.FindChild("init_from");
node1.SetText(imgFilepath.GetName() );
### set the position mesh count and array
numeric 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$ + "/" + imgName$ + ".dae";
printf("\nDAE filepath = %s\n", daefile$);
daedoc.Write(daefile$);
class FILEPATH daeFilepath(daefile$);
########################
# add COLLADA file to the KMZ
kmz.AddFile(kmzFilepath, daeFilepath, "files/" + imgName$ + ".dae");
# delete the source DAE file
DeleteFile(daefile$);
####################################################
#####
##### add "document" to the doc.kml file to reference and place this model
document = folder.NewChild("Document");
node = document.NewChild("name");
node.SetText(imgName$);
placemark = document.NewChild("Placemark");
# create name element for placemark
placemark.NewChild("name", imgName$ );
# create style element for placemark
node = placemark.NewChild("Style");
node.SetAttribute("id", "default");
# create model element for placemark
model = placemark.NewChild("Model");
model.SetAttribute("id", "model");
# create altitude mode element for model
model.NewChild("altitudeMode", "absolute");
# create location element for model
node = model.NewChild("Location");
node.NewChild("longitude", sprintf("%.8f", originPt.x) );
node.NewChild("latitude", sprintf("%.8f", originPt.y) );
node.NewChild("altitude", sprintf("%.1f", originPt.z + zOffset) );
# create orientation element for model
node = model.NewChild("Orientation");
node.NewChild("heading", NumToStr(angleToNorth) );
node.NewChild("tilt", "0");
node.NewChild("roll", "0");
# create link element for model
node = model.NewChild("Link");
node.NewChild("href", "files/" + imgName$ + ".dae");
# create ResourceMap element for model
node = model.NewChild("ResourceMap");
node1 = node.NewChild("Alias");
node1.NewChild("targetHref", "files/" + imgName$ +".jpg");
node1.NewChild("sourceHref", "files/" + imgName$ +".jpg");
############## Write the KML structure to a file and add to KMZ
kmldoc.Write(kmlfile$);
kmz.AddFile(kmzFilepath, kmlFilepath, "doc.kml");
DeleteFile(kmlfile$);
if (tempfile3exists) then tempfile3.Close();
print("End.");
Download