Using The EMAPF Functions (FORTRAN)

The EMAPF functions are an extension of the CMAPF functions to accommodate Ellipsoidal Geoids. The CMAPF functions are a library of FORTRAN-callable routines which process geographic and dynamic data between grids and geographical coordinates on a family of conformal map projections. The CMAPF functions assume a spherical shape to the Earth, and map this shape onto Lambert Conformal, Polar Stereographic or Mercator Projections. The EMAPF functions do the same to an ellipsoidal shape for the Earth, which may be chosen from a collection of geoids that have been used for the Earth, including WGS84 currently used by the United States, or the user may specify his own ellipsoidal shape.

The figure of the Earth (its “geoid”) is generally considered to be an ellipsoid, formed by rotating an ellipse, whose minor axis lies along the Earth’s axis, around that same axis. The equator is then a circle, whose radius ARAD is half the major axis of the ellipse, and the half-length of the minor axis, BRAD, may be thought of as the polar radius.

Note: in all cases, North latitudes will be given as positive degrees, South latitudes negative degrees, East longitudes positive degrees and West longitudes negative degrees.

Parameters determining the shape of the geoid, and the scale and orientation of the projection itself, are set in a maparam structure in the form of a dimension 14 array

      DIMENSION STCPRM(14)

by means of the following three-stage initialization procedure:

Initializing (configuring) the maparam structure:

First Stage – Specifying a Geoid.
NOTE: It is VITAL to FIRST select a geoid. If other routines are given maparam structures with inconsistent geoids, they will return error codes and invalid results.

We can specify the geoid in either one of two ways:

  1. Predefined Geoids:
    We have provided a hardwired standard set of Geoids. A listing of the Geoids can be obtained by using the command “listGeoid”, which is a part of the EMAPF package. It will list on stdout all the geoids in the EMAPF package by name and by EMAPF Index. To specify in your program which geoid to use,

       CALL GEOIDIX(STCPRM, NGEOID)

    will place in STCPRM the dimensions and eccentricity of the geoid whose index in the listGeoid listing is NGEOID. 0 is the index for WGS84, i.e. the standard geoid established as the World Geodetic System in 1984, currently the international standard geoid. Other standard geoids are (8) the International, or HAYFORD, adopted internationally from 1924 to 1967, and (12) the IUGG of 1967, adopted from 1967 to 1984.

  2. User-defined Geoids:
    Alternatively, a user can define his own geoid by providing these two items:

    1. ARAD – Radius of the Equator
    2. BRAD – Radius from Earth center to Pole
          CALL SETGEOID(STCPRM, ARAD, BRAD)

    will place the dimensions and eccentricity for a geoid of the indicated size in STCPRM. Note: Often, published geoid parameters will be ARAD and flattening, the flattening being specified in the form 1.0/###.####. Then BRAD = ARAD * (1. – FLATTEN). Thus:

         ARAD = 6378.1234
     FLATTEN = 1./298.2
     CALL SETGEOID(STCPRM, ARAD, ARAD*(1. - FLATTEN))

    To specify a sphere of radius R, use

         CALL SETGEOID(STCPRM, R, R)
Second Stage – Selecting a Map Projection
After selecting a geoid, a map projection may be selected. To initialize a maparam structure with a Lambert Conformal Projection:

      DIMENSION STCPRM(9)

      CALL GEOIDIX(STCPRM, 0)
      TNGLAT = 40.
      REFLON = -75.
								
      CALL STLMBR(STCPRM, TNGLAT, REFLON)

will initialize STCPRM to a Lambert Conformal Projection on a cone which is tangent at latitude 40°N to a WGS 84 Geoid and is centered at longitude 75° West.

Types of projection are determined by the value of TNGLAT, which ranges from -90° to +90°. If TNGLAT = 0, the projection is a Mercator Projection. If TNGLAT = 90 or ‑90, the projection is a Polar Stereographic Projection, centered at the North or South Pole, respectively. For other values of TNGLAT, the projection is a Lambert Conic Projection which is tangent to the geoid at latitude TNGLAT. TNGLAT also defines the width of the cone: it is the angle between the cone and its axis.

TNGLAT also determines the spread of lines of longitude. The spread of such lines is proportional to the sine of TNGLAT. Thus, two longitude lines, say 90°W and 100°W, may be 10° apart at the poles on the globe, but they are drawn 7.071° apart on the map (assuming, e.g. that TNGLAT = 45°).

