   Chapter 13: Geographical Projections and Plotting Maps

This chapter presents different  methods  to project geogra-
phical coordinates  onto a plane surface.  Several base maps
are stored in DISLIN for plotting maps.

13.1 Axis Systems and Secondary Axes

                         G R A F M P

The routine GRAFMP plots a geographical axis system.

The call is:  CALL GRAFMP (XA, XE, XOR, XSTP,
                           YA, YE, YOR, YSTP)        level 2

XA, XE        are the lower and upper limits of the X-axis.
XOR, XSTP     are the  first  X-axis label  and the step be-
              tween labels.
YA, YE        are the lower and upper limits of the Y-axis.
YOR, YSTP     are the  first  Y-axis label  and the step be-
              tween labels.

Notes:     -  GRAFMP  must be called  from level 1  and sets
              the level to 2.
           -  The axes must be linear  and  scaled in ascen-
              ding order. In general,  X-axes must be scaled
              between  -180  and 180 degrees and  Y-axes be-
              tween -90 and 90 degrees.
           -  For elliptical projections, the plotting of an
              axis system will be suppressed. This will also
              be done for azimuthal projections with YE - YA
              > 90.
           -  The statement  CALL GRIDMP (I, J)  overlays an
              axis system with a longitude and latitude grid
              where I and J are the number of grid lines be-
              tween labels in the X- and Y-direction.

                         X A X M A P

The routine XAXMAP plots a secondary X-axis.

The call is:  CALL XAXMAP (A, B, OR, STEP, CSTR, NT, NY)
                                                     level 2

A, B          are the lower and upper limits of the X-axis.
OR, STEP      are the  first label  and the step between la-
              bels.
CSTR          is a  character  string  containing  the  axis
              name.
NT            indicates how ticks,  labels and the axis name
              are plotted. If NT = 0,  they are plotted in a
              clockwise direction. If NT = 1, they are plot-
              ted in a counter-clockwise direction.
NY            defines the horizontal position of the X-axis.
              A secondary axis  must be  located inside  the
              axis system.

                        Y A X M A P

The routine YAXMAP plots a secondary Y-axis.

The call is:  CALL YAXMAP (A, B, OR, STEP, CSTR, NT, NX)
                                                     level 2

A, B          are the lower and upper limits of the Y-axis.
OR, STEP      are the  first label  and the step between la-
              bels.
CSTR          is a  character  string  containing  the  axis
              name.
NT            indicates how ticks,  labels and the axis name
              are plotted. If NT = 0,  they are plotted in a
              clockwise direction. If NT = 1, they are plot-
              ted in a counter-clockwise direction.
NX            defines the  vertical  position of the Y-axis.
              A secondary axis  must be  located inside  the
              axis system.

13.2 Defining the Projection

Since a  globe cannot be unfolded into a plane surface, many
different  methods have been  developed to represent a globe
on a plane surface.  In cartography,  there are  4 basic me-
thods differentiated by attributes  such as  equal distance,
area and angle.

The 4 basic methods are:

a) Cylindrical Projections

   The surface  of the globe  is projected  onto  a cylinder
   which can be unfolded  into a  plane surface  and touches
   the globe at the equator. The latitudes and longitudes of
   the globe are projected as straight lines.

b) Conical Projections

   The surface  of the globe is projected  onto a cone which
   can  also be  unfolded  into  a plane  surface.  The cone
   touches or intersects the  globe  at two  latitudes.  The
   longitudes are projected  as straight lines  intersecting
   at the top of the cone and the latitudes are projected as
   concentric circles around the top of the cone.

c) Azimuthal Projections

   For azimuthal projections, a hemisphere is projected onto
   a plane which touches  the hemisphere  at a point  called
   the map pole.  The longitudes and latitudes are projected
   as circles.

d) Elliptical Projections

   Elliptical projections  project the entire surface of the
   globe onto an elliptical region.

                       P R O J C T

The routine PROJCT selects the geographical projection.

The call is:  CALL PROJCT (CTYPE)                    level 1

