nad2bin products NTv2, nad2nad removed

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@2118 4e78687f-474d-0410-85f9-8d5e500ac6b2

6 files changed, 259 insertions, 73 deletions

diff --git a/ChangeLog b/ChangeLog
index 06b20525..42c14102 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,7 @@
 2011-11-13  Frank Warmerdam  <warmerdam@pobox.com>
 
+	* src/nad2bin.c: Modified to write NTv2 format files.
+
 	* src/pj_init.c: avoid casting warning with old_locale.
 
 2011-09-28  Frank Warmerdam  <warmerdam@pobox.com>
diff --git a/man/man1/Makefile.am b/man/man1/Makefile.am
index 46b4fe4c..b175aaa0 100644
--- a/man/man1/Makefile.am
+++ b/man/man1/Makefile.am
@@ -1,4 +1,4 @@
-man_MANS	=	geod.1 nad2nad.1 proj.1 cs2cs.1
+man_MANS	=	geod.1 proj.1 cs2cs.1
 
 EXTRA_DIST	=	$(man_MANS)
 
diff --git a/man/man1/Makefile.in b/man/man1/Makefile.in
index 6ef17879..e04701f4 100644
--- a/man/man1/Makefile.in
+++ b/man/man1/Makefile.in
@@ -187,7 +187,7 @@ target_alias = @target_alias@
 top_build_prefix = @top_build_prefix@
 top_builddir = @top_builddir@
 top_srcdir = @top_srcdir@
-man_MANS = geod.1 nad2nad.1 proj.1 cs2cs.1
+man_MANS = geod.1 proj.1 cs2cs.1
 EXTRA_DIST = $(man_MANS)
 all: all-am
 
diff --git a/src/Makefile.am b/src/Makefile.am
index 7a6b356d..9d57053d 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,4 +1,4 @@
-bin_PROGRAMS =	proj nad2nad nad2bin geod cs2cs
+bin_PROGRAMS =	proj nad2bin geod cs2cs
 
 INCLUDES =	-DPROJ_LIB=\"$(pkgdatadir)\" \
 		-DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@
@@ -9,13 +9,11 @@ EXTRA_DIST = makefile.vc proj.def
 
 proj_SOURCES = proj.c gen_cheb.c p_series.c
 cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c
-nad2nad_SOURCES = nad2nad.c 
 nad2bin_SOURCES = nad2bin.c
 geod_SOURCES = geod.c geod_set.c geod_for.c geod_inv.c geodesic.h
 
 proj_LDADD = libproj.la
 cs2cs_LDADD = libproj.la
-nad2nad_LDADD = libproj.la
 nad2bin_LDADD = libproj.la
 geod_LDADD = libproj.la
 
diff --git a/src/Makefile.in b/src/Makefile.in
index 5e02a52c..83fceb88 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -36,8 +36,8 @@ PRE_UNINSTALL = :
 POST_UNINSTALL = :
 build_triplet = @build@
 host_triplet = @host@
-bin_PROGRAMS = proj$(EXEEXT) nad2nad$(EXEEXT) nad2bin$(EXEEXT) \
-	geod$(EXEEXT) cs2cs$(EXEEXT)
+bin_PROGRAMS = proj$(EXEEXT) nad2bin$(EXEEXT) geod$(EXEEXT) \
+	cs2cs$(EXEEXT)
 subdir = src
 DIST_COMMON = $(include_HEADERS) $(srcdir)/Makefile.am \
 	$(srcdir)/Makefile.in $(srcdir)/proj_config.h.in
@@ -121,9 +121,6 @@ geod_DEPENDENCIES = libproj.la
 am_nad2bin_OBJECTS = nad2bin.$(OBJEXT)
 nad2bin_OBJECTS = $(am_nad2bin_OBJECTS)
 nad2bin_DEPENDENCIES = libproj.la
-am_nad2nad_OBJECTS = nad2nad.$(OBJEXT)
-nad2nad_OBJECTS = $(am_nad2nad_OBJECTS)
-nad2nad_DEPENDENCIES = libproj.la
 am_proj_OBJECTS = proj.$(OBJEXT) gen_cheb.$(OBJEXT) p_series.$(OBJEXT)
 proj_OBJECTS = $(am_proj_OBJECTS)
 proj_DEPENDENCIES = libproj.la
