C C $Id: mpex13.f,v 1.1 2007/09/24 23:25:48 haley Exp $ C PROGRAM TESTIT C C This program shows how to use EZMAP routines providing access to C the RANGS database ("Regionally Accessible Nested Global Shorelines") C created by Reiner Feistel, which make use of the GSHHS database C ("Global Self-consistent Heirarchical High-resolution Shorelines") C created by Paul Wessel and Walter H. F. Smith. The original data were C from the World Vector Shoreline and World Databank II databases. C C It also demonstrates the use of routines for labelling latitudes or C longitudes along a specified line on a map. C C Note that, as it stands, this program will fail to produce any plots C using the RANGS/GSHHS data, because that data will not be available C to it. The user must add a routine called MDRGDI, which looks like C this: C C SUBROUTINE MDRGDI (DINM) C C This is a user-replaceable routine that returns the name of the C directory in which the RANGS/GSHHS data files have been placed. C C CHARACTER*(*) DINM C C Return the name of the directory where the RANGS/GSHHS data reside. C C DINM='/fs/scd/home1/kennison/Mapping/RANGSData' C C Done. C C RETURN C C END C C The specified directory must contain the following files: C C gshhs(0).rim 89496428 bytes C gshhs(1).rim 20775224 bytes C gshhs(2).rim 5099196 bytes C gshhs(3).rim 1109192 bytes C gshhs(4).rim 170948 bytes C rangs(0).cat 259200 bytes C rangs(0).cel 6518664 bytes C rangs(1).cat 259200 bytes C rangs(1).cel 5945492 bytes C rangs(2).cat 259200 bytes C rangs(2).cel 4062215 bytes C rangs(3).cat 259200 bytes C rangs(3).cel 3398642 bytes C rangs(4).cat 259200 bytes C rangs(4).cel 3045848 bytes C C These files may be obtained from Rainer Feistel's Web site: C C http://www.io-warnemuende.de/homepages/rfeistel/index.html C C A description of the files is available here: C C http://www.agu.org/eos_elec/99063e.html C C Hopefully, these files will also be available from the NCAR Graphics C Web site. C C Define error file, Fortran unit number, and workstation type, C and workstation ID. C PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1, IWKID=1) C C Declare arrays to hold center-point positions, angular half-widths, C grid spacings, resolutions, and titles for four different maps to C be drawn. For each of the four, the map is drawn three different C times, once using the old 'CO' database, once using the database C "Earth..2", and again using the RANGS/GSHHS database. C DOUBLE PRECISION PLAT(4),PLON(4),ANGL(4),GRSP(4) INTEGER IRGL(4) CHARACTER*16 TITL(4) C C Declare arrays to hold color indices to be passed to MDRGSC, defining C the colors to be used on maps drawn from the RANGS/GSHHS database C DIMENSION ICOL(5),ICSF(5) C C Declare a character-size variable. C REAL CHSZ C C Declare stuff required to do color fill. C PARAMETER (LAMA=200000,NCRA=10000,NGPS=10,LRWK=2*NCRA) C DIMENSION IAMA(LAMA),XCRA(NCRA),YCRA(NCRA),IAAI(NGPS),IAGI(NGPS) C DIMENSION RWRK(LRWK) EQUIVALENCE (RWRK(1),XCRA(1)),(RWRK(NCRA+1),YCRA(1)) C EXTERNAL COLRCO,COLRE2 C C Define the center-point positions and angular half-width for each of C the four plots. C DATA PLAT(1),PLON(1),ANGL(1) / 23.0D0 , -76.0D0 , 5.0D0 / DATA PLAT(2),PLON(2),ANGL(2) / -53.5D0 , -70.0D0 , 3.0D0 / DATA PLAT(3),PLON(3),ANGL(3) / 54.5D0 , -2.0D0 , 4.5D0 / DATA PLAT(4),PLON(4),ANGL(4) / 37.0D0 , 126.5D0 , 4.0D0 / C C Define the resolution of the grid and the level of the RANGS data C to be used on each of the four plots. C DATA GRSP(1),IRGL(1) / 1.D0 , 0 / DATA GRSP(2),IRGL(2) / 1.D0 , 0 / DATA GRSP(3),IRGL(3) / 1.D0 , 0 / DATA GRSP(4),IRGL(4) / 1.D0 , 0 / C C Define titles for the tops of the plots, describing the place mapped. C DATA TITL(1) / 'Caribbean' / DATA TITL(2) / 'Tierra del Fuego' / DATA TITL(3) / 'Great Britain' / DATA TITL(4) / 'Korea' / C C Define the character size to be used for lat/lon labels. C DATA CHSZ / .014 / C C Open GKS. C CALL GOPKS (IERRF, ISZDM) CALL GOPWK (IWKID, LUNIT, IWTYPE) CALL GACWK (IWKID) C C Turn off clipping by GKS (so labels don't get clipped). C CALL GSCLIP (0) C C Define some colors to use. C CALL GSCR (1,0,0.,0.,0.) ! Default background CALL GSCR (1,1,1.,1.,1.) ! Default foreground C CALL GSCR (1,2,0.,0.,1.) ! RANGS/GSHHS - Ocean CALL GSCR (1,3,0.,1.,0.) ! RANGS/GSHHS - Land CALL GSCR (1,4,0.,1.,1.) ! RANGS/GSHHS - Lake CALL GSCR (1,5,0.,1.,0.) ! RANGS/GSHHS - Island in lake CALL GSCR (1,6,0.,1.,1.) ! RANGS/GSHHS - Pond on island C CALL GSCR (1,7,1.,0.,0.) ! Debug C C Tell EZMAP to draw the map in a bit smaller square, leaving room for C lat/lon labels around the edge. C CALL MAPPOS (.08,.92,.08,.92) C C Tell EZMAP what colors to use for the outlines and filled areas C produced from the RANGS/GSHHS data. C ICOL(1)=7 ICOL(2)=7 ICOL(3)=7 ICOL(4)=7 ICOL(5)=7 C ICSF(1)=2 ICSF(2)=3 ICSF(3)=4 ICSF(4)=5 ICSF(5)=6 C CALL MDRGSC (ICOL,ICSF) C C Increase the number of points to be interpolated along the edges of C the 1-degree cells (to get rid of some edge effects that otherwise C occur). C CALL MPSETI ('II',64064) C CALL MPSETD ('GD',.1D0) C C Set the PLOTCHAR font, character color, and outlining. C CALL PCSETI ('FN',26) CALL PCSETI ('CC', 1) CALL PCSETI ('OF', 1) CALL PCSETI ('OC', 1) C C Loop to do views of four different spots on the globe. C DO 101 ISPT=1,4 C C Print the loop index, so the user will know something is happening. C (This program runs pretty slowly.) C PRINT * , ' ' PRINT * , 'POSITION NUMBER ',ISPT,' ',TITL(ISPT) C C Define a satellite altitude (for use if 'SV' is selected below). C CALL MAPSTD ('SA',2.5D0) C C Define the projection to be used. (If only a small portion of the C projection plane is being shown, it doesn't make too much difference C which projection is chosen.) C C CALL MDPROJ ('ST',PLAT(ISPT),PLON(ISPT),0.D0) CALL MDPROJ ('OR',PLAT(ISPT),PLON(ISPT),0.D0) C CALL MDPROJ ('CE',PLAT(ISPT),PLON(ISPT),0.D0) C CALL MDPROJ ('MO',PLAT(ISPT),PLON(ISPT),0.D0) C CALL MDPROJ ('ME',PLAT(ISPT),PLON(ISPT),0.D0) C CALL MDPROJ ('RO',PLAT(ISPT),PLON(ISPT),0.D0) C CALL MDPROJ ('AE',PLAT(ISPT),PLON(ISPT),0.D0) C CALL MDPROJ ('LE',PLAT(ISPT),PLON(ISPT),0.D0) C CALL MDPROJ ('GN',PLAT(ISPT),PLON(ISPT),0.D0) C CALL MDPROJ ('SV',PLAT(ISPT),PLON(ISPT),0.D0) C C Define the limits of the map. C CALL MDPSET ('AN',ANGL(ISPT),ANGL(ISPT),ANGL(ISPT),ANGL(ISPT)) C C Here's an interesting combination: C C CALL MAPSTD ('SA',1.03125D0) C CALL MDPROJ ('SV',PLAT(ISPT),PLON(ISPT),0.D0) C CALL MDPSET ('AN', 65.D0, 65.D0, 65.D0, 65.D0) C C Tell EZMAP what spacing to use on the grid. C CALL MAPSTD ('GR',GRSP(ISPT)) C C Draw a default EZMAP map, using the outline database 'CO', and then C add a grid, lat/lon labels, and informational labels. C PRINT * , ' USING DATABASE "CO"' CALL ARINAM (IAMA,LAMA) CALL MDPBLA (IAMA) CALL ARSCAM (IAMA,XCRA,YCRA,NCRA,IAAI,IAGI,NGPS,COLRCO) CALL GSPLCI (7) CALL MDPLOT CALL GSPLCI (1) CALL MDPGRD CALL MDLBLT (.08,.08,.08,.92,-CHSZ, 0.,CHSZ,0.,+1.) CALL MDLBLT (.92,.08,.92,.92, CHSZ, 0.,CHSZ,0.,-1.) CALL MDLBLN (.08,.08,.92,.08,0.,-1.5*CHSZ,CHSZ,0., 0.) CALL MDLBLN (.08,.92,.92,.92,0.,+1.5*CHSZ,CHSZ,0., 0.) CALL PLCHHQ (CFUX(.5),CFUY(.975), + TITL(ISPT)(1:MDILNB(TITL(ISPT))),.015,0.,0.) CALL PLCHHQ (CFUX(.5),CFUY(.025),'Old Database "CO"', + .015,0.,0.) CALL NGSETR('OY',0.03) CALL NGEZLOGO CALL FRAME C C Draw a map using outlines from the database "Earth..2", and then C add a grid, lat/lon labels, and informational labels. C PRINT * , ' USING DATABASE "Earth..2"' CALL ARINAM (IAMA,LAMA) CALL MDLNAM ('Earth..2',1,IAMA) CALL ARSCAM (IAMA,XCRA,YCRA,NCRA,IAAI,IAGI,NGPS,COLRE2) CALL GSPLCI (7) CALL MDLNDR ('Earth..2',1) CALL GSPLCI (1) CALL MDPGRD CALL MDLBLT (.08,.08,.08,.92,-CHSZ, 0.,CHSZ,0.,+1.) CALL MDLBLT (.92,.08,.92,.92, CHSZ, 0.,CHSZ,0.,-1.) CALL MDLBLN (.08,.08,.92,.08,0.,-1.5*CHSZ,CHSZ,0., 0.) CALL MDLBLN (.08,.92,.92,.92,0.,+1.5*CHSZ,CHSZ,0., 0.) CALL PLCHHQ (CFUX(.5),CFUY(.975), + TITL(ISPT)(1:MDILNB(TITL(ISPT))),.015,0.,0.) CALL PLCHHQ (CFUX(.5),CFUY(.025),'Earth..2',.015,0.,0.) CALL NGEZLOGO CALL FRAME C C Now, draw a solid-filled map using the RANGS/GSHHS data. C PRINT * , ' USING RANGS/GSHHS DATABASE' C CALL MDRGSF (IRGL(ISPT),RWRK,LRWK,IAMA,LAMA) C C Add outlines using the RANGS/GSHHS data. C CALL MDRGOL (IRGL(ISPT),RWRK,LRWK) C C Overlay a lat/lon grid. C CALL MDPGRD C C Add latitude and longitude labels. C CALL MDLBLT (.08,.08,.08,.92,-CHSZ, 0.,CHSZ,0.,+1.) CALL MDLBLT (.92,.08,.92,.92, CHSZ, 0.,CHSZ,0.,-1.) CALL MDLBLN (.08,.08,.92,.08,0.,-1.5*CHSZ,CHSZ,0., 0.) CALL MDLBLN (.08,.92,.92,.92,0.,+1.5*CHSZ,CHSZ,0., 0.) C C Put a label at the top of the plot. C CALL PLCHHQ (CFUX(.5),CFUY(.975), + TITL(ISPT)(1:MDILNB(TITL(ISPT))),.015,0.,0.) C C Put a label at the bottom of the plot. C CALL PLCHHQ (CFUX(.5),CFUY(.025),'RANGS/GSHHS',.015,0.,0.) C C Advance the plotter frame. C CALL NGEZLOGO CALL FRAME C C End of loop. C 101 CONTINUE C C Reset the PLOTCHAR font, character color, and outlining. C CALL PCSETI ('FN', 0) CALL PCSETI ('CC',-1) CALL PCSETI ('OF', 0) CALL PCSETI ('OC', 1) C C Close GKS. C CALL GDAWK (IWKID) CALL GCLWK (IWKID) CALL GCLKS C C Done. C STOP C END SUBROUTINE COLRCO (XCRA,YCRA,NCRA,IAAI,IAGI,NGPS) DIMENSION XCRA(NCRA),YCRA(NCRA),IAAI(NGPS),IAGI(NGPS) IAI1=-1 DO 101 IGRP=1,NGPS IF (IAGI(IGRP).EQ.1) IAI1=IAAI(IGRP) 101 CONTINUE IF (IAI1.GT.0) THEN IF (MDPACI(IAI1).EQ.1) THEN CALL GSFACI (2) ! (OCEAN) ELSE CALL GSFACI (3) ! (LAND) END IF CALL GFA (NCRA-1,XCRA,YCRA) END IF RETURN END SUBROUTINE COLRE2 (XCRA,YCRA,NCRA,IAAI,IAGI,NGPS) DIMENSION XCRA(NCRA),YCRA(NCRA),IAAI(NGPS),IAGI(NGPS) IAI1=-1 DO 101 IGRP=1,NGPS IF (IAGI(IGRP).EQ.1) IAI1=IAAI(IGRP) 101 CONTINUE IF (IAI1.GT.0) THEN IF (MDIPAN(IAI1,'Water').NE.0) THEN CALL GSFACI (2) ! (OCEAN) ELSE CALL GSFACI (3) ! (LAND) END IF CALL GFA (NCRA-1,XCRA,YCRA) END IF RETURN END