CTYPE         is a character string defining the projection.
  = 'CYLI'    defines a cylindrical equidistant projection.
  = 'MERC'    selects the Mercator projection.
  = 'EQUA'    defines a cylindrical equal-area projection.
  = 'HAMM'    defines the elliptical projection of Hammer.
  = 'AITO'    defines the elliptical projection of Aitoff.
  = 'WINK'    defines the elliptical projection of Winkel.
  = 'SANS'    defines the elliptical Mercator-Sanson projec-
              tion.
  = 'CONI'    defines a conical equidistant projection.
  = 'ALBE'    defines a conical equal-area projection
              (Albert).
  = 'CONF'    defines a conical conformal projection.
  = 'AZIM'    defines an azimuthal equidistant projection.
  = 'LAMB'    defines an azimuthal equal-area projection.
  = 'STER'    defines an azimuthal stereographic projection.
  = 'ORTH'    defines an azimuthal orthographic projection.
                                    Default: CTYPE = 'CYLI'.

Notes:     -  For cylindrical equidistant  projections,  the
              scaling  parameters in  GRAFMP  must be in the
              range:

                         -540 <= XA <= XE <= 540
                         -180 <= YA <= YE <= 180

              For Mercator projections:

                         -540 <= XA <= XE <= 540
                         - 85 <= YA <= YE <=  85

              For cylindrical equal-area projections:

                         -540 <= XA <= XE <= 540
                         - 90 <= YA <= YE <=  90

              For elliptical projections:

                         -180 <= XA <= XE <= 180
                         - 90 <= YA <= YE <=  90

              For conical projections:

                         -180 <= XA <= XE <= 180
                            0 <= YA <= YE <=  90  or
                         - 90 <= YA <= YE <=   0

              For azimuthal projections  with  YE - YA > 90,
              the  hemisphere  around  the map pole  is pro-
              jected  onto  a  circle. Therefore,  the hemi-
              sphere  can be selected with the map pole. The
              plotting  of the  axis system will be suppres-
              sed.  If YE - YA <= 90,  the part of the globe
              defined by the axis scaling  is projected onto
              a rectangle. The map pole will be set by GRAF-
              MP to ((XA+XE)/2, (YE+YA)/2).  The scaling pa-
              rameters must be in the range:

                         -180 <= XA <= XE <= 180 and
                                  XE - XA <= 180
                         - 90 <= YA <= YE <=  90
           -  For all projections except the default projec-
              tion, longitude and latitude coordinates  will
              be projected with the same scaling factor  for
              the  X- and  Y-axis. The scaling factor is de-
              termined  by the  scaling  of the Y-axis while
              the scaling of the  X-axis  is used  to centre
              the map.  The longitude  (XA+XE)/2 always lies
              in the centre of the axis system.
           -  Geographical projections  can only be used for
              routines described in this chapter or routines
              that plot contours.

13.3 Plotting Maps

                         W O R L D

The routine WORLD plots coastlines and lakes.

The call is:  CALL WORLD                             level 2

Note:         The routine WORLD  supports also some external
              map coordinates (see MAPBAS).

                         S H D M A P

The routine SHDMAP plots shaded continents.

The call is:  CALL SHDMAP (CMAP)                     level 2
CMAP          is a character string defining the continent.
  = 'AFRI'    means Africa.
  = 'ANTA'    means the Antarctic.
  = 'AUST'    means Australia.
  = 'EURA'    means Europe and Asia.
  = 'NORT'    means North America.
  = 'SOUT'    means South America.
  = 'LAKE'    means lakes.
  = 'ALL'     means all continents and lakes.

Note:         Shading patterns  can be selected with  SHDPAT
              and MYPAT.  Colours can be defined with  COLOR
              and SETCLR.

                        S H D E U R

The routine SHDEUR plots shaded European countries.

The call is:  CALL SHDEUR (INRAY, IPRAY, ICRAY, N)   level 2

INRAY         is an  integer array  containing the countries
              to be shaded. INRAY can have the following va-
              lues:

           1: Albania     13: Iceland        24: Portugal
           2: Andorra     14: Italy          25: Romania
           3: Belgium     15: Yugoslavia     26: Sweden
           4: Bulgaria    16: Liechtenstein  27: Switzerland
           5: Germany     17: Luxembourg     28: Spain
           6: Denmark     18: Malta          29: CSFR
           8: England/GB  19: Netherlands    30: Turkey
           9: Finland     20: North Ireland  31: USSR
          10: France      21: Norway         32: Hungary
          11: Greece      22: Austria
          12: Ireland     23: Poland

IPRAY         is an  integer array  containing  shading pat-
              terns.