@@ -141,9 +138,9 @@ LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \
 	--mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \
 	$(LDFLAGS) -o $@
 SOURCES = $(libproj_la_SOURCES) $(cs2cs_SOURCES) $(geod_SOURCES) \
-	$(nad2bin_SOURCES) $(nad2nad_SOURCES) $(proj_SOURCES)
+	$(nad2bin_SOURCES) $(proj_SOURCES)
 DIST_SOURCES = $(libproj_la_SOURCES) $(cs2cs_SOURCES) $(geod_SOURCES) \
-	$(nad2bin_SOURCES) $(nad2nad_SOURCES) $(proj_SOURCES)
+	$(nad2bin_SOURCES) $(proj_SOURCES)
 HEADERS = $(include_HEADERS)
 ETAGS = etags
 CTAGS = ctags
@@ -271,12 +268,10 @@ include_HEADERS = proj_api.h org_proj4_Projections.h
 EXTRA_DIST = makefile.vc proj.def
 proj_SOURCES = proj.c gen_cheb.c p_series.c
 cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c
-nad2nad_SOURCES = nad2nad.c 
 nad2bin_SOURCES = nad2bin.c
 geod_SOURCES = geod.c geod_set.c geod_for.c geod_inv.c geodesic.h
 proj_LDADD = libproj.la
 cs2cs_LDADD = libproj.la
-nad2nad_LDADD = libproj.la
 nad2bin_LDADD = libproj.la
 geod_LDADD = libproj.la
 lib_LTLIBRARIES = libproj.la
@@ -457,9 +452,6 @@ geod$(EXEEXT): $(geod_OBJECTS) $(geod_DEPENDENCIES)
 nad2bin$(EXEEXT): $(nad2bin_OBJECTS) $(nad2bin_DEPENDENCIES) 
 	@rm -f nad2bin$(EXEEXT)
 	$(LINK) $(nad2bin_OBJECTS) $(nad2bin_LDADD) $(LIBS)
-nad2nad$(EXEEXT): $(nad2nad_OBJECTS) $(nad2nad_DEPENDENCIES) 
-	@rm -f nad2nad$(EXEEXT)
-	$(LINK) $(nad2nad_OBJECTS) $(nad2nad_LDADD) $(LIBS)
 proj$(EXEEXT): $(proj_OBJECTS) $(proj_DEPENDENCIES) 
 	@rm -f proj$(EXEEXT)
 	$(LINK) $(proj_OBJECTS) $(proj_LDADD) $(LIBS)
@@ -577,7 +569,6 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/jniproj.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/mk_cheby.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/nad2bin.Po@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/nad2nad.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/nad_cvt.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/nad_init.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/nad_intr.Plo@am__quote@
diff --git a/src/nad2bin.c b/src/nad2bin.c
index db4bb076..17927f75 100644
--- a/src/nad2bin.c
+++ b/src/nad2bin.c
@@ -1,65 +1,260 @@
-/* Convert bivariate ASCII NAD27 to NAD83 tables to binary structure */
+/* Convert bivariate ASCII NAD27 to NAD83 tables to NTv2 binary structure */
 #include <stdio.h>
 #include <stdlib.h>
+#include <assert.h>
+
 #define PJ_LIB__
 #include <projects.h>
 #define U_SEC_TO_RAD 4.848136811095359935899141023e-12
 	static char
 *usage = "<ASCII_dist_table local_bin_table";
 
