home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume4 / wdb-ii < prev    next >
Encoding:
Text File  |  1989-02-03  |  17.7 KB  |  706 lines
  1. Path: xanth!mcnc!rutgers!apple!bionet!ig!agate!pasteur!ames!necntc!ncoast!allbery
  2. From: awpaeth@watcgl.waterloo.edu (Alan Wm Paeth)
  3. Newsgroups: comp.sources.misc
  4. Subject: v04i087: WDBII mapping software (10:1 compress, Unix PLOT)
  5. Keywords: digital cartography world databank database mapping
  6. Message-ID: <5886@watcgl.waterloo.edu>
  7. Date: 20 Sep 88 01:23:24 GMT
  8. Sender: allbery@ncoast.UUCP
  9. Reply-To: awpaeth@watcgl.waterloo.edu (Alan Wm Paeth)
  10. Organization: U of Waterloo, Ontario
  11. Lines: 692
  12. Approved: allbery@ncoast.UUCP
  13.  
  14. Posting-number: Volume 4, Issue 87
  15. Submitted-by: "Alan Wm Paeth" <awpaeth@watcgl.waterloo.edu>
  16. Archive-name: wdb-ii
  17.  
  18. I have had a number of requests for my World Data Bank II mapping
  19. software (C Unix). For sites that have all five mag tape volumes of
  20. this dataset (don't come crying to me!), this software package provides
  21. two C programs. The first does a 10:1 compression on the datasets and
  22. additionally maps the straight vector detail into an index file of offsets,
  23. plus a companion vector file (two bytes per vector). The second program
  24. draws a rectangular projection of the dataset for specified lat,long,scale
  25. and produces Unix PLOT output. It can be very easily upgraded. Enjoy.
  26.  
  27.    /Alan Paeth
  28.    Computer Graphics Laboratory
  29.    University of Waterloo
  30.  
  31. -----------------------------------------------------------------------------
  32.  
  33. # This is a shell archive.  Remove anything before this line,
  34. # then unpack it by saving it in a file and typing "sh file".
  35. #
  36. # Wrapped by watcgl!awpaeth on Sat Sep 17 18:28:33 EDT 1988
  37. # Contents:  README Makefile wdbread.c wdbplot.c
  38.  
  39. echo x - README
  40. sed 's/^@//' > "README" <<'@//E*O*F README//'
  41. README    -- this file
  42.  
  43. MAKEFILE  -- compile both sources
  44.  
  45. wdbread.c -- (see source code) converts an EBCDIC WDBII dataset (20 char recs)
  46.          into two output files: map.index and map.data, while producing
  47.          a large amount of line and summary detail. The first file provides
  48.          bounding box and # of vertex information into the second file
  49.          (map.data), this latter file contains (dx,dy) byte pairs. wdbread
  50.          will fragment long vectors as needed to provide this format, and
  51.          will report this on the log it produces on stdout.
  52.  
  53.          Multiple runs of this program will place (overwrite, if you're
  54.          not careful) map/index pairs for the 16 or so volumes which
  55.          comprise all of WDBII (arranged by continents and detail).
  56.          All .index and .data files may then be "cat"ed into one giant
  57.          file which represents the entire world. This is possible by
  58.          design: the skip offsets in the .index file are relative byte
  59.          counts in each companion .data file, so each file can be
  60.          appended to in tandem.
  61.          
  62. wdbplot.c -- given location and scale produces UNIX PLOT output. The last
  63.          section of source code provides a simple means to customize
  64.          the software to direct drive other devices, although we most
  65.          often use Tektronix output (via plot). Other suggested fixes
  66.          would be to "trigger" based on vector rank, which changes as
  67.          a function of political boundary, river, coastline, etc. More
  68.          ambitious users might wish to encode more general map projections:
  69.          Mercator and Rectangular equal area would be the most simple
  70.          to employ because the clipping borders would remain roughly
  71.          unchanged.
  72.  
  73. Author's disclamer: I have no intent in modifying or upgrading this ancient
  74. code: if you can make big bucks on it, go for it. If you want to discuss
  75. digital cartography from a less nuts and bolts level, drop me a line sometime!
  76.  
  77. -- Alan Wm Paeth (awpaeth@watcgl.waterloo.edu, awpaeth@watcgl.waterloo.cdn).
  78. @//E*O*F README//
  79. chmod u=rw,g=r,o=r README
  80.  
  81. echo x - Makefile
  82. sed 's/^@//' > "Makefile" <<'@//E*O*F Makefile//'
  83. CFLAGS = -g
  84.  
  85.  
  86. wdbplot: wdbplot.c
  87.     cc $(CFLAGS) wdbplot.c -lplot -lm -o $@
  88.  
  89. wdbread: wdbread.c
  90.     cc $(CFLAGS) wdbread.c -o $@
  91. @//E*O*F Makefile//
  92. chmod u=rw,g=r,o=r Makefile
  93.  
  94. echo x - wdbread.c
  95. sed 's/^@//' > "wdbread.c" <<'@//E*O*F wdbread.c//'
  96. /* wdbread - converts raw world data bank II files into vax words
  97.  
  98. /* NOTE: like old wdbiiread, but limits consec. points to +-127 units */
  99.  
  100. /* written by Alan Paeth, University of Waterloo, January, 1984 */
  101.  
  102. /* wdbread <in >out 
  103.    reads in hunks of twenty byte EBCDIC records, as per the FORTRAN 
  104.    spec of WDBII, and writes them out as signed 8-bit integers (signed
  105.    bytes), as delta values. No conversion of characters is done for other
  106.    than the digits and few other expected sparse characters.
  107.  
  108. NOTE: the output is named map.data and map.index, both placed on the
  109. current directory. The file >out contains a summary of processed records.
  110. */
  111.  
  112. /* DEFINITIONS */
  113.  
  114. #include <stdio.h>
  115.  
  116. #define PROGNAME "wdbread"
  117. #define EBCDICFLAG 0        /* 0=ascii input, 1=ebcdic */
  118. #define INBUFSIZE 20
  119. #define PIPEOUT 1
  120. #define DATANAME "map.data"
  121. #define INDEXNAME "map.index"
  122. #define WRITEMODE "w"
  123. #define OPENFAIL 0 
  124. #define LON360 1296000
  125. #define FRAGLIMIT 20
  126. #define maxint  2147483647
  127. #define minint -2147483648 
  128. #define v(arg) (inbuf[arg]-'0')
  129.  
  130. #define ABS(a) ((a)<0?(-a):(a))
  131.  
  132. FILE *mapdata, *mapindex;
  133. char inbuf[INBUFSIZE+1], remap[256];
  134. int rank, count;
  135. int fragcount, totfrag, readrec, writerec, readpoint, writepoint;
  136. int i, j, k, n, s, e, w;
  137. short shortn, shorts, shorte, shortw;
  138. int olat, olon, privlat, privlon, dlat, dlon, oprivlat, oprivlon;
  139.  
  140. /* THE MAIN */
  141.  
  142. main(argc,argv)
  143.     int  argc;
  144.     char **argv;
  145.     {
  146.  
  147. /* open the data and index output */
  148.  
  149.     if ((mapdata = fopen(DATANAME, WRITEMODE)) == OPENFAIL)
  150.         abort("open fail on %s", DATANAME);
  151.  
  152.     if ((mapindex = fopen(INDEXNAME, WRITEMODE)) == OPENFAIL)
  153.         abort("open fail on %s", INDEXNAME);
  154.     
  155. /* set remap to map into ASCII, digit 0 for all other characters */
  156.  
  157.     if (EBCDICFLAG) 
  158. /* EBCDIC REMAPPING STUFF */
  159.         {
  160.     for (i=0; i<256; i++) remap[i] = '0';        /* def. 0 */
  161.     for (i=0; i<10; i++) remap[240+i] = 48+i;    /* digits */
  162.     remap[103] = 'W';                /* w,W,s,S */
  163.     remap[231] = 'W';
  164.     remap[99] = 'S';
  165.     remap[227] = 'S';
  166.     }
  167.     else
  168. /* ASCII REMAPPING (identity xform, space -> digit 0 */
  169.     {
  170.     for (i=0; i<256; i++) remap[i] = i;        /* identity */
  171.     remap[' '] = '0';                /* space -> 0 */
  172.     }
  173.  
  174.     inbuf[INBUFSIZE] = '\0';
  175.     readrec = 0;
  176.     writerec = 0;
  177.     totfrag = 0;
  178.  
  179. /* now the top level count loop */
  180.  
  181.     while( fread(&inbuf[0], 1, INBUFSIZE, stdin) == INBUFSIZE)
  182.     {
  183.     readrec++;
  184.     readpoint = 0;
  185.     writepoint = 0;
  186.     
  187.     for (i=0; i<INBUFSIZE; i++) inbuf[i] = remap[inbuf[i]];
  188.     printf("/%s/\n", inbuf);
  189.     sscanf( &inbuf[9], "%6ld", &count);
  190.     sscanf( &inbuf[7], "%2ld", &rank);
  191.     printf("count %d rank %d\n", count, rank);
  192.     fflush(stdout);
  193.     
  194.     n = minint;
  195.     s = maxint;
  196.     e = minint;
  197.     w = maxint;
  198.  
  199. /* NEW STUFF */
  200.  
  201.     fragcount = 0;
  202.     
  203. /* the high-speed data loop */
  204.     for (j=0; j<count; j++)
  205.         {
  206. /* optimizations */
  207.         register int lat, lon, rSIX, rTEN;
  208.         register char *base;
  209.         rTEN = 10;
  210.         rSIX = 6;
  211.  
  212.         olat = lat;
  213.         olon = lon;
  214.         
  215.             fread(&inbuf[0], 1, INBUFSIZE, stdin);
  216.         for (i=0; i<INBUFSIZE; i++) inbuf[i] = remap[inbuf[i]];
  217.         
  218.         readpoint++;
  219. /*
  220.  * optimization of the following:
  221.  *
  222.  *        lat = v(0)*36000 + v(1)*3600 + v(2)*600 +
  223.  *          v(3)*60 + v(4)*10 + v(5);
  224.  *        if (inbuf[6] == 'S') lat = -lat;
  225.  */
  226.         base = &inbuf[0];
  227.         lat = *base++;
  228.         lat *= rTEN;
  229.         lat += *base++;
  230.         lat *= rSIX;
  231.         lat += *base++;
  232.         lat *= rTEN;
  233.         lat += *base++;
  234.         lat *= rSIX;
  235.         lat += *base++;
  236.         lat *= rTEN;
  237.         lat += *base++;
  238.         lat -= '0'*(36000+3600+600+60+10+1); /* correct for ascii 000000 */
  239.         if (*base == 'S') lat = -lat;
  240. /*
  241.  * optimization of the following:
  242.  *
  243.  *        lon = v(7)*360000 + v(8)*36000 + v(9)*3600 +
  244.  *          v(10)*600 + v(11)*60 + v(12)*10 + v(13);
  245.  *          if (inbuf[14] == 'W') lon = -lon;
  246.  */
  247.         base = &inbuf[7];
  248.         lon = *base++;
  249.         lon *= rTEN;
  250.         lon += *base++;
  251.         lon *= rTEN;
  252.         lon += *base++;
  253.         lon *= rSIX;
  254.         lon += *base++;
  255.         lon *= rTEN;
  256.         lon += *base++;
  257.         lon *= rSIX;
  258.         lon += *base++;
  259.         lon *= rTEN;
  260.         lon += *base++;
  261.         lon -= '0'*(360000+36000+3600+600+60+10+1); /* ascii 0000000 */
  262.         if (*base == 'W') lon = -lon;
  263. /*
  264.  * hemisphere check
  265.  */
  266.         if (lon >= LON360) lon -= LON360;
  267.         if ((j != 0) && (ABS(olon-lon)>(LON360/2)))
  268.         {
  269.         warning("date-line crossing: %d\n", lon);
  270.         }
  271. /*
  272.  * adjust bounding box
  273.  */
  274.         if (lat > n) n = lat;
  275.         if (lat < s) s = lat;
  276.         if (lon > e) e = lon;
  277.         if (lon < w) w = lon;
  278.         
  279.             if (j == 0)
  280.         {
  281.         dlat = lat;
  282.         dlon = lon;
  283.         fwrite(&dlat, 1, 4, mapdata);
  284.         fwrite(&dlon, 1, 4, mapdata);
  285.         writepoint++;
  286.         }
  287.         else
  288.         {
  289.         dlat = lat - olat;
  290.             dlon = lon - olon;
  291.             fragcount = roundupdiv( lmax( labs(dlat), labs(dlon) ) , 127l);
  292.         if (fragcount ==1)
  293.             {
  294.             fwrite(&dlat, 1, 1, mapdata);
  295.             fwrite(&dlon, 1, 1, mapdata);
  296.             writepoint++;
  297.             }
  298.         else
  299.             {
  300.             if (fragcount == 0)
  301.                 warning("fragcount=0 => consec. duplicate points",NULL);
  302.             if (fragcount > FRAGLIMIT)
  303.                 {
  304.             warning("excessive fragments %ld",fragcount);
  305.             warning("occurring at segment: %d\n", j);
  306.             if (fragcount > 2000)
  307.                 {
  308.                 warning("trace vars: \n");
  309.                 warning(" lat = %d\n",  lat);
  310.                 warning("olat = %d\n", olat);
  311.                 warning("dlat = %d\n", dlat);
  312.                 warning(" lon = %d\n",  lon);
  313.                 warning("olon = %d\n", olon);
  314.                 warning("dlon = %d\n", dlon);
  315.                 }
  316.             }
  317.             oprivlat = 0;
  318.             oprivlon = 0;
  319.             for (k = 0; k < fragcount; k++)
  320.                 {
  321.                 privlat = rounddiv( dlat*(k+1), fragcount) - oprivlat;
  322.                 privlon = rounddiv( dlon*(k+1), fragcount) - oprivlon;
  323.                 oprivlat = privlat;
  324.                 oprivlon = privlon;
  325.  
  326.                 fwrite( &privlat, 1, 1, mapdata);
  327.                 fwrite( &privlon, 1, 1, mapdata);
  328.  
  329.                 writepoint++;
  330.                 }
  331.             printf(" frg=%ld dx=%ld dy=%ld\n", fragcount,dlat,dlon);
  332.             totfrag += (fragcount - 1);
  333.             }
  334.             } 
  335.         }
  336.         fwrite( &writepoint, 1, 2, mapindex);
  337.         fwrite( &rank, 1, 2, mapindex);
  338.     writerec++;
  339.  
  340. /* round n,e +, round s,w - */
  341.  
  342.     shortn = round20up(n);
  343.     shorts = round20down(s);
  344.     shorte = round20up(e);
  345.     shortw = round20down(w);
  346.  
  347.     fwrite( &shortn, 1, 2, mapindex);
  348.         fwrite( &shorts, 1, 2, mapindex);
  349.         fwrite( &shorte, 1, 2, mapindex);
  350.         fwrite( &shortw, 1, 2, mapindex);
  351.  
  352.     printf("read=%ld write=%ld n=%ld s=%ld e=%ld w=%ld\n\n", readpoint,
  353.        writepoint, n, s, e, w);
  354.     fflush(stdout);
  355.     }
  356.     printf("\n\nTOTALS: records in=%ld records out=%ld new fragment=%ld\n",
  357.         readrec,writerec,totfrag);
  358.     fclose(mapdata);
  359.     fclose(mapindex);
  360.     exit(0);
  361.     }
  362.  
  363. /* ROUTINES */
  364.  
  365. abort(a,b)
  366.     char *a, *b;
  367.     {
  368.     fprintf(stderr,"%s error: ",PROGNAME);
  369.     fprintf(stderr,a,b);
  370.     fprintf(stderr,"\n");
  371.     exit(1);
  372.     }
  373.  
  374. warning(a,b)
  375.     char *a, *b;
  376.     {
  377.     printf("\nWARNING: ");
  378.     printf(a,b);
  379.     printf("\n");
  380.     fflush(stdout);
  381.     }
  382.  
  383. roundupdiv(a,b)
  384.     int a,b;
  385.     {
  386.     return (a >= 0) ? ((a - 1) / b + 1) : -roundupdiv(-a,b);
  387.     }
  388.  
  389. rounddiv(a,b)
  390.     int a,b;
  391.     {
  392.     return (a >= 0) ? (2*a + b) / (2*b) : -rounddiv(-a,b);
  393.     }
  394.  
  395. lmax(a,b)
  396.     int a,b;
  397.     {
  398.     return ((a >= b) ? a : b);
  399.     }
  400.  
  401. labs(a)
  402.     int a;
  403.     {
  404.     return ((a >= 0) ? a : -a);
  405.     }
  406.  
  407. round20up(a)
  408.     int a;
  409.     {
  410.     return ( (a >= 0) ? (a+19)/20 : (-round20down(-a)) );
  411.     }
  412.  
  413. round20down(a)
  414.     int a;
  415.     {
  416.     return ( (a >= 0) ? a/20 : (-round20up(-a)) );
  417.     }
  418. @//E*O*F wdbread.c//
  419. chmod u=rwx,g=rx,o=rx wdbread.c
  420.  
  421. echo x - wdbplot.c
  422. sed 's/^@//' > "wdbplot.c" <<'@//E*O*F wdbplot.c//'
  423. /* wdbplot -  draws wdb maps on unix plot files, using converted files */
  424.  
  425. /* written by Alan Paeth, University of Waterloo, January, 1984 */
  426.  
  427. /* wdbplot [clat] [clon] [latscale]  
  428.    draws a map centered at clat,clon degrees which spans a total of
  429.    lonscale degrees (all values are floats).
  430.  
  431.    NOTE: this version adjusts the vertical scale by secant(lat), ie, it
  432.    creates maps which are Mercator-like, especially at small scalings.
  433.    Here, the scale parameter (user specified) controls lon. (width) only.
  434.  
  435. SPECIAL NOTE: y coordinates are left-handed, ie, increase as one moves
  436.    downward from the northnmost edge of the map. Many display devices
  437.    work this way.
  438.    
  439.    wN - lat  (but will revert back to the proper)   lat - wS
  440.     
  441.    Color (based on rank) could also use some work
  442. */
  443.  
  444.  
  445. /* DEFINITIONS */
  446.  
  447. #define XRES     512    /* max output width in integer device units */
  448. #define YRES     512    /* max output height in integer device units */
  449.  
  450.  
  451. #include <stdio.h>
  452. #include <graphics/const.h>
  453. #include <graphics/ik_type.h>
  454. #include <graphics/ik_const.h>
  455.  
  456. #define PROGNAME "wdbplot"
  457. #define MASKREG 1
  458. #define DATANAME "map.data"
  459. #define INDEXNAME "map.index"
  460. #define READMODE 0 
  461. #define OPENFAIL -1
  462. #define UNITPERDEG 3600.0
  463.  
  464. #define GRIDCOLOR  8
  465. #define COASTCOLOR 1
  466. #define RIVERCOLOR 2
  467. #define BOUNDCOLOR 8
  468.  
  469.  
  470. static int lastcolor = 0;
  471. static int pengrabbed = 0;
  472. static int lastx = -1;
  473. static int lasty = -1;
  474.  
  475. double atof(), cos();
  476.  
  477. double xscale, yscale, clat, clon, xs, ys;
  478. double latgridwidth, latgridcount, latbase;
  479. double longridwidth, longridcount, lonbase;
  480. int mapcolor[64];
  481. int mapdata, mapindex;
  482.  
  483. short rank, count, shortn, shorts, shorte, shortw;
  484. long n, s, e, w;
  485. long i, totalrec;
  486. long lat, lon;
  487. char dlat, dlon;
  488. long wN, wS, wE, wW;
  489.  
  490.  
  491. /* ROUTINES */
  492.  
  493. double intpart(a)
  494. double a;
  495.     {
  496.     return( (double) ( (int) a ) );
  497.     }
  498.  
  499. lattopoint(latd)
  500. double latd;
  501.     {
  502.     return (   (double)  (  ( ((float) wN) - latd*UNITPERDEG) * ys)   );
  503.     }
  504.     
  505. lontopoint(lond)
  506. double lond;
  507.     {
  508.     return (   (double)  (   ( lond*UNITPERDEG - ((float) wW) ) * xs)   );
  509.     }
  510.     
  511. double pointtolat(y)
  512. int y;
  513.     {
  514.     return (  ( ((float) wN) - ((float) y) / ys) / UNITPERDEG );
  515.     }
  516.     
  517. double pointtolon(x)
  518. int x;
  519.     {
  520.     return ( ( ((float) x) / xs + ((float) wW) )/ UNITPERDEG );
  521.     }
  522.     
  523. abort(a,b)
  524.     char *a, *b;
  525.     {
  526.     fprintf(stderr,"%s error: ",PROGNAME);
  527.     fprintf(stderr,a,b);
  528.     fprintf(stderr,"\n");
  529.     exit(1);
  530.     }
  531.  
  532. /* THE MAIN */
  533.  
  534. main(argc,argv)
  535.     int argc;
  536.     char *argv[];
  537.     {
  538.     if (argc != 4)
  539.     abort("usage: [centerlat] [centerlon] [scale]", NULL);
  540.  
  541. /* open the data and index input, plotter for output */
  542.  
  543.     if ((mapdata = open(DATANAME, READMODE)) == OPENFAIL)
  544.         abort("open fail on %s", DATANAME);
  545.  
  546.     if ((mapindex = open(INDEXNAME, READMODE)) == OPENFAIL)
  547.         abort("open fail on %s", INDEXNAME);
  548.     
  549.     openplot();
  550.  
  551. /* now the initial sizing and reticle construction */
  552.  
  553.     clat =   atof(argv[1]);
  554.     clon =   atof(argv[2]);
  555.     xscale = atof(argv[3]);
  556.     
  557.     for (i =  0; i <  3; i++) mapcolor[i] = COASTCOLOR;
  558.     for (i =  3; i < 32; i++) mapcolor[i] = RIVERCOLOR;
  559.     for (i = 32; i < 48; i++) mapcolor[i] = BOUNDCOLOR;
  560.     for (i = 48; i < 64; i++) mapcolor[i] = BOUNDCOLOR;
  561.    
  562.     yscale = xscale * cos(clat*0.01745329252);
  563.  
  564.     wN = (long)( (clat + yscale/2.0) * UNITPERDEG );
  565.     wS = (long)( (clat - yscale/2.0) * UNITPERDEG );
  566.     wE = (long)( (clon + xscale/2.0) * UNITPERDEG );
  567.     wW = (long)( (clon - xscale/2.0) * UNITPERDEG );
  568.     ys = YRES / (float)(wN - wS);    /* 512 -> 0..512 on calls to draw */
  569.     xs = XRES / (float)(wE - wW);
  570.  
  571.     latgridcount = pointtolat(0) - pointtolat(511);
  572.     latgridwidth = 1.0;
  573.     while (latgridcount > 5.0)
  574.     {
  575.     latgridwidth *= 2.0;
  576.     latgridcount /= 2.0;
  577.     }
  578.     while (latgridcount < 3.0)
  579.     {
  580.     latgridwidth /= 2.0;
  581.     latgridcount *= 2.0;
  582.     }
  583.  
  584.     longridcount = pointtolon(511) - pointtolon(0);
  585.     longridwidth = 1.0;
  586.     while (longridcount > 5.0)
  587.     {
  588.     longridwidth *= 2.0;
  589.     longridcount /= 2.0;
  590.     }
  591.     while (longridcount < 3.0)
  592.     {
  593.     longridwidth /= 2.0;
  594.         longridcount *= 2.0;
  595.     }
  596.     
  597.     latbase = intpart( pointtolat(511)*latgridwidth ) / latgridwidth;
  598.     lonbase = intpart( pointtolon( 0 )*longridwidth ) / longridwidth;
  599.     
  600.     setcolor(GRIDCOLOR);
  601.     for (i = -1.0; i < latgridcount + 1.0; ++i)
  602.     {
  603.     moveto( 0.0   , (double) lattopoint(latbase + i*latgridwidth) );
  604.     drawto( 511.0 , (double) lattopoint(latbase + i*latgridwidth) );
  605.     }
  606.     for (i = -1.0; i < longridcount + 1.0; ++i)
  607.     {
  608.     moveto( (double) lontopoint(lonbase + i*longridwidth),   0.0 );
  609.     drawto( (double) lontopoint(lonbase + i*longridwidth), 511.0 );
  610.     }
  611.  
  612. /* now the top level count loop */
  613.  
  614.     while( read(mapindex, &count, 2) != 0)
  615.     {
  616.     read(mapindex, &rank, 2);
  617.     read(mapindex, &shortn, 2);
  618.     read(mapindex, &shorts, 2);
  619.     read(mapindex, &shorte, 2);
  620.     read(mapindex, &shortw, 2);
  621.     
  622.     setcolor( mapcolor[rank] );
  623.     
  624.     n = shortn * 20;
  625.     s = shorts * 20;
  626.     e = shorte * 20;
  627.     w = shortw * 20;
  628.     
  629.     if (s <= wN && n >= wS && w <= wE && e >= wW)
  630.         {
  631.         read(mapdata, &lat, 4);
  632.         read(mapdata, &lon, 4);
  633.         moveto( ((float)(lon - wW)) * xs, ((float)(lat - wS)) * ys);
  634.         for (i=0; i < (count - 1); i++)
  635.         {
  636.         read(mapdata, &dlat, 1);
  637.         read(mapdata, &dlon, 1);
  638.         lat += dlat;
  639.         lon += dlon;
  640.         drawto( ((float)(lon - wW)) * xs, ((float)(lat - wS)) * ys);
  641.         }
  642.         }
  643. /* skip forward from cur pos */
  644.     else lseek(mapdata, (long)(2*count+6), 1l);
  645.     }
  646.  
  647.     close(mapdata);
  648.     close(mapindex);
  649.     closeplot();
  650.     exit(0);
  651.     }
  652.  
  653. /*
  654.  * INTERFACE TO UNIX PLOT UTILITIES
  655.  */
  656.  
  657. openplot()
  658.     {
  659.     openpl();
  660.     erase();
  661.     space(0, 0, XRES, YRES);
  662.     }
  663.     
  664. closeplot()
  665.     {
  666.     closepl();
  667.     }
  668.     
  669. /* move and draw do simple clipping because grid is larger than viewport */
  670.  
  671. moveto(x, y)
  672.     double x, y;
  673.     {
  674.     int ix, iy;
  675.     ix = (int)(x+0.5);
  676.     iy = (int)(y+0.5);
  677.     if (ix < 0) ix = 0; else if (ix >= XRES) ix = XRES-1;
  678.     if (iy < 0) iy = 0; else if (iy >= YRES) iy = YRES-1;
  679.     if ((lastx != ix) || (lasty != iy)) move(ix, iy);
  680.     lastx = ix;
  681.     lasty = iy;
  682.     }
  683.  
  684. drawto(x, y)
  685.     double x, y;
  686.     {
  687.     int ix, iy;
  688. /* one could filter based on the present vector "color", if desired */
  689.     ix = (int)(x+0.5);
  690.     iy = (int)(y+0.5);
  691.     if (ix < 0) ix = 0; else if (ix >= XRES) ix = XRES-1;
  692.     if (iy < 0) iy = 0; else if (iy >= YRES) iy = YRES-1;
  693.     if ((lastx != ix) || (lasty != iy)) cont(ix, iy);
  694.     lastx = ix;
  695.     lasty = iy;
  696.     }
  697.  
  698. setcolor(c)
  699.     {
  700.     /* globalcol = c;    /*    /* one could record present vector color */
  701.     }
  702. @//E*O*F wdbplot.c//
  703. chmod u=rw,g=r,o=r wdbplot.c
  704.  
  705. exit 0
  706.