ICRAY         is an integer array containing colour numbers.
N             is the number of countries to be shaded.

Note:         The  plotting  of outlines  can be  suppressed
              with CALL NOARLN.

13.4 Plotting Data Points

                         C U R V M P

The routine CURVMP plots curves through data points or marks
them with symbols.

The call is:  CALL CURVMP (XRAY, YRAY, N)            level 2

XRAY, YRAY    are real arrays containing the data points.
N             is the number of data points.

Notes:     -  CURVMP is similar to  CURVE except that only a
              linear interpolation can be used.
           -  In general,  a line between two points  on the
              globe  will  not  be  projected  as a straight
              line. Therefore, CURVMP interpolates lines li-
              nearly by small steps.  An alternate  plotting
              mode can be set with MAPMOD.

13.5 Parameter Setting Routines

                         M A P B A S

The routine  MAPBAS  defines the map data file  used  in the
routine WORLD. A DISLIN map file and some external map files
compiled by Paul Wessel can be used.  The external map files
can be copied via FTP anonymous from the servers

   ftp://ftp.ngdc.noaa.gov/MGG/shorelines/
   ftp://kiawe.soest.hawaii.edu/pub/wessel/gshhs/

The external map files 'gshhs_l.b', 'gshhs_i.b', 'gshhs_h.b'
and  'gshhs_f.b'  must be copied  to the map subdirectory of 
the DISLIN directory.

The call is:  CALL MAPBAS (CBAS)                  level 1, 2

CBAS          is a character string defining the map data
              file:

 = 'DISLIN'   defines the DISLIN base map.
 = 'GSHL'     defines 'gshhs_l.b' as base map.
 = 'GSHI'     defines 'gshhs_i.b' as base map.
 = 'GSHH'     defines 'gshhs_h.b' as base map.
 = 'GSHF'     defines 'gshhs_f.b' as base map.
                                  Default: CBAS = 'DISLIN'.

                         M A P L E V

The routine  MAPLEV  defines  land or lake  coordinates  for 
WORLD if the external map files from Paul Wessel are used.

The call is:  CALL MAPLEV (COPT)                  level 1, 2

COPT          is a character string that can have the values
              'BOTH', 'LAND' and 'LAKE'.
                                    Default:  COPT = 'BOTH'.

                         M A P P O L

MAPPOL defines the map pole used for azimuthal projections.

The call is:  CALL MAPPOL (XPOL, YPOL)               level 1

XPOL, YPOL    are the longitude and latitude coordinates  in
              degrees where:
                         -180 <= XPOL <= 180 and
                          -90 <= YPOL <=  90.
                                           Default: (0., 0.)

Note:         For an azimuthal  projection  with  YE - YA <=
              90,  the  map pole  will be  set by  GRAFMP to
              ((XA+XE)/2, (YA+YE)/2).

                         M A P S P H

For an azimuthal projection with YE - YA > 90,  DISLIN auto-
matically projects a hemisphere  around the map pole  onto a
circle. The hemisphere can be reduced using MAPSPH.

The call is:  CALL MAPSPH (XRAD)                     level 1

XRAD          defines the  region  around  the map pole that
              will be projected onto a circle   (0 < XRAD <=
              90).
                                         Default: XRAD = 90.

                         M A P R E F

The routine MAPREF defines, for conical projections, two la-
titudes where the cone intersects or touches the globe.

The call is:  CALL MAPREF (YLW, YUP)                 level 1

YLW, YUP      are the lower and upper latitudes where:

                      0 <= YLW <= YUP <= 90  or
                   - 90 <= YLW <= YUP <=  0
                       Default:  YLW = YA + 0.25 * (YE - YA)
                                 YUP = YA + 0.75 * (YE - YA)

Note:         YLW and  YUP can have identical values and lie
              outside of the axis scaling.

                         M A P M O D

The routine MAPMOD determines how data  points will be  con-
nected by CURVMP.

The call is:  CALL MAPMOD (CMODE)                 level 1, 2

CMODE         is a character string defining  the connection
              mode.
 = 'STRAIGHT' defines straight lines.
 = 'INTER'    means  that lines will be interpolated linear-
              ly.
                                  Default:  CMODE = 'INTER'.

13.6 Conversion of Coordinates

                         P O S 2 P T

The routine  POS2PT converts map coordinates to plot coordi-
nates.

