cavern: Support base location in Compass MAK

This is how Compass specifies a location to calculate magnetic
declination at.

We now handle this like Survex's native "*declination auto X Y Z".

doc/manual.sgml

@@ -4250,7 +4250,7 @@ can use a *equate with *case preserve on:
 </Para>
 
 <Para>
-    Survex understands basic MAK file features.  Known current limitations and
+    Survex understands most MAK file features.  Known current limitations and
     assumptions:
 
 <ItemizedList>
@@ -4308,10 +4308,15 @@ the MAK version you can change the A16 equate to:
 *equate file1.A16 file2.A16 file3.A16</programlisting>
     </Para></ListItem>
 
+    <ListItem><Para>
+	The @ command which specifes a base location to calculate magnetic
+	declinations at is handled, provided the specified datum is WGS 1984.
+	The UTM convergence angle specified as part of this command is ignored
+	as Survex knows how to calculate it.
+    </Para></ListItem>
     <ListItem><Para>
 	The following commands (and any other unknown commands) are currently
 	ignored:
-	    @ (Data about location to calculate declination at),
 	    % (Convergence angle (file-level)),
 	    * (Convergence angle (non file-level)),
 	    ! (Project parameters)

src/commands.c

@@ -747,6 +747,62 @@ report_declination(settings *p)
     }
 }
 
+void
+set_declination_location(real x, real y, real z, const char *proj_str)
+{
+    /* Convert to WGS84 lat long. */
+    PJ *transform = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
+					   proj_str,
+					   WGS84_DATUM_STRING,
+					   NULL);
+    if (transform) {
+	/* Normalise the output order so x is longitude and y latitude - by
+	 * default new PROJ has them switched for EPSG:4326 which just seems
+	 * confusing.
+	 */
+	PJ* pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX,
+						       transform);
+	proj_destroy(transform);
+	transform = pj_norm;
+    }
+
+    if (proj_angular_input(transform, PJ_FWD)) {
+	/* Input coordinate system expects radians. */
+	x = rad(x);
+	y = rad(y);
+    }
+
+    PJ_COORD coord = {{x, y, z, HUGE_VAL}};
+    coord = proj_trans(transform, PJ_FWD, coord);
+    x = coord.xyzt.x;
+    y = coord.xyzt.y;
+    z = coord.xyzt.z;
+
+    if (x == HUGE_VAL || y == HUGE_VAL || z == HUGE_VAL) {
+       compile_diagnostic(DIAG_ERR, /*Failed to convert coordinates: %s*/436,
+			  proj_errno_string(proj_errno(transform)));
+       /* Set dummy values which are finite. */
+       x = y = z = 0;
+    }
+    proj_destroy(transform);
+
+    report_declination(pcs);
+
+    double lon = rad(x);
+    double lat = rad(y);
+    pcs->z[Q_DECLINATION] = HUGE_REAL;
+    pcs->dec_lat = lat;
+    pcs->dec_lon = lon;
+    pcs->dec_alt = z;
+    pcs->dec_filename = file.filename;
+    pcs->dec_line = file.line;
+    pcs->dec_context = grab_line();
+    /* Invalidate cached declination. */
+    pcs->declination = HUGE_REAL;
+    /* Invalidate cached grid convergence. */
+    pcs->convergence = HUGE_REAL;
+}
+
 void
 pop_settings(void)
 {
@@ -1775,57 +1831,7 @@ cmd_declination(void)
 	    compile_diagnostic(DIAG_ERR, /*Input coordinate system must be specified for “*DECLINATION AUTO”*/301);
 	    return;
 	}
