The NGCD software Geomag, http://www.ngdc.noaa.gov/ , computes the magnetic declinations. There are two models of the earth magnetic field. The first from Internation Association of Geomagnetism and Aeronomy (IAGA), the second for U.S. Dept. of Dedfense and NATO.

The program computes the vector components of the magnetic field and outputs also declination, inclination, and field strength (both total and in the horizontal plane). It is interactive. Upon start it asks which model to use; enter the model filename, eg, "IGRF10.unx". Next it asks the date. Then the coordinate system, either WGS84 (latitude, longitude and height) or spherical (earth centered), and the coordinates.

The output lists the values of declination, inclination, horizontal field strength, North, East and vertical components, and total field strength, and their rates of variation.

Digital elevation models (DEM) are useful to get the height grids for the surface command.

Height data for the world at a resolution of 3x3 degrees can be downloaded from http://www.panoramaview.com . This site has also some data at 1x1 degrees, which means roughly 30x30 m. Each file contains a square array of 1201x1201 values (3601x3601 for the 1x1 degree files). The coordinate of the point at the lower-left corner are in the filename; for example the 1x1 degree file N45E08.hgt contains the data from 45 to 46°North, and from 8 to 9°East. The height values are listed rowwise from the top (northest row) down (southest row). Each row contains the data from west to east.

To get the grid data useful for the Therion surface command, you can use gdal, as described in the Therion wiki [wollez].

Another program to get grid data is TerrainTool. It uses the data from SRTM (Shuttle Radar Topography Mission) and ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer). The SRTM data are dowmloaded automatically by the program from the NASA server. They are mostly at 3 arc-sec resolution and contains gaps. The ASTER data are at 1 arc-sec, and must be installe manually. TerrainTool handles several coordinate systems, has color display of the DEM, and saves the data in Survex and Therion formats.

The proj.4 library is used in Therion. It can be downloaded from http://trac.osgeo.org/proj/ .

The distribution contains also the executables, proj and cs2cs. For example,

• proj -le lists the ellipsoids
• proj -ld lists the datum
• proj -lp lists the projections
• proj -l=id lists info for the projection "id"

The coordinate transformation is interactive, taking input fron the stdin and writing the output to stdout (redirect if neceessary). On the command line you must specify the input and the output coordinate frames separated by +to. You can also add further command options ("-r" reverse the input x-y, "-s" reverse the output x-y).

The proj.4 coordinate system definitions let us specify several parameters,

• +proj=projection_name
• +zone=zone_number
• +datum=datum_name
• +ellps=ellipsoid_name
• +R=spherical_radius
• +a=major_axis
• +es=square_eccentricity
• +R_A states that the ellipsoid has the equivalent surface of a sphere
• +R_a states that the arithmetic mean of the ellipsoid axes is equal to the radius of the sphere
• +x_0=false_easting added to the x coordinate to have only positive values
• +y_0=false_northing
• +lon_0=central meridian
• +lat_0=central latitude
• +k=...
• +geoc states that the coordinates are geocentric

The coordinate can be specified also with +init=file:key where "file" is a file listing the coordinate system specifications, and "key" is the tag of the wanted coordinate system. Here are a few examples of coordinate specifications. The sample data refers to the same point; the first coordinate is the X (east), the second one is Y (north). The input can have also the Z coordinate (ellipsoid height). The output of the conversion contains a different height is the two coordinate systems have different ellipsoids.

• WGS84: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs. Data format example: 9d04'20.45"E 45d50'18.89"N.
• UTM32: +proj=utm +zone=32 +datum=WGS84. Data format example: 505618.3 5076115.5.
• Gauss-Boaga: +proj=tmerc +lat_0=0 +long_0=-3.45233333333333 +k=0.9996 +x_0=1500000 +y_0=0 +ellps=intl +pm=rome +units=m +no_defs +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68. For the east zone change +x_0=2520000 +lon_0=15. Data format example: 1505646.3 5076136.3.
• ED50: +proj=utm +zone=32 +ellps=intl +units=m +towgs84=-87,-98,-120,0,0,0,0 for data in Km, 505701.5 5076313.5, and +proj=longlat +ellps=intl +no_defs +towgs84=-87,-98,-120,0,0,0,0 for data in long-lat, 9d4'24.299"E 45d50'22.201"N.
• Rome40: +proj=longlat +ellps=intl +pm=rome +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68. Data format example: 3d22'46.7"W 45d50'16.5"N.

It is convenient to put the input data in an input data file, one coordinate pair per line. Example

  cs2cs +proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs +to 
        +proj=utm +zone=32 +datum=WGS84S < input_data_file

Some useful options:

• -I (inverse tranform): converts from the "+to" coordinate system
• -r (reserved): input coordinates are swapped, eg, y-x instead of x-y
• -s (reverse output): swapped output coordinate, eg, lat-long instead of long-lat
• -f format: specify output format, eg. DMS (degrees, minutes, seconds) or "%.6f" (decimal)

Another conversion program is geotrans, "W 45d50'16.5"N