   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
              (Albers).
  = '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.
  = 'MYPR'    defines an user-defined 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 A F R

The routine SHDAFR plots shaded African countries.

The call is:  CALL SHDAFR (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: Algeria           19: Gabon          37: Nigeria
    2: Angola            20: Gambia         38: Rwanda
    3: Benin             21: Ghana          39: Senegal
    4: Botswana          22: Guinea         40: Seychelles
    5: Burkina Faso      23: Guinea Bissau  41: Sierra Leone
    6: Burundi           24: Kenya          42: Somalia
    7: Cameroon          25: Lesotho        43: South Africa
    8: Central Afr. Rep. 26: Liberia        44: Sudan 
    9: Chad              27: Libya          45: Swaziland
   10: Comoros           28: Madagascar     46: Tanzania
   11: Congo, Dem. Rep.  29: Malawi         47: Togo        
   12: Congo, Rep.       30: Mali           48: Tunesia
   13: Cote d'Ivoire     31: Mauritania     49: Uganda
   14: Djibouti          32: Mauritius      50: West Sahara
   15: Egypt             33: Morocco        51: Zambia 
   16: Equatorial Guinea 34: Mozambique     52: Zimbawe
   17: Eritrea           35: Namibia       
   18: Ethiopia          36: Niger           0: Africa

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.

                        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         17: Luxembourg    33: Belarus
         2: Andorra         18: Malta         34: Bosnia 
         3: Belgium         19: Netherlands   35: Croatia
         4: Bulgaria        20: North Ireland 36: Czech 
         5: Germany         21: Norway        37: Estonia
         6: Denmark         22: Austria       38: Latvia
         7: Cyprus          23: Poland        39: Lithuania
         8: United Kingdom  24: Portugal      40: Macedonia
         9: Finland         25: Romania       41: Moldava
        10: France          26: Sweden        42: Russia
        11: Greece          27: Switzerland   43: Serbia         
        12: Ireland         28: Spain         44: Slovakia
        13: Iceland         29: CSFR          45: Slowenia
        14: Italy           30: Turkey        46: Ukraine
        15: Yugoslavia      31: USSR          
        16: Liechtenstein   32: Hungary        0: Europe

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.

Notes:      - The  plotting  of outlines  can be  suppressed
              with CALL NOARLN.
            - To stay  compatible  with older programs,  the 
              number 15 (Yugoslavia) plots Bosnia,  Croatia,
              Macedonia, Serbia and Slowenia,  the number 29
              (CSFR) plots  Czech Republic and Slovakia  and
              the number 31 (USSR)  plots  Belarus, Estonia,
              Latvia,  Lithuania,  Moldava, Russia and Ukra-
              ine.

                        S H D U S A

The routine SHDUSA plots shaded USA states.

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

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

        1: Alabama     19: Maine          37: Oregon 
        2: Alaska      20: Maryland       38: Pennsylvania
        3: Arizona     21: Massachusetts  39: Rhode Island
        4: Arkansas    22: Michigan       40: South Carolina
        5: California  23: Minnesota      41: South Dakota
        6: Colorado    24: Mississippi    42: Tennessee
        7: Connecticut 25: Missouri       43: Texas
        8: Delaware    26: Montana        44: Utah
        9: Floria      27: Nbraska        45: Vermont 
       10: Georgia     28: Nevada         46: Virginia
       11: Hawaii      29: New Hampshire  47: Washington
       12: Idaho       30: New Jersey     48: West Virginia
       13: Ilinois     31: New Mexico     49: Wisconsin
       14: Indiana     32: New York       50: Wyoming
       15: Iowa        33: North Carolina 51: Washington DC 
       16: Kansas      34: North Dakota   
       17: Kentucky    35: Ohio           
       18: Louisiana   36: Oklahoma        0: USA

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

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.  An internal  DISLIN  map file, some external
map files compiled by  Paul Wessel  and map files in  Mapgen
format can be used.  The map files  compiled by  Paul Wessel
can be copied via FTP anonymous from the servers

   ftp://ftp.ngdc.noaa.gov/MGG/shorelines/
   ftp://gmt.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.

Map files in  Mapgen format are available from the Coastline
Extractor:

   http://rimmer.ngdc.noaa.gov/

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.
 = 'MAPFIL'   defines an external map file as base map 
              that is specified with the routine MAPFIL.
                                   Default: CBAS = 'DISLIN'.

                         M A P F I L

The routine  MAPFIL  defines  an external map file.  The map
file can be used as base map if the routine MAPBAS is called
with the parameter 'MAPFIL'.

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

CFIL          is a character string  containing the filename
              of the external map file.
COPT          is a character string describing the format of
              the map coordinates.  COPT can have the values
              'GSHHS' and 'MAPGEN'.

                         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
              'ALL', 'LAND' and 'LAKE'.
                                    Default:  COPT = 'ALL'.

                         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 User-defined Projections

An user-defined  projection can be enabled  with the keyword
'MYPR' in the routine PROJCT. For a user-defined projection,
the user must write a routine that  converts  longitude  and
latitude coordinates  to axis coordinates  (plot coordinates
relative to the origin of the axis system).  The name of the
user written routine  must  then passed  to  DISLIN with the 
routine SETCBK.

                         S E T C B K

The routine SETCBK defines a user written callback routine.

The call is:  CALL SETCBK (ROUTINE, 'MYPR') level 0, 1, 2, 3

ROUTINE       is the name of a routine defined by the user. 
              In Fortran, the routine must be declared as 
              EXTERNAL.

In the following example, a cylindrical projection is imple-
mented as an user-defined projection:

      PROGRAM MYPR
      EXTERNAL MYFUNC
      COMMON /MYCOMM/ XA,XE,YA,YE,NXL,NYL

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

      NXL = 2400
      NYL = 1200

      CALL METAFL ('cons')
      CALL DISINI
      CALL PAGERA 
      CALL COMPLX 
      CALL AXSLEN (NXL, NYL)

      CALL PROJCT ('MYPR') 
      CALL SETCBK (MYFUNC, 'MYPR')

      CALL GRAFMP (XA, XE, XA, 90., YA, YE, YA, 30.)
      CALL GRIDMP (1,1)
      CALL WORLD 
      CALL DISFIN 
      END

      SUBROUTINE MYFUNC (XP, YP)  
      COMMON /MYCOMM/ XA,XE,YA,YE,NXL,NYL
      XP = (XP - XA)/(XE - XA) * (NXL - 1)
      YP = (YP - YA)/(YE - YA) * (NYL - 1)
      END

13.8 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
            DIMENSION INRAY(1),IPRAY(1),ICRAY(1)

            INRAY(1)=0
            IPRAY(1)=0
            ICRAY(1)=1

            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,1)

            CALL HEIGHT(50)
            CALL TITLE
            CALL DISFIN
            END

