[Meshdb-cvs] src/geo/make-hgt-qt main.c,1.1,1.2

[<il...@pr...>]

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

Modified Files:
	main.c 
Log Message:
Merge from leonard-dev: working programme to import SRTM tiles to meshdb format.

Index: main.c
===================================================================
RCS file: /cvsroot/meshdb/src/geo/make-hgt-qt/main.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -d -r1.1 -r1.2
--- main.c	20 Sep 2004 13:37:19 -0000	1.1
+++ main.c	27 Dec 2004 04:09:10 -0000	1.2
@@ -11,6 +11,7 @@
 
 #include <stdio.h>
 #include <string.h>
+#include <stdlib.h>
 #include <errno.h>
 #include <fcntl.h>
 #include <err.h>
@@ -19,7 +20,9 @@
 #include <unistd.h>
 
 #include "compat/compat.h"
-#include "qt.h"
+#include "libgeo/qt.h"
+
+#define RES 1200
 
 int
 main(argc, argv)
@@ -30,8 +33,8 @@
 	int qtfd;
 	unsigned long w, h, i;
 	int16_t nodata;
-	float v, xll, yll, csz;
-	int first;
+	float v, csz;
+	int first, *ew, *sn;
 	off_t mark;
 	char *infile;
 	char *qtdata;
@@ -43,30 +46,49 @@
 	int n;
 	char *s, SN[2], EW[2];
 	int sndeg, ewdeg;
+	int ewmin, ewmax, snmin, snmax;
 
-	if (argc != 4) {
-		fprintf(stderr, "usage: %s infile qt.data qt.desc\n", argv[0]);
+	if (argc <= 3) {
+		fprintf(stderr, "usage: %s qt.data qt.desc file...\n", argv[0]);
 		exit(1);
 	}
-	infile = argv[1];
-	qtdata = argv[2];
-	qtdesc = argv[3];
+	qtdata = argv[1];
+	qtdesc = argv[2];
 
-	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;
+	ew = (int *)alloca((argc - 3) * sizeof (int));
+	sn = (int *)alloca((argc - 3) * sizeof (int));
 
-	if ((f = fopen(infile, "r")) == NULL)
-		err(1, "%s", infile);
+	/* Compute the total enclosing rectangle */
+	for (i = 3; i < argc; i++) {
+		infile = argv[i];
+		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);
+		if (*EW == 'W') ewdeg = -ewdeg;
+		if (*SN == 'S') sndeg = -sndeg;
+		fprintf(stderr, "%s -> %d,%d\n", infile, ewdeg, sndeg);
+		if (i == 3) {
+		    ewmin = ewmax = ewdeg;
+		    snmin = snmax = sndeg;
+		} else {
+		    if (ewdeg < ewmin) ewmin = ewdeg;
+		    if (ewdeg > ewmax) ewmax = ewdeg;
+		    if (sndeg < snmin) snmin = sndeg;
+		    if (sndeg > snmax) snmax = sndeg;
+		}
+		ew[i - 3] = ewdeg;
+		sn[i - 3] = sndeg;
+	}
+	ewmax += 1;
+	snmax += 1;
 
-	fprintf(stderr, "xll = %f, yll = %f\n", xll, yll);
+	w = RES * (ewmax - ewmin + 1) + 1;
+	h = RES * (snmax - snmin + 1) + 1;
+	csz = 1./RES;
+	nodata = -32768;
+
+	fprintf(stderr, "xll = %d, yll = %d\n", ewmin, snmin);
 	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)
@@ -81,8 +103,8 @@
 	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,
+		(float)ewmin + csz/2, 
+		(float)snmin + csz/2,
 		csz);
 	fflush(desc);
 
@@ -103,25 +125,40 @@
 		((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) {
+	for (i = 3; i < argc; i++) {
+	    int xo = (ew[i - 3] - ewmin) * RES;
+	    int yo = (sn[i - 3] - snmin) * RES;
+	    const unsigned long hh = RES + 1;
+	    const unsigned long ww = RES + 1;
+
+	    infile = argv[i];
+	    fprintf(stderr, "filling from %s @ %d,%d -> %d,%d:\n", 
+		infile, ew[i-3], sn[i-3], xo, yo);
+
+	    if ((f = fopen(infile, "r")) == NULL)
+		    err(1, "%s", infile);
+
+	    progress = 0;
+	    for (y = 0; y < hh; y++)  {
+		int yy = hh - y;
+		if (y * 100 / hh >= progress) {
 			fprintf(stderr, "%d%% ", progress);
-			progress  = 10 + y * 100 / h;
+			progress  = 10 + y * 100 / hh;
 		}
-		for (x = 0; x < w; x++)  {
+		for (x = 0; x < ww; x++)  {
 			int16_t a,b;
 			float *fp;
 			fread(&a, sizeof a, 1, f);
 			b = ntohs(a);
-			fp = qtgetpoint_byindex(qt, x, y);
+			fp = qtgetpoint_byindex(qt, x + xo, yo + yy);
 			if (b != nodata)
 				*fp = b;
 		}
+	    }
+	    fprintf(stderr, " ok\n");
 	}
-	fprintf(stderr, "\nsync\n");
+        fprintf(stderr, "sync\n");
 
 	fclose(f);
 	qtfree(qt);