import_srtm.sml

# import_srtm.sml
# Written By: Dan Glasser MI - 01Aug03
#
# Revised 16Sep03; no longer requires a format file.
#
# This script imports and sets the georeference for the 'headerless' *.hgt SRTM-3 files.
# These are raw elevation grids (DEMs) produced by the Shuttle Radar Topography Mission.
# Each *.hgt file is a grid with 1201 lines and 1201 columns sampled at 3 arc-seconds
# of latitude and longitude, covering one square degree. 

# Georeference information is derived from the parsed filenames, which have the form
# N00W050.hgt, where the first three characters give the latitude (integer degrees) and 
# the last four characters give the longitude of the lower left (southwest) corner of 
# the tile. The top row and right column of each tile overlap the adjacent tile. 

# Note: the *.hgt data is unprocessed and may have gaps or holes.
#
# For more information on SRTM data see: http://www.jpl.nasa.gov/srtm/
# ftp download available from: ftp://edcsgs9.cr.usgs.gov/pub/data/srtm/
#

clear();

class FILEPATH filepath;
class STRINGLIST filenames;
class MieUSERDEFINEDRASTER srtm;

# get the latitude value from the *.hgt filename
func getLatitude( string hgtfile$ )
{
	local numeric latvalue = StrToNum( GetToken( hgtfile$, "NSEW.hgt", 1, 1 ) );
	string hemisphere = GetToken( hgtfile$, "EW.hgt0123456789", 1, 1 );

	if ( hemisphere == "N" )		# if northern hemisphere
	{
		return latvalue;
	}
	else if ( hemisphere == "S" )	# if southern hemisphere
	{
		return -latvalue;
	}
	else
	{
		print( "Error parsing filename for latitude" );
		print( "Bad file was: " + hgtfile$ );
		return -9999;
	}
}

# get the longitude value from the *.hgt filename
func getLongitude( string hgtfile$ )
{
	local numeric longvalue = StrToNum( GetToken( hgtfile$, "NSEW.hgt", 2, 1 ) );
	string hemisphere = GetToken( hgtfile$, "NS.hgt0123456789", 1, 1 );

	if ( hemisphere == "E" )		# if eastern hemisphere
	{
		return longvalue;
	}
	else if ( hemisphere == "W" )	# if western hemisphere
	{
		return -longvalue;
	}
	else
	{
		print( "Error parsing filename for longitude" );
		print( "Bad file was: " + hgtfile$ );
		return -9999;
	}
}

# set the parameters for USERDEFINEDRASTER import
proc setImportParameters()
{
	srtm.NumLins = 1201;
	srtm.NumCols = 1201;
	srtm.DataType = "16-bit signed";
	srtm.ByteOrder = "High-Low";
	srtm.UseFile = 0;
	srtm.DoPyramid = 1;
	srtm.DoCompress = 1;
}

# get the SRTM data directory
string defaultpath$ = _context.ScriptDir;
filepath.SetName( GetDirectory( defaultpath$, "Please select the directory containing the SRTM files for import" ) );

setImportParameters();

print( filepath.GetPath() );

filenames = filepath.GetFileList( "*.hgt" );

numeric filecount = filenames.GetNumItems();
print( filecount, "files found" );

# get the output file name
string outputfile$ = GetOutputFileName(filepath.GetName(), "Enter a file name for the imported srtm raster rvc", ".rvc" );

numeric lat, long;

numeric i;
for i = 0 to filecount-1
{
	print( "Now importing: " + filenames.GetString( i ) + " raster number: " + NumToStr( i + 1 ) );

	lat = getLatitude( filenames.GetString( i ) );
	long = getLongitude( filenames.GetString( i ) );

	if ( lat != -9999 && long != -9999 )
	{
		string inputfile$ = filepath.GetPath() + "\\" + filenames.GetString( i );
		string objname$ = FileNameGetName( filenames.GetString( i ) );

		# need to import raster
		ImportRaster( srtm, inputfile$, outputfile$, objname$, "file imported from: " + filenames.GetString( i ) );

		# 'open' the raster so that it can be georeferenced
		Raster importedRaster;
		OpenRaster( importedRaster, outputfile$, objname$ );

		SetNull( importedRaster, -32768 );

		# set the raster extents (in object coordinates) based on centers of corner cells
		numeric xmin, ymin, xmax, ymax;
		xmin = 0.5;			# left edge
		xmax = 1200.5;		# right edge
		ymin = 0.5;			# top edge
		ymax = 1200.5;		# bottom edge

		# save the extents as the 'source' and the lat long values as 'dest' arrays
		# to be used to create control points at the raster corners
		# [1] = ul corner, [2] = ll corner, [3] = ur corner
		array numeric sourceX[3]; 
		array numeric sourceY[3];
		array numeric destX[3];
		array numeric destY[3];

		sourceX[1] = xmin;
		sourceY[1] = ymin;
		sourceX[2] = xmin;
		sourceY[2] = ymax;
		sourceX[3] = xmax;
		sourceY[3] = ymin;

		destX[1] = long;
		destY[1] = lat + 1;
		destX[2] = long;
		destY[2] = lat;
		destX[3] = long + 1;
		destY[3] = lat + 1;

		# set projection parameters
		class MAPPROJ mapproj;
		mapproj.SetSystemLatLon();
		mapproj.Datum = "World Geodetic System 1984";

		# use 'source' and 'dest' arrays to create control point georeference for the object
		CreateControlPointGeorefDefaultAccuracy( importedRaster, mapproj, 3, sourceX, sourceY, destX, destY );

		CloseRaster( importedRaster );
	}
	else
	{
		print( "Check file name and structure, lat/long not computed correctly" );
	}
}

# take care of plural
string s$ = "s";

if ( filecount == 1 )
{
	s$ = "";
}

# let user know that it's is done.  Enjoy!
print( "Done. ", filecount, "raster" + s$ + " imported." );