How to generate a background map.

[Steffen Weber]

Here is a small programm written by Dr. Rupa Kumar Kolli.
It was posted on the list on Thu, 17 Feb 2005 17:48:46

Hope this helps:


       PROGRAM GRADS_MAP
c
c     Program to prepare customized background map for GrADS from a text
c     file containing line segments in (longitude,latitude) coordinates.
c
c     By Rupa Kumar Kolli, Indian Institute of Tropical Meteorology
c        kolli at tropmet.res.in
c        created:  December 10, 1997
c        modified: March 12, 1998
c        PC Fortran support added: January 26, 1999
c        minor adjustments: March 20, 2002
c     Please do contact me for comments/feedback etc. on this program.
c
c     This program converts map files in ASCII into GrADS readable format.
c     The resolution of the map background, of course, depends on the
c     basic data that you can get.  The format of the GrADS map data set
c     prepared by this program conforms to that used in 'lowres'.  While it
c     should work for any kind of resolution, it may be a bit slower than
c     the standard 'hires' for huge files.
c
c     Input File Format:
c     The input file can contain any number of line segments.
c     Each segment starts with a header record containing the
c         the number of pairs of coordinates and associated map type
c         (see format statement 1)
c         (the map type is used to selectively display map lines using
c          the 'set mpt ...' command)
c     The coordinates should be listed from the subsequent record,
c         as longitude,latitude pair, one pair per line
c         (see format statement 2).
c
c     Example:
c
c           5       2
c     209.830 -76.580
c     211.500 -76.750
c     212.830 -76.420
c     211.830 -76.170
c     209.830 -76.580
c
c     VERY IMPORTANT: NOTE THE FOLLOWING !!!
c
c     1. Longitudes should be in the range 0 to 360 (No negative values !).
c
c     2. Latitudes should be in the range -90 (South Pole) to 90 (North Pole).
c
c     3. DON'T CROSS THE GREENWICH MERIDIAN (0.0 or 360.0) IN ANY SEGMENT.
c        For lines across this meridian, split them into two parts, with
c        a common point between them.  Use the value 360.0 for Greenwich
c        meridian if you are moving towards it from west to east, and the
c        value 0.0 if you are moving from east to west.  That is, you
c        cannot go from 359.0 to 0.0; this will wrap the line across the
c        whole globe !
c
c     4. The maximum number of points (pairs of coordinates) in a single
c        segment is 255.
c
c     5. Don't use GrADS standard file names (lowres, mres, hires, nam) for
c        the output file (e.g., rupres). Preserve them for standard use.
c
c     6. After generating the GrADS map background file, copy it to a path
c        where GrADS looks for its support files (generally the path defined
c        in $GADDIR)
c
c     7. Finally, start GrADS and give the command 'set mpdset filename':
c        e.g., set mpdset rupres.  Note that you can keep multiple map
c        data sets open in GrADS 1.7; thus, you can create separate map
c        data files and selectively use subsets of the files.  For this
c        purpose, use the command like 'set mpdset map1 map2', whereupon
c        GrADS uses both the files for drawing the map.
c
c     8. This program was initially developed on a HP 9000/720 machine with
c        HP-UX 9.05 operating system.  I have successfully tested it in dos/
c        Windows as well as in RedHat Linux (g77). It should work on other
c        platforms as well.  For using this program on a PC (dos/windows/
c        linux), a slight change is necessary, as indicated in the comments
c        preceding the output statements (marked by --->) in the program.
c
c     9. The GrADS map file generated by this program is platform-independent,
c        and it can be used across platforms (I didn't check Mac, though).
c
c    10. If you have successfully generated a specialized GrADS map that may
c        be of interest to other GrADS fans, please consider offering it for
c        sharing, possibly by making a suitable announcement in GrADS mailing
c        list.
c
c     ENJOY HAVING YOUR OWN CUSTOMIZED GrADS MAP !!!



       byte ihead(3),ilon(4),ilat(4)
       integer lon,lat
       equivalence (lon,ilon), (lat,ilat)
c
c     Open the input and output files
c
       open(10,file='map.txt',
      >     status='old')
       open(20,file='map',
      >     form='unformatted',access='direct',recl=3)
c
c     Initialize the first header-record byte, ihead(1): fixed to be 1
c     A value of 2 can be used to skip unwanted sections of the file,
c     (mres and hires were created like this), but this facility
c     is not implemented in this program.  So don't change this value
c     unless you want to implement skipping (could be fairly complicated !)
c
       ihead(1)=1
       irec=0
       nsegs=0
c
c     Read the number of pairs and map type for this segment.
c     Leave the loop if end of file is reached.
c
    99 read(10,1,end=100)npairs,ihead(2)
     1 format(2i6)
       print *,npairs,ihead(2)
       nsegs=nsegs+1
       if (npairs.gt.255) then
         print *,'Number of pairs in Segment No.',nsegs,' exceeds 255 !'
         stop 'Job Aborted'
       endif
c
c     Write the header record (3 bytes)
c     The third byte indicates the number of pairs
c     of values following (maximum=255)
c
       ihead(3)=npairs
       irec=irec+1
       write(20,rec=irec)ihead
c
c     Loop over the current segment
c
       do 10 i=1,npairs
c
c     Read one pair of coordinates (Longitude,Latitude)
c
         read(10,2)alon,alat
     2   format(2f6.2)
c
c     Convert the values to integers, in GrADS readable format
c
         lon=int(alon*10000.)+0.5
         lat=int((alat+90.)*10000.)+0.5
c
c     For use with PCs running DOS/Windows/Linux and Alpha workstations:
c     Write the first three-bytes of lon and lat on the output file,
c     in reverse order.
c     Comment out the statements meant for Other Workstations.
c
c     For use with Other Workstations (e.g., HP-UX, IRIX, SOLARIS):
c     Write the last three-bytes of lon and lat on the output file.
c     Comment out the statements meant for PCs and Alphas.
c
c     Note that ilon and ilat are byte equivalents of lon and lat.
c
         irec=irec+1
c --->  For PC (dos/windows/linux) and Alpha workstations:
         write(20,rec=irec)(ilon(k),k=3,1,-1)
c --->  For Other Workstations:
c        write(20,rec=irec)(ilon(k),k=2,4)
         irec=irec+1
c --->  For PC (dos/windows/linux) and Alpha workstations:
         write(20,rec=irec)(ilat(k),k=3,1,-1)
c --->  For Other Workstations:
c        write(20,rec=irec)(ilat(k),k=2,4)
    10 continue
c
c     Go back to read next line segment
c
       goto 99
   100 print *,'Total Number of Segments Processed=',nsegs
       close (10)
       close (20)
       stop 'Normal Termination'
       end