v04i087: WDBII mapping software (10:1 compress, Unix PLOT)

Alan Wm Paeth awpaeth at watcgl.waterloo.edu
Tue Sep 20 11:23:24 AEST 1988


Posting-number: Volume 4, Issue 87
Submitted-by: "Alan Wm Paeth" <awpaeth at watcgl.waterloo.edu>
Archive-name: wdb-ii

I have had a number of requests for my World Data Bank II mapping
software (C Unix). For sites that have all five mag tape volumes of
this dataset (don't come crying to me!), this software package provides
two C programs. The first does a 10:1 compression on the datasets and
additionally maps the straight vector detail into an index file of offsets,
plus a companion vector file (two bytes per vector). The second program
draws a rectangular projection of the dataset for specified lat,long,scale
and produces Unix PLOT output. It can be very easily upgraded. Enjoy.

   /Alan Paeth
   Computer Graphics Laboratory
   University of Waterloo

-----------------------------------------------------------------------------

# This is a shell archive.  Remove anything before this line,
# then unpack it by saving it in a file and typing "sh file".
#
# Wrapped by watcgl!awpaeth on Sat Sep 17 18:28:33 EDT 1988
# Contents:  README Makefile wdbread.c wdbplot.c
 
echo x - README
sed 's/^@//' > "README" <<'@//E*O*F README//'
README    -- this file

MAKEFILE  -- compile both sources

wdbread.c -- (see source code) converts an EBCDIC WDBII dataset (20 char recs)
	     into two output files: map.index and map.data, while producing
	     a large amount of line and summary detail. The first file provides
	     bounding box and # of vertex information into the second file
	     (map.data), this latter file contains (dx,dy) byte pairs. wdbread
	     will fragment long vectors as needed to provide this format, and
	     will report this on the log it produces on stdout.

	     Multiple runs of this program will place (overwrite, if you're
	     not careful) map/index pairs for the 16 or so volumes which
	     comprise all of WDBII (arranged by continents and detail).
	     All .index and .data files may then be "cat"ed into one giant
	     file which represents the entire world. This is possible by
	     design: the skip offsets in the .index file are relative byte
	     counts in each companion .data file, so each file can be
	     appended to in tandem.
	     
wdbplot.c -- given location and scale produces UNIX PLOT output. The last
	     section of source code provides a simple means to customize
	     the software to direct drive other devices, although we most
	     often use Tektronix output (via plot). Other suggested fixes
	     would be to "trigger" based on vector rank, which changes as
	     a function of political boundary, river, coastline, etc. More
	     ambitious users might wish to encode more general map projections:
	     Mercator and Rectangular equal area would be the most simple
	     to employ because the clipping borders would remain roughly
	     unchanged.

Author's disclamer: I have no intent in modifying or upgrading this ancient
code: if you can make big bucks on it, go for it. If you want to discuss
digital cartography from a less nuts and bolts level, drop me a line sometime!

-- Alan Wm Paeth (awpaeth at watcgl.waterloo.edu, awpaeth at watcgl.waterloo.cdn).
@//E*O*F README//
chmod u=rw,g=r,o=r README
 
echo x - Makefile
sed 's/^@//' > "Makefile" <<'@//E*O*F Makefile//'
CFLAGS = -g


wdbplot: wdbplot.c
	cc $(CFLAGS) wdbplot.c -lplot -lm -o $@

wdbread: wdbread.c
	cc $(CFLAGS) wdbread.c -o $@
@//E*O*F Makefile//
chmod u=rw,g=r,o=r Makefile
 
echo x - wdbread.c
sed 's/^@//' > "wdbread.c" <<'@//E*O*F wdbread.c//'
/* wdbread - converts raw world data bank II files into vax words

/* NOTE: like old wdbiiread, but limits consec. points to +-127 units */

/* written by Alan Paeth, University of Waterloo, January, 1984 */

/* wdbread <in >out 
   reads in hunks of twenty byte EBCDIC records, as per the FORTRAN 
   spec of WDBII, and writes them out as signed 8-bit integers (signed
   bytes), as delta values. No conversion of characters is done for other
   than the digits and few other expected sparse characters.

NOTE: the output is named map.data and map.index, both placed on the
current directory. The file >out contains a summary of processed records.
*/

/* DEFINITIONS */

#include <stdio.h>

#define PROGNAME "wdbread"
#define EBCDICFLAG 0		/* 0=ascii input, 1=ebcdic */
#define INBUFSIZE 20
#define PIPEOUT 1
#define DATANAME "map.data"
#define INDEXNAME "map.index"
#define WRITEMODE "w"
#define OPENFAIL 0 
#define LON360 1296000
#define FRAGLIMIT 20
#define maxint  2147483647
#define minint -2147483648 
#define v(arg) (inbuf[arg]-'0')

#define ABS(a) ((a)<0?(-a):(a))

FILE *mapdata, *mapindex;
char inbuf[INBUFSIZE+1], remap[256];
int rank, count;
int fragcount, totfrag, readrec, writerec, readpoint, writepoint;
int i, j, k, n, s, e, w;
short shortn, shorts, shorte, shortw;
int olat, olon, privlat, privlon, dlat, dlon, oprivlat, oprivlon;

/* THE MAIN */

main(argc,argv)
    int  argc;
    char **argv;
    {

/* open the data and index output */

    if ((mapdata = fopen(DATANAME, WRITEMODE)) == OPENFAIL)
        abort("open fail on %s", DATANAME);

    if ((mapindex = fopen(INDEXNAME, WRITEMODE)) == OPENFAIL)
        abort("open fail on %s", INDEXNAME);
    
/* set remap to map into ASCII, digit 0 for all other characters */

    if (EBCDICFLAG) 
/* EBCDIC REMAPPING STUFF */
        {
	for (i=0; i<256; i++) remap[i] = '0';		/* def. 0 */
	for (i=0; i<10; i++) remap[240+i] = 48+i;	/* digits */
	remap[103] = 'W';				/* w,W,s,S */
	remap[231] = 'W';
	remap[99] = 'S';
	remap[227] = 'S';
	}
    else
/* ASCII REMAPPING (identity xform, space -> digit 0 */
	{
	for (i=0; i<256; i++) remap[i] = i;		/* identity */
	remap[' '] = '0';				/* space -> 0 */
	}

    inbuf[INBUFSIZE] = '\0';
    readrec = 0;
    writerec = 0;
    totfrag = 0;

/* now the top level count loop */

    while( fread(&inbuf[0], 1, INBUFSIZE, stdin) == INBUFSIZE)
	{
	readrec++;
	readpoint = 0;
	writepoint = 0;
	
	for (i=0; i<INBUFSIZE; i++) inbuf[i] = remap[inbuf[i]];
	printf("/%s/\n", inbuf);
	sscanf( &inbuf[9], "%6ld", &count);
	sscanf( &inbuf[7], "%2ld", &rank);
	printf("count %d rank %d\n", count, rank);
	fflush(stdout);
	
	n = minint;
	s = maxint;
	e = minint;
	w = maxint;

/* NEW STUFF */

	fragcount = 0;
	
/* the high-speed data loop */
	for (j=0; j<count; j++)
	    {
/* optimizations */
	    register int lat, lon, rSIX, rTEN;
	    register char *base;
	    rTEN = 10;
	    rSIX = 6;

	    olat = lat;
	    olon = lon;
	    
    	    fread(&inbuf[0], 1, INBUFSIZE, stdin);
	    for (i=0; i<INBUFSIZE; i++) inbuf[i] = remap[inbuf[i]];
	    
	    readpoint++;
/*
 * optimization of the following:
 *
 *	    lat = v(0)*36000 + v(1)*3600 + v(2)*600 +
 *	      v(3)*60 + v(4)*10 + v(5);
 *	    if (inbuf[6] == 'S') lat = -lat;
 */
	    base = &inbuf[0];
	    lat = *base++;
	    lat *= rTEN;
	    lat += *base++;
	    lat *= rSIX;
	    lat += *base++;
	    lat *= rTEN;
	    lat += *base++;
	    lat *= rSIX;
	    lat += *base++;
	    lat *= rTEN;
	    lat += *base++;
	    lat -= '0'*(36000+3600+600+60+10+1); /* correct for ascii 000000 */
	    if (*base == 'S') lat = -lat;
/*
 * optimization of the following:
 *
 *	    lon = v(7)*360000 + v(8)*36000 + v(9)*3600 +
 *	      v(10)*600 + v(11)*60 + v(12)*10 + v(13);
 *  	    if (inbuf[14] == 'W') lon = -lon;
 */
	    base = &inbuf[7];
	    lon = *base++;
	    lon *= rTEN;
	    lon += *base++;
	    lon *= rTEN;
	    lon += *base++;
	    lon *= rSIX;
	    lon += *base++;
	    lon *= rTEN;
	    lon += *base++;
	    lon *= rSIX;
	    lon += *base++;
	    lon *= rTEN;
	    lon += *base++;
	    lon -= '0'*(360000+36000+3600+600+60+10+1); /* ascii 0000000 */
	    if (*base == 'W') lon = -lon;
/*
 * hemisphere check
 */
	    if (lon >= LON360) lon -= LON360;
	    if ((j != 0) && (ABS(olon-lon)>(LON360/2)))
		{
		warning("date-line crossing: %d\n", lon);
		}
/*
 * adjust bounding box
 */
	    if (lat > n) n = lat;
	    if (lat < s) s = lat;
	    if (lon > e) e = lon;
	    if (lon < w) w = lon;
	    
    	    if (j == 0)
		{
		dlat = lat;
		dlon = lon;
		fwrite(&dlat, 1, 4, mapdata);
		fwrite(&dlon, 1, 4, mapdata);
		writepoint++;
		}
	    else
		{
		dlat = lat - olat;
	        dlon = lon - olon;
	        fragcount = roundupdiv( lmax( labs(dlat), labs(dlon) ) , 127l);
		if (fragcount ==1)
		    {
		    fwrite(&dlat, 1, 1, mapdata);
		    fwrite(&dlon, 1, 1, mapdata);
		    writepoint++;
		    }
		else
		    {
		    if (fragcount == 0)
		        warning("fragcount=0 => consec. duplicate points",NULL);
		    if (fragcount > FRAGLIMIT)
		        {
			warning("excessive fragments %ld",fragcount);
			warning("occurring at segment: %d\n", j);
			if (fragcount > 2000)
			    {
			    warning("trace vars: \n");
			    warning(" lat = %d\n",  lat);
			    warning("olat = %d\n", olat);
			    warning("dlat = %d\n", dlat);
			    warning(" lon = %d\n",  lon);
			    warning("olon = %d\n", olon);
			    warning("dlon = %d\n", dlon);
			    }
			}
		    oprivlat = 0;
		    oprivlon = 0;
		    for (k = 0; k < fragcount; k++)
		        {
		        privlat = rounddiv( dlat*(k+1), fragcount) - oprivlat;
		        privlon = rounddiv( dlon*(k+1), fragcount) - oprivlon;
		        oprivlat = privlat;
		        oprivlon = privlon;

		        fwrite( &privlat, 1, 1, mapdata);
		        fwrite( &privlon, 1, 1, mapdata);

		        writepoint++;
		        }
		    printf(" frg=%ld dx=%ld dy=%ld\n", fragcount,dlat,dlon);
		    totfrag += (fragcount - 1);
		    }
	        } 
	    }
    	fwrite( &writepoint, 1, 2, mapindex);
    	fwrite( &rank, 1, 2, mapindex);
	writerec++;

/* round n,e +, round s,w - */

	shortn = round20up(n);
	shorts = round20down(s);
	shorte = round20up(e);
	shortw = round20down(w);

	fwrite( &shortn, 1, 2, mapindex);
    	fwrite( &shorts, 1, 2, mapindex);
    	fwrite( &shorte, 1, 2, mapindex);
    	fwrite( &shortw, 1, 2, mapindex);

	printf("read=%ld write=%ld n=%ld s=%ld e=%ld w=%ld\n\n", readpoint,
	   writepoint, n, s, e, w);
	fflush(stdout);
	}
    printf("\n\nTOTALS: records in=%ld records out=%ld new fragment=%ld\n",
		readrec,writerec,totfrag);
    fclose(mapdata);
    fclose(mapindex);
    exit(0);
    }

/* ROUTINES */

abort(a,b)
    char *a, *b;
    {
    fprintf(stderr,"%s error: ",PROGNAME);
    fprintf(stderr,a,b);
    fprintf(stderr,"\n");
    exit(1);
    }

warning(a,b)
    char *a, *b;
    {
    printf("\nWARNING: ");
    printf(a,b);
    printf("\n");
    fflush(stdout);
    }

roundupdiv(a,b)
    int a,b;
    {
    return (a >= 0) ? ((a - 1) / b + 1) : -roundupdiv(-a,b);
    }

rounddiv(a,b)
    int a,b;
    {
    return (a >= 0) ? (2*a + b) / (2*b) : -rounddiv(-a,b);
    }

lmax(a,b)
    int a,b;
    {
    return ((a >= b) ? a : b);
    }

labs(a)
    int a;
    {
    return ((a >= 0) ? a : -a);
    }

round20up(a)
    int a;
    {
    return ( (a >= 0) ? (a+19)/20 : (-round20down(-a)) );
    }

round20down(a)
    int a;
    {
    return ( (a >= 0) ? a/20 : (-round20up(-a)) );
    }
@//E*O*F wdbread.c//
chmod u=rwx,g=rx,o=rx wdbread.c
 
echo x - wdbplot.c
sed 's/^@//' > "wdbplot.c" <<'@//E*O*F wdbplot.c//'
/* wdbplot -  draws wdb maps on unix plot files, using converted files */

/* written by Alan Paeth, University of Waterloo, January, 1984 */

/* wdbplot [clat] [clon] [latscale]  
   draws a map centered at clat,clon degrees which spans a total of
   lonscale degrees (all values are floats).

   NOTE: this version adjusts the vertical scale by secant(lat), ie, it
   creates maps which are Mercator-like, especially at small scalings.
   Here, the scale parameter (user specified) controls lon. (width) only.

SPECIAL NOTE: y coordinates are left-handed, ie, increase as one moves
   downward from the northnmost edge of the map. Many display devices
   work this way.
   
   wN - lat  (but will revert back to the proper)   lat - wS
    
   Color (based on rank) could also use some work
*/


/* DEFINITIONS */

#define XRES	 512	/* max output width in integer device units */
#define YRES	 512	/* max output height in integer device units */


#include <stdio.h>
#include <graphics/const.h>
#include <graphics/ik_type.h>
#include <graphics/ik_const.h>

#define PROGNAME "wdbplot"
#define MASKREG 1
#define DATANAME "map.data"
#define INDEXNAME "map.index"
#define READMODE 0 
#define OPENFAIL -1
#define UNITPERDEG 3600.0

#define GRIDCOLOR  8
#define COASTCOLOR 1
#define RIVERCOLOR 2
#define BOUNDCOLOR 8


static int lastcolor = 0;
static int pengrabbed = 0;
static int lastx = -1;
static int lasty = -1;

double atof(), cos();

double xscale, yscale, clat, clon, xs, ys;
double latgridwidth, latgridcount, latbase;
double longridwidth, longridcount, lonbase;
int mapcolor[64];
int mapdata, mapindex;

short rank, count, shortn, shorts, shorte, shortw;
long n, s, e, w;
long i, totalrec;
long lat, lon;
char dlat, dlon;
long wN, wS, wE, wW;


/* ROUTINES */

double intpart(a)
double a;
    {
    return( (double) ( (int) a ) );
    }

lattopoint(latd)
double latd;
    {
    return (   (double)  (  ( ((float) wN) - latd*UNITPERDEG) * ys)   );
    }
    
lontopoint(lond)
double lond;
    {
    return (   (double)  (   ( lond*UNITPERDEG - ((float) wW) ) * xs)   );
    }
    
double pointtolat(y)
int y;
    {
    return (  ( ((float) wN) - ((float) y) / ys) / UNITPERDEG );
    }
    
double pointtolon(x)
int x;
    {
    return ( ( ((float) x) / xs + ((float) wW) )/ UNITPERDEG );
    }
    
abort(a,b)
    char *a, *b;
    {
    fprintf(stderr,"%s error: ",PROGNAME);
    fprintf(stderr,a,b);
    fprintf(stderr,"\n");
    exit(1);
    }

/* THE MAIN */

main(argc,argv)
    int argc;
    char *argv[];
    {
    if (argc != 4)
	abort("usage: [centerlat] [centerlon] [scale]", NULL);

/* open the data and index input, plotter for output */

    if ((mapdata = open(DATANAME, READMODE)) == OPENFAIL)
        abort("open fail on %s", DATANAME);

    if ((mapindex = open(INDEXNAME, READMODE)) == OPENFAIL)
        abort("open fail on %s", INDEXNAME);
    
    openplot();

/* now the initial sizing and reticle construction */

    clat =   atof(argv[1]);
    clon =   atof(argv[2]);
    xscale = atof(argv[3]);
    
    for (i =  0; i <  3; i++) mapcolor[i] = COASTCOLOR;
    for (i =  3; i < 32; i++) mapcolor[i] = RIVERCOLOR;
    for (i = 32; i < 48; i++) mapcolor[i] = BOUNDCOLOR;
    for (i = 48; i < 64; i++) mapcolor[i] = BOUNDCOLOR;
   
    yscale = xscale * cos(clat*0.01745329252);

    wN = (long)( (clat + yscale/2.0) * UNITPERDEG );
    wS = (long)( (clat - yscale/2.0) * UNITPERDEG );
    wE = (long)( (clon + xscale/2.0) * UNITPERDEG );
    wW = (long)( (clon - xscale/2.0) * UNITPERDEG );
    ys = YRES / (float)(wN - wS);	/* 512 -> 0..512 on calls to draw */
    xs = XRES / (float)(wE - wW);

    latgridcount = pointtolat(0) - pointtolat(511);
    latgridwidth = 1.0;
    while (latgridcount > 5.0)
	{
	latgridwidth *= 2.0;
	latgridcount /= 2.0;
	}
    while (latgridcount < 3.0)
	{
	latgridwidth /= 2.0;
	latgridcount *= 2.0;
	}

    longridcount = pointtolon(511) - pointtolon(0);
    longridwidth = 1.0;
    while (longridcount > 5.0)
	{
	longridwidth *= 2.0;
	longridcount /= 2.0;
	}
    while (longridcount < 3.0)
	{
	longridwidth /= 2.0;
        longridcount *= 2.0;
	}
	
    latbase = intpart( pointtolat(511)*latgridwidth ) / latgridwidth;
    lonbase = intpart( pointtolon( 0 )*longridwidth ) / longridwidth;
    
    setcolor(GRIDCOLOR);
    for (i = -1.0; i < latgridcount + 1.0; ++i)
	{
	moveto( 0.0   , (double) lattopoint(latbase + i*latgridwidth) );
	drawto( 511.0 , (double) lattopoint(latbase + i*latgridwidth) );
	}
    for (i = -1.0; i < longridcount + 1.0; ++i)
	{
	moveto( (double) lontopoint(lonbase + i*longridwidth),   0.0 );
	drawto( (double) lontopoint(lonbase + i*longridwidth), 511.0 );
	}

/* now the top level count loop */

    while( read(mapindex, &count, 2) != 0)
	{
	read(mapindex, &rank, 2);
	read(mapindex, &shortn, 2);
	read(mapindex, &shorts, 2);
	read(mapindex, &shorte, 2);
	read(mapindex, &shortw, 2);
	
	setcolor( mapcolor[rank] );
	
	n = shortn * 20;
	s = shorts * 20;
	e = shorte * 20;
	w = shortw * 20;
	
	if (s <= wN && n >= wS && w <= wE && e >= wW)
	    {
	    read(mapdata, &lat, 4);
	    read(mapdata, &lon, 4);
	    moveto( ((float)(lon - wW)) * xs, ((float)(lat - wS)) * ys);
	    for (i=0; i < (count - 1); i++)
		{
		read(mapdata, &dlat, 1);
		read(mapdata, &dlon, 1);
		lat += dlat;
		lon += dlon;
		drawto( ((float)(lon - wW)) * xs, ((float)(lat - wS)) * ys);
		}
	    }
/* skip forward from cur pos */
	else lseek(mapdata, (long)(2*count+6), 1l);
	}

    close(mapdata);
    close(mapindex);
    closeplot();
    exit(0);
    }

/*
 * INTERFACE TO UNIX PLOT UTILITIES
 */

openplot()
    {
    openpl();
    erase();
    space(0, 0, XRES, YRES);
    }
    
closeplot()
    {
    closepl();
    }
    
/* move and draw do simple clipping because grid is larger than viewport */

moveto(x, y)
    double x, y;
    {
    int ix, iy;
    ix = (int)(x+0.5);
    iy = (int)(y+0.5);
    if (ix < 0) ix = 0; else if (ix >= XRES) ix = XRES-1;
    if (iy < 0) iy = 0; else if (iy >= YRES) iy = YRES-1;
    if ((lastx != ix) || (lasty != iy)) move(ix, iy);
    lastx = ix;
    lasty = iy;
    }

drawto(x, y)
    double x, y;
    {
    int ix, iy;
/* one could filter based on the present vector "color", if desired */
    ix = (int)(x+0.5);
    iy = (int)(y+0.5);
    if (ix < 0) ix = 0; else if (ix >= XRES) ix = XRES-1;
    if (iy < 0) iy = 0; else if (iy >= YRES) iy = YRES-1;
    if ((lastx != ix) || (lasty != iy)) cont(ix, iy);
    lastx = ix;
    lasty = iy;
    }

setcolor(c)
    {
    /* globalcol = c;	/*	/* one could record present vector color */
    }
@//E*O*F wdbplot.c//
chmod u=rw,g=r,o=r wdbplot.c
 
exit 0



More information about the Comp.sources.misc mailing list