-	/* Convert to WGS84 lat long. */
-	PJ *transform = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
-					       pcs->proj_str,
-					       WGS84_DATUM_STRING,
-					       NULL);
-	if (transform) {
-	    /* Normalise the output order so x is longitude and y latitude - by
-	     * default new PROJ has them switched for EPSG:4326 which just seems
-	     * confusing.
-	     */
-	    PJ* pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX,
-							   transform);
-	    proj_destroy(transform);
-	    transform = pj_norm;
-	}
-
-	if (proj_angular_input(transform, PJ_FWD)) {
-	    /* Input coordinate system expects radians. */
-	    x = rad(x);
-	    y = rad(y);
-	}
-
-	PJ_COORD coord = {{x, y, z, HUGE_VAL}};
-	coord = proj_trans(transform, PJ_FWD, coord);
-	x = coord.xyzt.x;
-	y = coord.xyzt.y;
-	z = coord.xyzt.z;
-
-	if (x == HUGE_VAL || y == HUGE_VAL || z == HUGE_VAL) {
-	   compile_diagnostic(DIAG_ERR, /*Failed to convert coordinates: %s*/436,
-			      proj_errno_string(proj_errno(transform)));
-	   /* Set dummy values which are finite. */
-	   x = y = z = 0;
-	}
-	proj_destroy(transform);
-
-	report_declination(pcs);
-
-	double lon = rad(x);
-	double lat = rad(y);
-	pcs->z[Q_DECLINATION] = HUGE_REAL;
-	pcs->dec_lat = lat;
-	pcs->dec_lon = lon;
-	pcs->dec_alt = z;
-	pcs->dec_filename = file.filename;
-	pcs->dec_line = file.line;
-	pcs->dec_context = grab_line();
-	/* Invalidate cached declination. */
-	pcs->declination = HUGE_REAL;
-	/* Invalidate cached grid convergence. */
-	pcs->convergence = HUGE_REAL;
+	set_declination_location(x, y, z, pcs->proj_str);
     } else {
 	/* *declination D UNITS */
 	int units = get_units(BIT(Q_DECLINATION), false);

src/commands.h

@@ -35,6 +35,7 @@ void default_calib(settings *s);
 void pop_settings(void);
 void invalidate_pj_cached(void);
 void report_declination(settings *p);
+void set_declination_location(real x, real y, real z, const char *proj_str);
 
 void copy_on_write_meta(settings *s);
 

src/datain.c

@@ -533,6 +533,20 @@ nextch_handling_eol(void)
    }
 }
 
+static char *
+compass_utm_proj_str(int datum, int utm_zone)
+{
+    char *proj_str = osmalloc(32);
+    // We only support WGS84 currently.
+    (void)datum;
+    if (utm_zone > 0) {
+	sprintf(proj_str, "EPSG:%d", 32600 + utm_zone);
+    } else {
+	sprintf(proj_str, "EPSG:%d", 32700 - utm_zone);
+    }
+    return proj_str;
+}
+
 #define LITLEN(S) (sizeof(S"") - 1)
 #define has_ext(F,L,E) ((L) > LITLEN(E) + 1 &&\
 			(F)[(L) - LITLEN(E) - 1] == FNM_SEP_EXT &&\
@@ -605,6 +619,7 @@ data_file(const char *pth, const char *fnm)
       pcs = pcsNew;
       default_units(pcs);
       default_calib(pcs);
+      pcs->z[Q_DECLINATION] = HUGE_REAL;
 
       pcs->recorded_style = pcs->style = STYLE_NORMAL;
       pcs->units[Q_LENGTH] = METRES_PER_FOOT;
@@ -732,8 +747,12 @@ data_file(const char *pth, const char *fnm)
 	 get_token();
 	 nextch(); /* : */
 	 skipblanks();
