Using The Dmapf Routines (FORTRAN)

The Dmapf routines are an extension of the Cmapf routines to a general orientation of projection surface to the Earth’s surface. The Cmapf routines are a library of FORTRAN-language subroutines and functions which process geographic and dynamic data between grids and geographical coordinates on a family of conformal map projections.

The Cmapf routines deal only with the cases in which the projection surfaces are either axially aligned with the Earth’s North-South axis, (in the case of the Mercator Cylindrical or the Lambert Conic Projection) or orthogonal to the Earth’s axis (in the case of the Polar Stereographic Projection.

The Dmapf routines extend the Cmapf routines by allowing projection surfaces whose axes lie in the equatorial plane (transverse Mercator and Transverse Lambert) or in other directions (Oblique Stereographic, orthogonal to a diameter through any specified point on the Earth’s surface). Oblique Mercator and Lambert projections are provided for, but rarely used.

Data characterizing the map projections and the grids on them is embedded in maparam structures, which are REAL*4 arrays of dimension 15, rather than 9 as in the Cmapf case. 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.

Initializing (configuring) the maparam structure:

Stage I – Selecting a Map Projection.
To initialize a map parameter (maparam) structure with a Conformal Polar Projection, use any one of the following subroutines:

  • Normal (aligned with Earth’s axis) Projections
    • dimension stcprm(15)
    • real reflat, reflon, tnglat
    • call stlmbr(stcprm, reflat, reflon)
      or
      equivalently
    • call stcmap(stcprm, tnglat, reflon)

    for a normal (polar) Lambert projection with reflat being the latitude at which the projection plane, cone, or cylinder is tangent to the Earth, while reflon is the longitude of some point 180 degrees away from the cut. This operates exactly the same as in the Cmapf routines.

  • Transverse (aligned with an equatorial diameter)
    Projections

    • dimension stcprm(15)
    • real r_lat, r_long
    • call stvmrc(stcprm, r_lat, r_long)

    will define a transverse Mercator projection onto a cylinder tangent to the Earth along a meridian. The meridian passes through the reference point r_lat, r_long, and has the least distortion near that meridian. The cylinder is cut half way from the reference point, and the cut is the projection of a great semi-circle, 180 degrees away from the r_lat latitude, extending to equatorial points 90 degrees removed from the r_long meridian.

  • Oblique (aligned with any other diameter) Projections:
    • dimension stcprm(15)
    • real r_lat, r_lon, an_lat, an_lon, a_lat1, a_lng1, a_lat2, a_lng2
    • call sobstr(stcprm, r_lat, r_lon)
    • call sobmrc(stcprm, r_lat, r_lon, an_lat, an_lon)
      or
    • call soblmb(stcprm, r_lat, r_lon, a_lat1, a_lng1, a_lat2, a_lng2)

    will define, respectively, the oblique Stereographic, oblique Mercator, and oblique Lambert Conformal projections. In each case, r_lat and r_lon are the latitude and longitude of the reference point, which is the point on the Earth at which this map’s distortion is least (because the gradient of the scale is zero) and which is furthest from any cut.

    an_lat, an_lon, a_lat1, a_lng1, a_lat2, and a_lng2 are latitudes and longitudes of anchor points, which inform the routines of the location of the tangent circles for the Mercator and Lambert projections, respectively. The reference point, together with one or two anchor points, define the circle at which the projection surface is tangent to the Earth.

    In the Stereographic case, the one point r_lat,r_long is the point where the projection plane is tangent to the Earth.

    In the Mercator case, two points define a great circle, at which the projection cylinder is tangent to the Earth. The reference and anchor points must be two different points.

    In the Lambert Conic Conformal case, three points are required to define a small circle at which the projection cone is tangent to the Earth’s surface. The reference point and the anchor point(s) must be three distinct points.

The above routines allow for an extensive selection of projection types, but in practice, most projections will consist of the normal projections, the transverse Mercator, and the Oblique Stereographic.

Whichever of the above routines is selected, it will fill the structure “stcprm” with default data specifying the selected projection on the appropriate surface tangent to the globe at the indicated point or circle. If a cut is needed to lay the cone or cylinder out flat on the plane, it will be made opposite to (the maximum distance from) the reference point.

Stage II – Laying out a grid on a selected map projection.
This stage configures the maparam structure with the information to locate the points of the grid on the selected map projection selected in Stage I above. Two alternative means are provided: Either One-Point specification or Two-Point specification.

  1. Two-Point Specification:
    • dimension stcprm(15)
    • real x1, y1, xlat1, xlon1, x2, y2, xlat2, xlon2
    • call stcm2p(stcprm, x1, y1, xlat1, xlon1, x2, y2, xlat2, xlon2)

    causes stcprm to be configured so that the grid point with grid coordinates x1, y1 will be located on the map at latitude xlat1, longitude xlon1, while x2, y2 is located on the map at xlat2, xlon2. The two points must be different. This specification locates all other grid points, assumed arranged on the projection with equal spacing in the x- and y- directions.

    The two points will often be opposite corners of the intended grid, but many times a different selection will be preferred. For example, in a mercator projection, it will often be better to select 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:
    • dimension stcprm(15)
    • real x1, y1, xlat1, xlon1, scalat, scalon, gsize, orient
    • call stcm1p(stcprm, x1, y1, xlat1, xlon1, scalat, scalon, gsize, orient)

    causes stcprm to be configured so that the grid point with grid coordinates x1, y1 will be located on the map at latitude xlat1, longitude xlon1. The remaining configuration information 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 in both the x- and y- directions. Finally, at the scaling point, the y-lines of the grid make an angle of “orient” degrees clockwise (East) of the scalon meridian.

    The x1, y1 coordinates can be thought of as “pinning” that grid point to the map at xlat1,xlon1, 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 y-axis through the scaling point has the desired bearing.

    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:

Coordinate transforms:
  • dimension stcprm(15)
  • real xlat,xlon, x,y
  • call cll2xy(stcprm, xlat, xlon, x, y)

will return grid coordinates x, y according to the grid configured as stcprm, for the point whose latitude and longitude are xlat,xlon. 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.

  • dimension stcprm(15)
  • real xlat, xlon, x, y
  • call cxy2ll(stcprm, x, y, xlat, xlon)

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, conforming to the WMO standard for reporting winds there. (In an earlier version of CMAPF, the value of reflon was returned at the poles.)

Wind Vector transforms:
  • dimension stcprm(15)
  • real xlat, xlon, ue, vn, ug, vg
  • call cc2gll(stcprm, xlat, xlon, ue, vn, ug, vg)
    and
  • call cg2cll(stcprm, xlat, xlon, 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 xlat,xlon.

At the poles (xlat=90. or xlat= -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 xlon meridian. Thus, if xlat,xlon = 90.,-70., “North” will be assigned to the direction of the 110E 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 xlat=90.,xlon=180. for the North Pole, and xlat=-90.,xlon=0. for the South Pole.

  • call cw2gll(stcprm, xlat, xlon, ue, vn, ug, vg)
  • call cg2wll(stcprm, xlat, xlon, 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 degree (111km or 60nm) or 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 routines, above.

Gridsize Evaluation:
  • dimension stcprm(15)
  • real gsize, xlat, xlon, x, y
  • gsize = cgszll(stcprm, xlat, xlon)
    and
  • gsize = cgszxy(stcprm, x, y)

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

Curvature Vector:
  • dimension stcprm(15)
  • real x, y, gx, gy, xlat, xlon
  • call ccrvxy (stcprm, x, y, gx, gy)
    and
  • call ccrvll (stcprm, xlat, xlon, gx, gy)

return, in gx and gy, the logarithmic gradient of the scale. This vector represents the curvature on the Earth of a straight line segment on the map. Units are in radians per kilometer on the Earth’s surface. See the paper for details on usage in equations.

North Polar Vector evaluations:
  • dimension stcprm(15)
  • real x, y, xlat, xlon, enx, eny, enz
  • call cpolxy (stcprm, x, y, enx, eny, enz)
    and
  • call cpolll (stcprm, xlat, xlon, enx, eny, enz)

return, in enx, eny, and enz the components of the direction of the Earth’s axis. Components enx and eny provide the local direction of North in grid coordinates, while enz is proportional to the local Coriolis parameter.

Other Routines:
In addition to the functions defined above for the end user, there are also a number of functions that assist the end-user routines, but are not generally used by the end user. These functions include

subroutines:
mpstrt, basm2g, basg2m, proj_3d, geo_ll, ll_geo,
xy_map, map_xy, map_xe, xe_map, xe_xy, xy_xe, cgrnll
and cgrnxy
real functions:
x_prod, cml2xy, csabva, snabva, lnabva, xpabva
and atnabv

These subroutine and function names should not be used in code created by the end user.

[/vc_column_text][/vc_column][/vc_row]