CRS.sml

# CRS.sml
# Sample script for tutorial Writing Scripts with SML.
# Coordinate Reference Systems and Coordinate Transformations.
# Transforms input latitude/longitude coordinates (in WGS84 Datum)
# to UTM coordinates in the appropriate zone 
clear();
numeric lat, lon;
class STRING prompt$;
prompt$ = "Enter a latitude value (-90 to +90 degrees)";
lat = PopupNum(prompt$, 0, -90, 90, 6);
prompt$ = "Enter a longitude value (-180 to +180 degrees)";
lon = PopupNum(prompt$, 0, -180, 180, 6);
printf("Entered: lat = %.6f, lon = %.6f\n", lat, lon); 
# set up class for input Coordinate Reference System: 
# WGS84 / Geographic; assign using its text ID string
class SR_COORDREFSYS latlonCRS = "Geographic2D_WGS84_Deg"; 
printf("Input CRS is: %s\n", latlonCRS.Name);
# compute the UTM zone number and longitude of its central
# meridian using the entered longitude
numeric zoneNum, centMerid;
zoneNum = ceil( (180 + lon) / 6 );
centMerid = (zoneNum * 6) - 183;
### set up a Coordinate Reference System for using this
### UTM zone and the WGS84 datum.  This requires a
### Coordinate System, Datum, and Coordinate Operation Definition (projection)
class SR_COORDREFSYS utmCRS;
# set projected 2D coordinate system (m) using its ID in the EPSG codespace
utmCRS.CoordSys = "EPSG:4400";
utmCRS.Datum = latlonCRS.Datum;  # get datum from the existing CRS
class SR_COORDOPDEF projectDef;	# projection definition class
# assign Transverse Mercator Projection (projection "method")
# using numeric ID in MicroImages (default) code space
projectDef.Method = "1909";
# 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 longitude of natural origin (central meridian of UTM zone) determined
# in calculation above
projectDef.SetParmValueNum("LongitudeOfNaturalOrigin", centMerid);
# assign false northing based on whether north or south of equator
if (lat < 0) then
	projectDef.SetParmValueNum("FalseNorthing", 10000);
else
	projectDef.SetParmValueNum("FalseNorthing", 0);
# set this projection for the UTM CRS
utmCRS.CoordOpDef = projectDef;
printf("Output CRS is: \n%s\n", utmCRS.GetDisplayStr());
# set up coordinate transformation from input to output CRS
class TRANS2D_MAPGEN transGeoToUTM;
transGeoToUTM.InputCoordRefSys = latlonCRS;
transGeoToUTM.OutputCoordRefSys = utmCRS;
 
class POINT2D pt;		# assign latitude and longitude as x and y of POINT2D
pt.x = lon;		pt.y = lat;
pt = transGeoToUTM.ConvertPoint2DFwd(pt);		# apply the coordinate tansformation
printf("Output UTM Zone %d coordinates:\n", zoneNum);
printf("Northing = %.6d m, Easting = %.6d m", pt.y, pt.x);