-	 pcs->z[Q_DECLINATION] = -read_numeric(false);
-	 pcs->z[Q_DECLINATION] *= pcs->units[Q_DECLINATION];
+	 if (pcs->dec_filename == NULL) {
+	     pcs->z[Q_DECLINATION] = -read_numeric(false);
+	     pcs->z[Q_DECLINATION] *= pcs->units[Q_DECLINATION];
+	 } else {
+	     (void)read_numeric(false);
+	 }
 	 get_token();
 	 pcs->ordering = compass_order;
 	 if (strcmp(buffer, "FORMAT") == 0) {
@@ -796,6 +815,10 @@ data_file(const char *pth, const char *fnm)
    } else if (fmt == FMT_MAK) {
       int datum = 0;
       int utm_zone = 0;
+      real base_x = 0.0, base_y = 0.0, base_z = 0.0;
+      int base_utm_zone = 0;
+      unsigned int base_line = 0;
+      long base_lpos = 0;
       char *path = path_from_fnm(file.filename);
       int path_len = strlen(path);
       struct mak_folder {
@@ -817,6 +840,20 @@ data_file(const char *pth, const char *fnm)
 	       nextch_handling_eol();
 	    }
 	    if (dat_fnm) {
+	       if (base_utm_zone) {
+		   // Process the previous @ command using the datum from &.
+		   char *proj_str = compass_utm_proj_str(datum, base_utm_zone);
+		   // Temporarily reset line and lpos so dec_context and
+		   // dec_line refer to the @ command.
+		   unsigned saved_line = file.line;
+		   file.line = base_line;
+		   long saved_lpos = file.lpos;
+		   file.lpos = base_lpos;
+		   set_declination_location(base_x, base_y, base_z, proj_str);
+		   file.line = saved_line;
+		   file.lpos = saved_lpos;
+		   osfree(proj_str);
+	       }
 	       ch_store = ch;
 	       data_file(path, dat_fnm);
 	       ch = ch_store;
@@ -906,14 +943,10 @@ data_file(const char *pth, const char *fnm)
 	    pcs->proj_str = NULL;
 	    if (datum && utm_zone && abs(utm_zone) <= 60) {
 		/* Set up coordinate system. */
-		pcs->proj_str = osmalloc(32);
-		if (utm_zone > 0) {
-		    sprintf(pcs->proj_str, "EPSG:%d", 32600 + utm_zone);
-		} else {
-		    sprintf(pcs->proj_str, "EPSG:%d", 32700 - utm_zone);
-		}
+		char *proj_str = compass_utm_proj_str(datum, utm_zone);
+		pcs->proj_str = proj_str;
 		if (!proj_str_out) {
-		    proj_str_out = osstrdup(pcs->proj_str);
+		    proj_str_out = osstrdup(proj_str);
 		}
 	    }
 	    invalidate_pj_cached();
@@ -984,6 +1017,45 @@ data_file(const char *pth, const char *fnm)
 	    if (ch == ';') nextch_handling_eol();
 	    break;
 	  }
+	  case '@': {
+	    /* "Base Location" to calculate magnetic declination at:
+	     * UTM East, UTM North, Elevation, UTM Zone, Convergence Angle
+	     * The first three are in metres.
+	     */
+	    nextch();
+	    real easting = read_numeric(false);
+	    skipblanks();
+	    if (ch != ',') break;
+	    nextch();
+	    real northing = read_numeric(false);
+	    skipblanks();
+	    if (ch != ',') break;
+	    nextch();
+	    real elevation = read_numeric(false);
+	    skipblanks();
+	    if (ch != ',') break;
+	    nextch();
+	    int zone = read_int(-60, 60);
+	    skipblanks();
+	    if (ch != ',') break;
+	    nextch();
+	    real convergence_angle = read_numeric(false);
+	    /* We've now read them all successfully so store them.  The Compass
+	     * documentation gives an example which specifies the datum *AFTER*
+	     * the base location, so we need to convert lazily.
+	     */
+	    base_x = easting;
+	    base_y = northing;
+	    base_z = elevation;
+	    base_utm_zone = zone;
+	    base_line = file.line;
+	    base_lpos = file.lpos;
+	    // We ignore the stored UTM grid convergence angle since we get
+	    // this from PROJ.
+	    (void)convergence_angle;
+	    if (ch == ';') nextch_handling_eol();
+	    break;
+	  }
 	  default:
 	    nextch_handling_eol();
 	    break;

tests/Makefile.am

@@ -98,7 +98,7 @@ flags.dat flags.dump lrud.dat lrud.out lrud.pos\
 nomeasure.dat nomeasure.out nomeasure.pos noteam.dat noteam.out noteam.pos\
 fixfeet.mak fixfeet.pos\
 folder.mak subdir/cave1a.dat subdir/cave1b.dat subdir/subsubdir/cave2.dat\
-utm.mak utm.dump\
+utm.mak utm.out utm.dump\
 dot17.svx dot17.pos 3dcorner.svx 3dcorner.pos surfequate.svx surfequate.dxf\
 nosurveyhanging.svx nosurveyhanging.out\
 cmd_solve_hanging.svx passage.svx hanging_lrud.svx\

tests/backread.dat

@@ -11,8 +11,20 @@ DECLINATION: 0 FORMAT: DDDDLRUDLADB CORRECTIONS: 0.00 0.00 0.00
  C1 C3 3.281 999 0 1 1 1 1 180 999   Karst uses 999 for "no value" which Compass accepts for compatibility
 
 Cave
+SURVEY NAME: C
+SURVEY DATE: 1 1 1  COMMENT:Test backsights with 13 character format
+SURVEY TEAM: 
+Them
+DECLINATION: 0 FORMAT: DDDDLRUDLADBF CORRECTIONS: 0.00 0.00 0.00
+
+ FROM TO LENGTH BEARING INC LEFT UP DOWN RIGHT BACKBEARING BACKINC FLAGS COMMENTS
+
+ C1 B1 3.281 90 1 1 1 1 1 270 -1
+ B1 B2 3.281 180 -1 1 1 1 1 0 1
+
+Cave
 SURVEY NAME: D
-SURVEY DATE: 1 1 1  COMMENT:Test newer format backsights
+SURVEY DATE: 12 01 1999  COMMENT:Test newer format backsights
 SURVEY TEAM: 
 Us, Them
 DECLINATION: 0 FORMAT: DDDWLRUDLADadBF CORRECTIONS: 0.00 0.00 0.00 CORRECTIONS2: 1.0 -1.0

tests/backread.dump

@@ -5,10 +5,14 @@ VERSION 8
 SEPARATOR '.'
 --
 LEG 0.00 0.00 0.00 0.00 1.00 0.00 [] STYLE=NORMAL 1986.10.13
+LEG 0.00 0.00 0.00 1.00 0.00 0.02 [] STYLE=NORMAL
+LEG 1.00 0.00 0.02 1.00 -1.00 0.00 [] STYLE=NORMAL
 LEG 0.00 0.00 0.00 0.00 1.00 0.00 [] STYLE=NORMAL 1986.10.13
-LEG 0.00 1.00 0.00 0.59 2.15 0.52 [] STYLE=DIVING
+LEG 0.00 1.00 0.00 0.59 2.15 0.52 [] STYLE=DIVING 1999.12.01
 NODE 0.59 2.15 0.52 [D1] UNDERGROUND
 NODE 0.00 1.00 0.00 [C2] UNDERGROUND
+NODE 1.00 -1.00 0.00 [B2] UNDERGROUND
+NODE 1.00 0.00 0.02 [B1] UNDERGROUND
 NODE 0.00 1.00 0.00 [C3] UNDERGROUND
 NODE 0.00 0.00 0.00 [C1] UNDERGROUND
 STOP

tests/badinc5.out

@@ -17,16 +17,17 @@ Calculating trailing traverses...
 
 Calculating statistics...
 
-Survey contains 5 survey stations, joined by 3 legs.
+Survey contains 7 survey stations, joined by 5 legs.
 There are 0 loops.
 Survey has 2 connected components.
-Total length of survey legs =    3.39m (   3.39m adjusted)
-Total plan length of survey legs =    3.29m
-Total vertical length of survey legs =    0.52m
+Total length of survey legs =    5.39m (   5.39m adjusted)
+Total plan length of survey legs =    5.29m
+Total vertical length of survey legs =    0.56m
 Vertical range = 0.52m (from D1 at 1235.08m to C2 at 1234.56m)
-North-South range = 2.15m (from D1 at 22.62m to C1 at 20.47m)
-East-West range = 0.59m (from D1 at 10.82m to C2 at 10.23m)
+North-South range = 3.15m (from D1 at 22.62m to B2 at 19.47m)
+East-West range = 1.00m (from B2 at 11.23m to C2 at 10.23m)
    1 0-node.
-   2 1-nodes.
+   3 1-nodes.
    2 2-nodes.
+   1 3-node.
 There were 1 warning(s) and 1 error(s) - no output files produced.

tests/cavern.tst

@@ -144,10 +144,14 @@ for file in $TESTS ; do
   pos=
 
   case $file in
-    backread.dat|flags.dat|folder.mak|utm.mak)
+    backread.dat|flags.dat|folder.mak)
       pos=dump
       warn=0
       ;;
+    utm.mak)
+      pos=dump
+      warn=1
+      ;;
     *.dat)
       # .dat files can't start with a comment.  All the other .dat tests
       # have the same settings.

tests/fixfeet.pos

@@ -1,4 +1,6 @@
 ( Easting, Northing, Altitude )
+(   11.23,    20.47,  1234.58 ) B1
+(   11.23,    19.47,  1234.56 ) B2
 (   10.23,    20.47,  1234.56 ) C1
 (   10.23,    21.47,  1234.56 ) C2
 (   10.23,    21.47,  1234.56 ) C3