[Meshdb-cvs] src/geo/make-hgt-qt BSDmakefile,NONE,1.1 GNUmakefile,NONE,1.1 Makefile.common,NONE,1.1 main.c,NONE,1.1

[<le...@pr...>]

Update of /cvsroot/meshdb/src/geo/make-hgt-qt
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv3783

Added Files:
	BSDmakefile GNUmakefile Makefile.common main.c 
Log Message:
something quick to convert .hgt files to quadtree alt files

--- NEW FILE: BSDmakefile ---
.include "Makefile.common"
.include <bsd.prog.mk>

--- NEW FILE: GNUmakefile ---
include Makefile.common
include ../mk/GNUmakefile.prog

--- NEW FILE: Makefile.common ---
PROG=	make-hgt-qt
SRCS=	main.c 
NOMAN=
DEBUG=	-ggdb -O2 -Wall
LDADD=	${LDADD_GEO}
DPADD=	${LIBGEO}

#default: ${PROG}
#
#install:
#	-${INSTALL} -m ${NONBINMODE} alt.data alt.desc ${LIBDIR}/
#
# this requires the AUSLIG SG56 DEM data - it costs about AUD$110.
#alt:
#	./make-alt-qt /cdrom/SG56_DEM/sg56.grd \
#		alt.data alt.desc

--- NEW FILE: main.c ---
/*
 * Convert altitude data from SRTM 3deg .HGT files into a quadtree file.
 *
 * Assumes the filename is in the form SnnnEnnn.hgt, and consists
 * of a 1201x1201 square of network-endian 16-bit signed integers indicating
 * height above sea level in metres, and arranged positive east, and 
 * then positive north.
 *
 * David Leonard, 2004. Public domain.
 */

#include <stdio.h>
#include <string.h>
#include <errno.h>
#include <fcntl.h>
#include <err.h>
#include <sys/types.h>
#include <sys/mman.h>
#include <unistd.h>

#include "compat/compat.h"
#include "qt.h"

int
main(argc, argv)
	int argc;
	char *argv[];
{
	FILE *f, *desc;
	int qtfd;
	unsigned long w, h, i;
	int16_t nodata;
	float v, xll, yll, csz;
	int first;
	off_t mark;
	char *infile;
	char *qtdata;
	char *qtdesc;
	int qtw, qth;
	int x, y;
	struct qt *qt;
	int progress;
	int n;
	char *s, SN[2], EW[2];
	int sndeg, ewdeg;

	if (argc != 4) {
		fprintf(stderr, "usage: %s infile qt.data qt.desc\n", argv[0]);
		exit(1);
	}
	infile = argv[1];
	qtdata = argv[2];
	qtdesc = argv[3];

	s = strrchr(infile, '/');
	s = s == NULL ? infile : s + 1;
	if (sscanf(s, "%1[SN]%u%1[EW]%u.", SN, &sndeg, EW, &ewdeg) != 4)
	    errx(1, "%s: not in the format '*/[SN]num[EW]num.hgt'", s);
	xll = (*EW == 'E') ? ewdeg : -ewdeg;
	yll = (*SN == 'N') ? sndeg : -sndeg;
	w = 1201;
	h = 1201;
	csz = 1./1200;
	nodata = -32768;

	if ((f = fopen(infile, "r")) == NULL)
		err(1, "%s", infile);

	fprintf(stderr, "xll = %f, yll = %f\n", xll, yll);
	fprintf(stderr, "w x h = %lu x %lu\nnodata = %d\n", w, h, nodata);

	if ((qtfd = open(qtdata, O_TRUNC|O_CREAT|O_RDWR, 0666)) == -1)
		err(1, "%s", qtdata);
	if ((desc = fopen(qtdesc, "w")) == NULL)
		err(1, "%s", qtdesc);

	/* Find powers of 2 to hold the data set */
	qtw = qth = 1;
	while (qtw < w) qtw <<= 1;
	while (qth < h) qth <<= 1;
	fprintf(desc, "width %d\nheight %d\nelement_size %d\nfilename %s\n", 
		qtw, qth, sizeof (float), qtdata);
	fprintf(desc, "xll %f\nyll %f\ncellsz %f\n",
		xll + csz/2, 
		yll + csz/2,
		csz);
	fflush(desc);

	fprintf(stderr, "mapping %lu bytes\n", qtw * qth * sizeof (float));
	if (ftruncate(qtfd, qtw * qth * sizeof (float)) == -1)
		err(1, "%s: %lu bytes", qtdata, qtw * qth * sizeof (float));
	close(qtfd);

	fclose(desc);

	/* Now, re-open the file with writing */

	qt = qtopen(qtdesc, 1, NULL);

	fprintf(stderr, "filling with NaN\n");
	v = 0.0/0.0; /* NaN */
	for (x = 0; x < qtw * qth; x++)
		((float *)qt->data)[x] = v;
	
	/* Fill in the quadtree */
	fprintf(stderr, "filling in quadtree: ");
	progress = 0;

	for (y = 0; y < h; y++)  {
		if (y * 100 / h >= progress) {
			fprintf(stderr, "%d%% ", progress);
			progress  = 10 + y * 100 / h;
		}
		for (x = 0; x < w; x++)  {
			int16_t a,b;
			float *fp;
			fread(&a, sizeof a, 1, f);
			b = ntohs(a);
			fp = qtgetpoint_byindex(qt, x, y);
			if (b != nodata)
				*fp = b;
		}
	}
	fprintf(stderr, "\nsync\n");

	fclose(f);
	qtfree(qt);

	exit(0);
}