Most published maps using the Lambert Conformal Projection do not specify TNGLAT. Rather, they specify a nominal scale and two reference latitudes. Contrary to popular misconception, these latitudes do not specify circles on the geoid through which the projection cone can be made to pass. Instead, these are the two latitudes at which the actual scale on the map matches the nominal scale, and consequently are two latitudes where the scale at one is the same as the scale at the other. For each pair of latitudes, there is only one projection cone which allows this to be true, and that cone is slightly wider (TNGLAT has greater magnitude) than the cone that passes through the two latitude circles.

The FUNCTION EQVLAT is provided to find the tangent latitude corresponding to the cone for which the scale at each reference latitude is equal to the scale at the other. We have not provided a similar function to determine the alternate cone that passes through the two latitude circles themselves. To compute the TNGLAT corresponding to two reference latitudes and the geoid configured in STCPRM, use:

       TNGLAT = EQVLAT(STCPRM, REFLAT1, REFLAT2)
       CALL STLMBR(STCPRM, TNGLAT, REFLON)

STCPRM is then configured to the projection cone appropriate to the two reference latitudes. It is important to have already partially configured STCPRM by selecting a geoid in Stage 1, as the actual value of EQVLAT depends on the flattening of the geoid.

Special cases – In the event REFLAT1 = REFLAT2, EQVLAT returns their common value, corresponding to a cone tangent at that latitude circle. If REFLAT1 = -REFLAT2, EQVLAT returns 0, corresponding to a cylinder, i.e. to a Mercator projection. If either REFLAT1 or REFLAT2 is +90. (or -90.), EQVLAT returns +90. (or -90.), corresponding to a polar stereographic projection tangent at the appropriate pole.

Third Stage – Laying out a grid on a selected map projection.
This stage configures the maparam structure with the information to locate the points of a grid on the selected map projection as unrolled onto a flat surface. This allows the EMAPF routines to relate the X,Y (or (I,J)) coordinates of the grid with the latitude-longitude coordinates of the Earth. Note that the first three configuration stages must be performed before this final stage.

Two alternative means are provided: Either

One-Point specification or Two-Point specification.

  1. Two-Point Specification: 
           CALL STCM2P(STCPRM, X1, Y1, LAT1, LON1, X2, Y2, LAT2, LON2)

    configures STCPRM so that the grid point with grid coordinates X1, Y1 will be located on the map at latitude LAT1, longitude LON1, while the grid point X2, Y2 is located on the map at LAT2, LON2. The two points must be different.

    This specification is sufficient to locate all other grid points, which are assumed to be arranged on the projection with equal spacing in the X- and Y- directions.

    The two points will often be chosen as opposite corners of the intended grid with a view to minimize rounding errors in the rest of the grid. But often, a different choice of grid points will be preferred. For example, in a Mercator projection, it is better to choose two points of the same latitude or the same longitude, to ensure exact N-S orientation of the grid. Alternatively, several subgrids of an overall grid may be guaranteed compatibility by using the same two points of the overall grid, even if neither of the points lie in the subgrid at all.

  2. One-Point Specification: 
           CALL STCM1P(STCPRM, X1, Y1, LAT1, LON1, SCALAT, SCALON, GSIZE, ORIENT)

    configures STCPRM so that the grid point with grid coordinates X1, Y1 will be located on the map at latitude LAT1, longitude LON1. The information needed to space and orient all other points on the grid is provided at the scaling point whose latitude and longitude are SCALAT, SCALON. At that point, and incidentally every point on the latitude circle SCALAT, the size of a grid cell is GSIZE kilometers in both the X- and Y- directions, from which the subroutines can compute the grid spacing at all other latitudes. Finally, to orient the other grid points on the map projection, the Y-lines of the grid make an angle of “ORIENT” degrees clockwise (East) of the SCALON meridian at the scaling point. From this, the orientation of all other grid lines at all other points can be calculated, fully determining the grid on the projection.

    The X1,Y1 coordinates can be thought of as “pinning” that grid point to the map at LAT1,LON1, the GSIZE value can be thought of as “zooming” the grid in and out until the desired size is attained, and the “ORIENT” specification can be thought of as rotating the grid about the pinned point until the grid has the wanted orientation.

    Hint:   frequently, a grid will be published as having an “LOV” value of a longitude intended to be vertical. For such grids, SCALON = LOV and ORIENT = 0.

Using the configured maparams:

Once STCPRM has been configured to the desired geoid, projection, and grid layout, the following routines will perform two-way transformations of coordinates and wind vectors, and illuminate useful geometry for modelers:

Coordinate transforms:
       CALL CLL2XY(STCPRM, LAT, LON, X, Y)

will return grid coordinates X,Y according to the grid configured as STCPRM, for the point whose latitude and longitude are LAT,LON. Longitude values are adjusted to fit between REFLON-180. and REFLON+180., so that nearby points on the globe remain nearby points on the map, unless they straddle the “cut” at longitude REFLON+180., and find themselves on opposite edges of an unrolled cylinder or cone.

       CALL CXY2LL(STCPRM, X,Y, LAT, LON)

will return latitude and longitude parameters according to the grid configured as STCPRM, for the point whose grid coordinates are X,Y. Longitude values will be returned in the range from -180. to +180.

If the point X, Y corresponds to the North or South Pole, where longitude is undefined, the returned longitude will be 180. for the North Pole, 0. for the South Pole, for compatibility with the WMO standard for reporting winds there. (In an earlier version of CMAPF, the value of REFLON was returned here.)

Wind Vector transforms:
     CALL CC2GLL(STCPRM, LAT, LON, UE, VN, UG, VG)			 
     and 			
     CALL CG2CLL(STCPRM, LAT, LON, UG, VG, UE, VN)

will convert the East- and North- components UE,VN (Cartographic components) of the wind to and from components UG, VG of the wind in the X- and Y- directions of the grid (Grid components). The first transfers from cartographic components to grid components, the second does the reverse. In both cases, the wind is that observed at the specified cartographic coordinates LAT, LON.

At the poles (LAT = 90. or LAT = -90.), the direction of “North”, which determines the meaning of UE and VN, is based on the nominal value of longitude. The direction of North will agree with that perceived by an observer stationed a short distance away on the LON meridian. Thus, if LAT,LON = 90.,-70., “North” will be assigned to the direction of the 110°E meridian.

     CALL CC2GXY(STCPRM, X,Y, UE,VN, UG,VG) 
     and
     CALL CG2CXY(STCPRM, X,Y, UG,VG, UE,VN)

will convert the East- and North- components UE, VN (Cartographic components) of the wind to and from components UG, VG of the wind in the X- and Y- directions of the grid (Grid components). The first transfers from cartographic components to grid components, the second does the reverse. In both cases, the wind is that observed at the specified grid coordinates X,Y.

If X,Y refers to the North or to the South pole, the orientation of “North” is defined according to World Meteorological Organization (WMO) code 878. The direction of “North” is that of the Greenwich meridian, which is the equivalent of using LAT = 90.,LON = 180. for the North Pole, and LAT = -90.,LON = 0. for the South Pole.

     CALL CW2GLL(STCPRM, LAT, LON, UE, VN, UG, VG)
     CALL CG2WLL(STCPRM, LAT, LON, UG, VG, UE, VN) 
     CALL CW2GXY(STCPRM, X, Y, UE, VN, UG, VG)
       and 
     CALL CG2WXY(STCPRM, X, Y, UG, VG, UE, VN)

perform the same services as the above transforms, except that within 1° (111km or 60nm) of either pole, UE and VN are interpreted according to WMO code 878 above. In an earlier version of CMAPF, this was the definition of the C2G and G2C functions, above.

Gridsize Evaluation:
The scale of the projection varies from one latitude to another, and therefore, the size of each cell of the grid laid out on the projection changes from one point to another. The functions

       GSIZE = CGSZLL(STCPRM, LAT,LON) 
       and 
       GSIZE = CGSZXY(STCPRM, X,Y)

return the size, in kilometers, of a grid cell at the indicated point.

Curvature Vector:
The subroutines

       CALL CCRVXY(STCPRM, X,Y, GX,GY) 
       CALL CCRVLL (STCPRM, LAT,LON, GX,GY)

return, in GX and GY, the logarithmic gradient of the scale.

Units are in radians per kilometer on the Earth’s surface. If a path is drawn on the map that appears to be a straight line, it represents a path on the Earth that may, in fact, be curved. The component of this vector, perpendicular to that path, represents the actual curvature of the path on the Earth. See the paper for details on this vector’s usage in equations.

North Direction Vector and Coriolis evaluations:
       CALL CPOLXY (STCPRM, X, Y, ENX,E NY, ENZ)
and 
       CALL CPOLLL (STCPRM, LAT, LON, ENX, ENY, ENZ)

return, in ENX, ENY, and ENZ the components of a unit vector in the direction of the Earth’s axis. Components ENX and ENY provide the local direction of North aligned with grid coordinates, while ENZ is proportional to the local Coriolis parameter.