awpaeth@watcgl.waterloo.edu (Alan Wm Paeth) (09/20/88)
Posting-number: Volume 4, Issue 87
Submitted-by: "Alan Wm Paeth" <awpaeth@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@watcgl.waterloo.edu, awpaeth@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