+/************************************************************************/
+/*                             swap_words()                             */
+/*                                                                      */
+/*      Convert the byte order of the given word(s) in place.           */
+/************************************************************************/
+
+static int  byte_order_test = 1;
+#define IS_LSB	(((unsigned char *) (&byte_order_test))[0] == 1)
+
+static void swap_words( unsigned char *data, int word_size, int word_count )
+
+{
+    int	word;
+
+    for( word = 0; word < word_count; word++ )
+    {
+        int	i;
+        
+        for( i = 0; i < word_size/2; i++ )
+        {
+            int	t;
+            
+            t = data[i];
+            data[i] = data[word_size-i-1];
+            data[word_size-i-1] = t;
+        }
+        
+        data += word_size;
+    }
+}
+
+/************************************************************************/
+/*                                main()                                */
+/************************************************************************/
 int main(int argc, char **argv) {
-	struct CTABLE ct;
-	FLP *p, t;
-	size_t tsize;
-	int i, j, ichk;
-	long lam, laml, phi, phil;
-	FILE *bin;
-
-	if (argc != 2) {
-		fprintf(stderr,"usage: %s %s\n", argv[0], usage);
-		exit(1);
-	}
-	fgets(ct.id, MAX_TAB_ID, stdin);
-	scanf("%d %d %*d %lf %lf %lf %lf", &ct.lim.lam, &ct.lim.phi,
-		&ct.ll.lam, &ct.del.lam, &ct.ll.phi, &ct.del.phi);
-	if (!(ct.cvs = (FLP *)malloc(tsize = ct.lim.lam * ct.lim.phi *
-		sizeof(FLP)))) {
-		perror("mem. alloc");
-		exit(1);
-	}
-	ct.ll.lam *= DEG_TO_RAD;
-	ct.ll.phi *= DEG_TO_RAD;
-	ct.del.lam *= DEG_TO_RAD;
-	ct.del.phi *= DEG_TO_RAD;
-	/* load table */
-	for (p = ct.cvs, i = 0; i < ct.lim.phi; ++i) {
-		scanf("%d:%ld %ld", &ichk, &laml, &phil);
-		if (ichk != i) {
-			fprintf(stderr,"format check on row\n");
-			exit(1);
-		}
-		t.lam = laml * U_SEC_TO_RAD;
-		t.phi = phil * U_SEC_TO_RAD;
-		*p++ = t;
-		for (j = 1; j < ct.lim.lam; ++j) {
-			scanf("%ld %ld", &lam, &phi);
-			t.lam = (laml += lam) * U_SEC_TO_RAD;
-			t.phi = (phil += phi) * U_SEC_TO_RAD;
-			*p++ = t;
-		}
-	}
-	if (feof(stdin)) {
-		fprintf(stderr, "premature EOF\n");
-		exit(1);
-	}
-	if (!(bin = freopen(argv[1], "wb", stdout))) {
-		perror(argv[1]);
-		exit(2);
-	}
-	if (fwrite(&ct, sizeof(ct), 1, stdout) != 1 ||
-		fwrite(ct.cvs, tsize, 1, stdout) != 1) {
-		fprintf(stderr, "output failure\n");
-		exit(2);
-	}
-	exit(0); /* normal completion */
+    struct CTABLE ct;
+    FLP *p, t;
+    size_t tsize;
+    int i, j, ichk;
+    long lam, laml, phi, phil;
+    FILE *fp;
+
+    const char *GS_TYPE  = "SECONDS";
+    const char *VERSION  = "";
+    const char *SYSTEM_F = "NAD27";
+    const char *SYSTEM_T = "NAD83";
+    const char *SUB_NAME = "";
+    const char *CREATED  = "";
+    const char *UPDATED  = "";
+
+    if (argc != 2) {
+        fprintf(stderr,"usage: %s %s\n", argv[0], usage);
+        exit(1);
+    }
+
+/* ==================================================================== */
+/*      Read the ASCII Table                                            */
+/* ==================================================================== */
+
+    fgets(ct.id, MAX_TAB_ID, stdin);
+    scanf("%d %d %*d %lf %lf %lf %lf", &ct.lim.lam, &ct.lim.phi,
+          &ct.ll.lam, &ct.del.lam, &ct.ll.phi, &ct.del.phi);
+    if (!(ct.cvs = (FLP *)malloc(tsize = ct.lim.lam * ct.lim.phi *
+                                 sizeof(FLP)))) {
+        perror("mem. alloc");
+        exit(1);
+    }
+    ct.ll.lam *= DEG_TO_RAD;
+    ct.ll.phi *= DEG_TO_RAD;
+    ct.del.lam *= DEG_TO_RAD;
+    ct.del.phi *= DEG_TO_RAD;
+    /* load table */
+    for (p = ct.cvs, i = 0; i < ct.lim.phi; ++i) {
+        scanf("%d:%ld %ld", &ichk, &laml, &phil);
+        if (ichk != i) {
+            fprintf(stderr,"format check on row\n");
+            exit(1);
+        }
+        t.lam = laml * U_SEC_TO_RAD;
+        t.phi = phil * U_SEC_TO_RAD;
+        *p++ = t;
+        for (j = 1; j < ct.lim.lam; ++j) {
+            scanf("%ld %ld", &lam, &phi);
+            t.lam = (laml += lam) * U_SEC_TO_RAD;
+            t.phi = (phil += phi) * U_SEC_TO_RAD;
+            *p++ = t;
+        }
+    }
+    if (feof(stdin)) {
+        fprintf(stderr, "premature EOF\n");
+        exit(1);
+    }
+
+/* ==================================================================== */
+/*      Write out the NTv2 format grid shift file.                      */
+/* ==================================================================== */
+    if (!(fp = fopen(argv[1], "wb"))) 
+    {
+        perror(argv[1]);
+        exit(2);
+    }
+        
+/* -------------------------------------------------------------------- */
+/*      Write the file header.                                          */
+/* -------------------------------------------------------------------- */
+    {    
+        char achHeader[11*16];
+
+        memset( achHeader, 0, sizeof(achHeader) );
+        
+        memcpy( achHeader +  0*16, "NUM_OREC", 8 );
+        achHeader[ 0*16 + 8] = 0xb;
+
+        memcpy( achHeader +  1*16, "NUM_SREC", 8 );
+        achHeader[ 1*16 + 8] = 0xb;
+
+        memcpy( achHeader +  2*16, "NUM_FILE", 8 );
+        achHeader[ 2*16 + 8] = 0x1;
+
+        memcpy( achHeader +  3*16, "GS_TYPE         ", 16 );
+        memcpy( achHeader +  3*16+8, GS_TYPE, MIN(16,strlen(GS_TYPE)) );
+
+        memcpy( achHeader +  4*16, "VERSION         ", 16 );
+        memcpy( achHeader +  4*16+8, VERSION, MIN(16,strlen(VERSION)) );
+
+        memcpy( achHeader +  5*16, "SYSTEM_F        ", 16 );
+        memcpy( achHeader +  5*16+8, SYSTEM_F, MIN(16,strlen(SYSTEM_F)) );
+
+        memcpy( achHeader +  6*16, "SYSTEM_T        ", 16 );
+        memcpy( achHeader +  6*16+8, SYSTEM_T, MIN(16,strlen(SYSTEM_T)) );
+
+        memcpy( achHeader +  7*16, "MAJOR_F ", 8);
+        memcpy( achHeader +  8*16, "MINOR_F ", 8 );
+        memcpy( achHeader +  9*16, "MAJOR_T ", 8 );
+        memcpy( achHeader + 10*16, "MINOR_T ", 8 );
+
+        fwrite( achHeader, 1, sizeof(achHeader), fp );
+    }
+        
+/* -------------------------------------------------------------------- */
+/*      Write the grid header.                                          */
+/* -------------------------------------------------------------------- */
+    {
+        unsigned char achHeader[11*16];
+        double dfValue;
+        int nGSCount = ct.lim.lam * ct.lim.phi;
+        LP ur;
+
+        ur.lam = ct.ll.lam + (ct.lim.lam-1) * ct.del.lam;
+        ur.phi = ct.ll.phi + (ct.lim.phi-1) * ct.del.phi;
+
+        assert( sizeof(nGSCount) == 4 );
+
+        memset( achHeader, 0, sizeof(achHeader) );
+
+        memcpy( achHeader +  0*16, "SUB_NAME        ", 16 );
+        memcpy( achHeader +  0*16+8, SUB_NAME, MIN(16,strlen(SUB_NAME)) );
+    
+        memcpy( achHeader +  1*16, "PARENT          ", 16 );
+        memcpy( achHeader +  1*16+8, "NONE", MIN(16,strlen("NONE")) );
+    
+        memcpy( achHeader +  2*16, "CREATED         ", 16 );
+        memcpy( achHeader +  2*16+8, CREATED, MIN(16,strlen(CREATED)) );
+    
+        memcpy( achHeader +  3*16, "UPDATED         ", 16 );
+        memcpy( achHeader +  3*16+8, UPDATED, MIN(16,strlen(UPDATED)) );
+
+        memcpy( achHeader +  4*16, "S_LAT   ", 8 );
+        dfValue = ct.ll.phi * 3600.0 / DEG_TO_RAD;
+        memcpy( achHeader +  4*16 + 8, &dfValue, 8 );
+
+        memcpy( achHeader +  5*16, "N_LAT   ", 8 );
+        dfValue = ur.phi * 3600.0 / DEG_TO_RAD;
+        memcpy( achHeader +  5*16 + 8, &dfValue, 8 );
+
+        memcpy( achHeader +  6*16, "E_LONG  ", 8 );
+        dfValue = -1 * ur.lam * 3600.0 / DEG_TO_RAD;
+        memcpy( achHeader +  6*16 + 8, &dfValue, 8 );
+
+        memcpy( achHeader +  7*16, "W_LONG  ", 8 );
+        dfValue = -1 * ct.ll.lam * 3600.0 / DEG_TO_RAD;
+        memcpy( achHeader +  7*16 + 8, &dfValue, 8 );
+
+        memcpy( achHeader +  8*16, "LAT_INC ", 8 );
+        dfValue = ct.del.phi * 3600.0 / DEG_TO_RAD;
+        memcpy( achHeader +  8*16 + 8, &dfValue, 8 );
+    
+        memcpy( achHeader +  9*16, "LONG_INC", 8 );
+        dfValue = ct.del.lam * 3600.0 / DEG_TO_RAD;
+        memcpy( achHeader +  9*16 + 8, &dfValue, 8 );
+    
+        memcpy( achHeader + 10*16, "GS_COUNT", 8 );
+        memcpy( achHeader + 10*16+8, &nGSCount, 4 );
+    
+        if( !IS_LSB ) 
+        {
+            swap_words( achHeader +  4*16 + 8, 8, 1 );
+            swap_words( achHeader +  5*16 + 8, 8, 1 );
+            swap_words( achHeader +  6*16 + 8, 8, 1 );
+            swap_words( achHeader +  7*16 + 8, 8, 1 );
+            swap_words( achHeader +  8*16 + 8, 8, 1 );
+            swap_words( achHeader +  9*16 + 8, 8, 1 );
+            swap_words( achHeader + 10*16 + 8, 4, 1 );
+        }
+
+        fwrite( achHeader, 1, sizeof(achHeader), fp );
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Write the actual grid cells.                                    */
+/* -------------------------------------------------------------------- */
+    {
+        float *row_buf;
+        int row;
+
+        row_buf = (float *) pj_malloc(ct.lim.lam * sizeof(float) * 4);
+        memset( row_buf, 0, sizeof(float)*4 );
+
+        for( row = 0; row < ct.lim.phi; row++ )
+        {
+            int	    i;
+
+            for( i = 0; i < ct.lim.lam; i++ )
+            {
+                FLP *cvs = ct.cvs + (row) * ct.lim.lam
+                    + (ct.lim.lam - i - 1);
+
+                /* convert radians to seconds */
+                row_buf[i*4+0] = cvs->phi * (3600.0 / (PI/180.0));
+                row_buf[i*4+1] = cvs->lam * (3600.0 / (PI/180.0));
+
+                /* We leave the accuracy values as zero */
+            }
+
+            if( !IS_LSB )
+                swap_words( (unsigned char *) row_buf, 4, ct.lim.lam * 4 );
+
+            if( fwrite( row_buf, sizeof(float), ct.lim.lam*4, fp ) 
+                != 4 * ct.lim.lam )
+            {
+                perror( "write()" );
+                exit( 2 );
+            }
+        }
+    }
+
+    fclose( fp );
+    exit(0); /* normal completion */
 }