The call is:  CALL POS2PT (XLONG, YLAT, XP, YP)      level 2

XLONG, YLAT   are the map coordinates in degrees.
XP, YP        are the plot coordinates calculated by POS2PT.

The corresponding functions are:

              XP = X2DPOS (XLONG, YLAT)
              YP = Y2DPOS (XLONG, YLAT)

13.7 Examples

            PROGRAM EX13_1

            CALL SETPAG('DA4L')
            CALL DISINI
            CALL PAGERA
            CALL COMPLX

            CALL FRAME(3)
            CALL AXSPOS(400,1850)
            CALL AXSLEN(2400,1400)

            CALL NAME('Longitude','X')
            CALL NAME('Latitude','Y')
            CALL TITLIN('World Coastlines and Lakes',3)

            CALL LABELS('MAP','XY')
            CALL GRAFMP(-180.,180.,-180.,90.,
           *             -90., 90., -90.,30.)

            CALL GRIDMP(1,1)
            CALL WORLD

            CALL HEIGHT(50)
            CALL TITLE
            CALL DISFIN
            END

            PROGRAM EX13_2
            CHARACTER*6 CPROJ(3),CTIT*60
            DATA CPROJ/'Sanson','Winkel','Hammer'/

            CALL SETPAG('DA4P')
            CALL DISINI
            CALL PAGERA
            CALL COMPLX

            CALL HEIGHT(40)
            CALL AXSLEN(1600,750)

            NYA=3850
            DO I=1,3
              NYA=NYA-950
              CALL AXSPOS(250,NYA)

              CALL PROJCT(CPROJ(I))
              CALL NOCLIP
              CALL GRAFMP(-180.,180.,-180.,30.,
          *                -90., 90., -90.,15.)

              WRITE(CTIT,'(2A)') 
          *     'Elliptical Projection of ', CPROJ(I)
              CALL TITLIN(CTIT,4)
              CALL TITLE

              CALL WORLD
              CALL GRIDMP(1,1)
              CALL ENDGRF
            END DO

            CALL DISFIN
            END

            PROGRAM EX13_3
            DIMENSION NXA(4),NYA(4),XPOL(4),YPOL(4)
            CHARACTER*60 CTIT
            DATA NXA/200,1150,200,1150/,
           *     NYA/1600,1600,2700,2700/
            DATA XPOL/0.,0.,0.,0./,YPOL/0.,45.,90.,-45./

            CTIT='Azimuthal Lambert Projections'

            CALL SETPAG('DA4P')
            CALL DISINI
            CALL PAGERA
            CALL COMPLX

            CALL HEIGHT(50)
            NL=NLMESS(CTIT)
            NX=(2250-NL)/2.
            CALL MESSAG(CTIT,NX,300)

            CALL AXSLEN(900,900)
            CALL PROJCT('LAMBERT')

            DO I=1,4
              CALL AXSPOS(NXA(I),NYA(I))
              CALL MAPPOL(XPOL(I),YPOL(I))
              CALL GRAFMP(-180.,180.,-180.,30.,
           *               -90., 90., -90.,30.)

              CALL WORLD
              CALL GRIDMP(1,1)
              CALL ENDGRF
            END DO

            CALL DISFIN
            END

            PROGRAM EX13_4
            PARAMETER (N = 32)
            DIMENSION INRAY(N),IPRAY(N),ICRAY(N)

            DO I=1,N
              INRAY(I)=I
              IPRAY(I)=0
              ICRAY(I)=1
            END DO

            CALL SETPAG('DA4P')
            CALL DISINI
            CALL PAGERA
            CALL COMPLX

            CALL INTAX
            CALL TICKS(1,'XY')
            CALL FRAME(3)
            CALL AXSLEN(1600,2200)
            CALL AXSPOS(400,2700)

            CALL NAME('Longitude','X')
            CALL NAME('Latitude','Y')
            CALL TITLIN('Conformal Conic Projection',3)

            CALL LABELS('MAP','XY')
            CALL PROJCT('CONF')
            CALL GRAFMP(-10.,30.,-10.,5.,35.,70.,35.,5.)

            CALL GRIDMP(1,1)
            CALL SHDEUR(INRAY,IPRAY,ICRAY,N)

            CALL HEIGHT(50)
            CALL TITLE
            CALL DISFIN
            END

