Ticket #150: 0001-cavern-Add-cs-file-to-read-the-coordinate-system-fro.patch

doc/3dformat.htm

@@ -60 +60 @@
 which can be passed to PROJ.  For a coordinate system with an EPSG code
 <code>EPSG:</code> followed by the code number can be used (we recommend
 using this if an EPSG code exists).  Similarly, an ESRI code can be specified
-with <code>ESRI:</code> followed by the code number.
+with <code>ESRI:</code> followed by the code number.  Any other description
+PROJ accepts may be used, such as a PROJ string, WKT or PROJJSON - this is how
+a coordinate system based on a datum PROJ has no <code>+datum=</code> keyword
+for gets recorded.  Since the metadata ends at a linefeed, the string must be
+written on a single line - WKT and PROJJSON are usually pretty printed over
+several lines, so a writer needs to join those lines (a single space in place
+of each line break works, as neither format allows a newline inside a quoted
+name).
   <li>Survey hierarchy separator character.  Survey station names form a
 hierarchy, and this character separates levels in the hierarchy.  E.g.
 <code>161.entrance.6</code>.  Defaults to <code>.</code> if this item is not

doc/datafile.rst

@@ -523 +523 @@
        ; Output in the coordinate system used in the Totes Gebirge in Austria
        *cs out custom "+proj=tmerc +lat_0=0 +lon_0=13d20 +k=1 +x_0=0 +y_0=-5200000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232"
 
+   ::
+
+       ; Output in a low distortion projection described by mammoth.prj
+       *cs out file "mammoth.prj"
+
 Description
    ``*cs`` allows the coordinate systems used for fixed points and for
    processed survey data to be specified.
@@ -560 +565 @@
    * ``EUR79Z30`` for UTM zone 30, EUR79 datum.  Supported since Survex
      1.2.15.
 
+   * ``FILE`` followed by the name of a file which describes the coordinate
+     system (like in the example above).  The description may be a PROJ string,
+     but the point of reading it from a file is that it can also be WKT or
+     PROJJSON, which are made up of double quoted strings written over several
+     lines and so are awkward to write in a ``.svx`` file.  It means a ``.prj``
+     file such as those which accompany ESRI shapefiles can be used directly,
+     so you can georeference a survey to match GIS data you already have.
+
+     WKT can also express datums which a PROJ string can't, such as
+     NAD83(2011), ETRS89 and GDA2020.
+
+     The file is resolved relative to the directory which the file containing
+     the ``*cs`` is in, and if the name as specified is not found cavern will
+     try adding a ``.prj`` extension.  If the filename contains spaces, it must
+     be enclosed in double quotes, and you can double the quote character to
+     include it in the filename.
+
+     The lines of the file are joined with a single space, since the coordinate
+     system is stored in the ``.3d`` file as a single line.  Neither WKT nor
+     PROJJSON allows a newline inside a quoted name, so only insignificant
+     whitespace is affected.
+
+     Supported since Survex 1.4.23.
+
    * ``IJTSK`` for the modified version of the Czechoslovak S-JTSK system
      where the axes point East and North.  Supported since Survex 1.2.15.
 

src/commands.c

@@ -33 +33 @@
 #include "datain.h"
 #include "date.h"
 #include "debug.h"
+#include "filelist.h"
 #include "filename.h"
 #include "message.h"
 #include "netbits.h"
@@ -2466 +2467 @@
     CS_EPSG,
     CS_ESRI,
     CS_EUR79Z30,
+    CS_FILE,
     CS_IJTSK,
     CS_IJTSK03,
     CS_JTSK,
@@ -2483 +2485 @@
      {"EPSG",     CS_EPSG},     /* EPSG:<number> */
      {"ESRI",     CS_ESRI},     /* ESRI:<number> */
      {"EUR79Z30", CS_EUR79Z30},
+     {"FILE",     CS_FILE},     /* FILE <filename> */
      {"IJTSK",    CS_IJTSK},
      {"IJTSK03",  CS_IJTSK03},
      {"JTSK",     CS_JTSK},
@@ -2497 +2500 @@
      {NULL,       CS_NONE}
 };
 
+/* Read a coordinate system description from the file FNM, as specified by
+ * `*cs file FILENAME`.  FP is the position of FILENAME, which any diagnostic
+ * is reported against.
+ *
+ * Returns the description, or NULL if the file couldn't be opened (in which
+ * case a diagnostic has been reported).
+ */
+static char *
+read_cs_from_file(const char *fnm, const filepos *fp)
+{
+   char *pth = path_from_fnm(file.filename);
+   char *fnm_used = NULL;
+   FILE *fh = fopen_portable(pth, fnm, EXT_PRJ, "rb", &fnm_used);
+   free(pth);
+   if (fh == NULL) {
+      set_pos(fp);
+      compile_diagnostic(DIAG_ERR|DIAG_STRING, /*Couldn’t open file “%s”*/1,
+                         fnm);
+      return NULL;
+   }
+
+   /* We store the coordinate system in the .3d file as a single line, so join
+    * the lines with a space, dropping blanks at the start and end of each.
+    * Neither WKT nor PROJJSON allows a newline inside a quoted name, so only
+    * insignificant whitespace is affected.
+    */
+   string cs = S_INIT;
+   /* s_clear() gives cs a buffer, which s_steal() needs even for an empty
+    * file. */
+   s_clear(&cs);
+   string blanks = S_INIT;
+   bool line_break = false;
+   for (int c = GETC(fh); c != EOF; c = GETC(fh)) {
+      if (c == '\n' || c == '\r') {
+         line_break = true;
+         s_clear(&blanks);
+         continue;
+      }
+      if (c == ' ' || c == '\t') {
+         /* Only keep blanks which turn out to be between two non-blanks on
+          * the same line. */
+         if (cs.len) s_appendch(&blanks, c);
+         continue;
+      }
+      if (cs.len) {
+         if (line_break) {
+            s_appendch(&cs, ' ');
+         } else {
+            s_appends(&cs, &blanks);
+         }
+      }
+      s_clear(&blanks);
+      line_break = false;
+      s_appendch(&cs, c);
+   }
+   s_free(&blanks);
+
+   if (FERROR(fh))
+      fatalerror_in_file(fnm_used, 0, /*Error reading file*/18);
+   fclose(fh);
+   free(fnm_used);
+
+   return s_steal(&cs);
+}
+
 static void
 cmd_cs(void)
 {
@@ -2564 +2632 @@
            proj_str = s_steal(&str);
            cs_sub = 0;
            break;
+         case CS_FILE: {
+           ok_for_output = MAYBE;
+           get_pos(&fp);
+           string fnm_str = S_INIT;
+           read_string(&fnm_str);
+           proj_str = read_cs_from_file(s_str(&fnm_str), &fp);
+           s_free(&fnm_str);
+           if (!proj_str) {
+              skipline();
+              return;
+           }
+           cs_sub = 0;
+           break;
+         }
          case CS_EPSG: case CS_ESRI:
            ok_for_output = MAYBE;
            if (ch == ':' && isdigit(nextch())) {
@@ -2635 +2717 @@
            break;
        }
    }
-   if (cs_sub == INT_MIN || (cs != CS_CUSTOM && isalnum(ch))) {
+   if (cs_sub == INT_MIN ||
+       (cs != CS_CUSTOM && cs != CS_FILE && isalnum(ch))) {
       set_pos(&fp);
       compile_diagnostic(DIAG_ERR|DIAG_WORD, /*Unknown coordinate system*/434);
       skipline();
@@ -2646 +2729 @@
       case CS_NONE:
          break;
       case CS_CUSTOM:
+      case CS_FILE:
          /* proj_str already set */
          break;
       case CS_EPSG:
@@ -2777 +2861 @@
          /* Same as the current output projection, so valid for input. */
       } else if (pcs->proj_str && strcmp(proj_str, pcs->proj_str) == 0) {
          /* Same as the current input projection, so nothing to do! */
+         free(proj_str);
          return;
       } else if (ok_for_output == MAYBE) {
          /* (ok_for_output == MAYBE) also happens to indicate whether we need
@@ -2790 +2875 @@
                                proj_context_errno_string(PJ_DEFAULT_CTX,
                                                          proj_context_errno(PJ_DEFAULT_CTX)));
             skipline();
+            free(proj_str);
             return;
          }
          proj_destroy(pj);

src/filelist.h

@@ -24 +24 @@
 #define EXT_SVX_MSG  "msg"
 #define EXT_INI      "ini"
 #define EXT_LOG      "log"
+#define EXT_PRJ      "prj"

tests/Makefile.am

@@ -178 +178 @@
 csbadsdfix.altout csbadsdfix.out csbadsdfix.svx\
 csfeet.out csfeet.pos csfeet.svx\
 cslonglat.out cslonglat.svx\
+csprj.out csprj.pos csprj.prj csprj.svx\
+csprj2.out csprj2.pos csprj2.prj csprj2.svx\
+csprjbad.out csprjbad.svx csprjbom.prj\
 omitfixaroundsolve.out omitfixaroundsolve.svx\
 repeatreading.svx repeatreading.out repeatreading.pos\
 mixedeols.out mixedeols.svx\

tests/cavern.tst

@@ -121 +121 @@
  passage hanging_lrud equatenosuchstn surveytypo\
  skipafterbadomit passagebad badreadingdotplus badcalibrate calibrate_clino\
  badunits badbegin anonstn anonstnbad anonstnrev doubleinc reenterlots\
- cs csbad csbadsdfix csfeet cslonglat omitfixaroundsolve repeatreading\
+ cs csbad csbadsdfix csfeet cslonglat csprj csprj2 csprjbad\
+ omitfixaroundsolve repeatreading\
  mixedeols utf8bom nonewlineateof suspectreadings cmd_data_default\
  cmd_data_ignore\
  quadrant_bearing bad_quadrant_bearing\

new file tests/csprj.out

@@ +1 @@
+
+Removing trailing traverses...
+
+Concatenating traverses...
+
+Simplifying network...
+
+Calculating network...
+
+Calculating traverses...
+
+Calculating trailing traverses...
+
+Calculating statistics...
+
+Approximate full range of grid convergence: 0.0dg at entrance to 0.0dg at b
+Survey contains 3 survey stations, joined by 2 legs.
+There are 0 loops.
+Total length of survey legs =   55.00m (  55.00m adjusted)
+Total plan length of survey legs =   54.87m
+Total vertical length of survey legs =    3.49m
+Vertical range = 2.61m (from entrance at 200.00m to a at 197.39m)
+North-South range = 17.67m (from a at 1109.86m to b at 1092.19m)
+East-West range = 47.55m (from b at 935.37m to entrance at 887.82m)

new file tests/csprj.pos

@@ +1 @@
+( Easting, Northing, Altitude )
+(  917.71,  1109.86,   197.39 ) a
+(  935.37,  1092.19,   198.26 ) b
+(  887.82,  1109.86,   200.00 ) entrance

new file tests/csprj.prj

@@ +1 @@
+PROJCRS["Mammoth Cave LDP",
+    BASEGEOGCRS["NAD83(2011)",
+        DATUM["NAD83 (National Spatial Reference System 2011)",
+            ELLIPSOID["GRS 1980",6378137,298.257222101,
+                LENGTHUNIT["metre",1]]],
+        PRIMEM["Greenwich",0,
+            ANGLEUNIT["degree",0.0174532925199433]],
+        ID["EPSG",6318]],
+    CONVERSION["Transverse Mercator",
+        METHOD["Transverse Mercator",
+            ID["EPSG",9807]],
+        PARAMETER["Latitude of natural origin",37.1866,
+            ANGLEUNIT["degree",0.0174532925199433],
+            ID["EPSG",8801]],
+        PARAMETER["Longitude of natural origin",-86.1005,
+            ANGLEUNIT["degree",0.0174532925199433],
+            ID["EPSG",8802]],
+        PARAMETER["Scale factor at natural origin",1,
+            SCALEUNIT["unity",1],
+            ID["EPSG",8805]],
+        PARAMETER["False easting",0,
+            LENGTHUNIT["metre",1],
+            ID["EPSG",8806]],
+        PARAMETER["False northing",0,
+            LENGTHUNIT["metre",1],
+            ID["EPSG",8807]]],
+    CS[Cartesian,2],
+    AXIS["(E)",east,
+        ORDER[1],
+        LENGTHUNIT["metre",1]],
+    AXIS["(N)",north,
+        ORDER[2],
+        LENGTHUNIT["metre",1]]]

new file tests/csprj.svx

@@ +1 @@
+; pos=yes warn=0
+; Test *cs file FILENAME, which reads the coordinate system from a separate
+; file.  Here that's a low distortion projection for a cave, written as WKT
+; over several lines - a PROJ string can't express the datum it's based on.
+*cs EPSG:6318
+*cs out FILE csprj.prj
+
+*fix entrance -86.0905 37.1966 200
+
+entrance a 30.0 090 -05
+a        b 25.0 135 +02

new file tests/csprj2.out

@@ +1 @@
+
+Removing trailing traverses...
+
+Concatenating traverses...
+
+Simplifying network...
+
+Calculating network...
+
+Calculating traverses...
+
+Calculating trailing traverses...
+
+Calculating statistics...
+
+Approximate full range of grid convergence: 0.0dg at entrance to 0.0dg at b
+Survey contains 3 survey stations, joined by 2 legs.
+There are 0 loops.
+Total length of survey legs =   55.00m (  55.00m adjusted)
+Total plan length of survey legs =   54.87m
+Total vertical length of survey legs =    3.49m
+Vertical range = 2.61m (from entrance at 200.00m to a at 197.39m)
+North-South range = 17.67m (from a at 41819.44m to b at 41801.78m)
+East-West range = 47.55m (from b at 10891.01m to entrance at 10843.46m)

new file tests/csprj2.pos

@@ +1 @@
+( Easting, Northing, Altitude )
+( 10873.35, 41819.44,   197.39 ) a
+( 10891.01, 41801.78,   198.26 ) b
+( 10843.46, 41819.44,   200.00 ) entrance

new file tests/csprj2.prj

@@ +1 @@
+PROJCRS["Mammoth Cave LDP", BASEGEOGCRS["NAD83(2011)", DATUM["NAD83 (National Spatial Reference System 2011)", ELLIPSOID["GRS 1980",6378137,298.257222101, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",6318]], CONVERSION["Mammoth Cave LDP TM", METHOD["Transverse Mercator", ID["EPSG",9807]], PARAMETER["Latitude of natural origin",37, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]], PARAMETER["Longitude of natural origin",-86.1, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]], PARAMETER["Scale factor at natural origin",1.0000383, SCALEUNIT["unity",1], ID["EPSG",8805]], PARAMETER["False easting",10000, LENGTHUNIT["metre",1], ID["EPSG",8806]], PARAMETER["False northing",20000, LENGTHUNIT["metre",1], ID["EPSG",8807]]], CS[Cartesian,2], AXIS["easting (E)",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["northing (N)",north, ORDER[2], LENGTHUNIT["metre",1]]]

new file tests/csprj2.svx

@@ +1 @@
+; pos=yes warn=0
+; Test that *cs file FILENAME assumes a .prj extension if the name as given
+; doesn't exist, that the filename can be quoted, and that a description
+; written on a single line (as ESRI .prj files are) works.
+*cs EPSG:6318
+*cs out FILE "csprj2"
+
+*fix entrance -86.0905 37.1966 200
+
+entrance a 30.0 090 -05
+a        b 25.0 135 +02

new file tests/csprjbad.out

@@ +1 @@
+./csprjbad.svx:3:14: error: Couldn't open file "nosuchfile.prj"
+ *cs out FILE nosuchfile.prj
+              ^~~~~~~~~~~~~~
+./csprjbad.svx:5:14: error: Invalid coordinate system: Invalid PROJ string syntax
+ *cs out FILE csprjbom.prj
+              ^~~~~~~~~~~~
+./csprjbad.svx:8:21: error: End of line not blank
+ *cs FILE "csprj.prj"x
+                     ^
+./csprjbad.svx:11: error: The input projection is set but the output projection isn't
+ *fix entrance -86.0905 37.1966 200
+
+Removing trailing traverses...
+
+Concatenating traverses...
+
+Simplifying network...
+
+Calculating network...
+
+Calculating traverses...
+
+Calculating trailing traverses...
+
+Calculating statistics...
+
+Survey contains 2 survey stations, joined by 1 leg.
+There are 0 loops.
+Total length of survey legs =   30.00m (  30.00m adjusted)
+Total plan length of survey legs =   29.89m
+Total vertical length of survey legs =    2.61m
+Vertical range = 2.61m (from entrance at 200.00m to a at 197.39m)
+North-South range = 0.00m (from a at 37.20m to a at 37.20m)
+East-West range = 29.89m (from a at -56.20m to entrance at -86.09m)
+
+There were 0 warning(s) and 4 error(s) - no output files produced.

new file tests/csprjbad.svx

@@ +1 @@
+; pos=fail warn=0 error=4
+; Test errors from *cs file FILENAME.
+*cs out FILE nosuchfile.prj
+; PROJ rejects a description which starts with a UTF-8 "BOM".
+*cs out FILE csprjbom.prj
+; Junk after the filename is reported as such, not as an unknown coordinate
+; system.
+*cs FILE "csprj.prj"x
+*cs EPSG:6318
+
+*fix entrance -86.0905 37.1966 200
+
+entrance a 30.0 090 -05

new file tests/csprjbom.prj

@@ +1 @@
+PROJCRS["Mammoth Cave LDP", BASEGEOGCRS["NAD83(2011)", DATUM["NAD83 (National Spatial Reference System 2011)", ELLIPSOID["GRS 1980",6378137,298.257222101, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",6318]], CONVERSION["Mammoth Cave LDP TM", METHOD["Transverse Mercator", ID["EPSG",9807]], PARAMETER["Latitude of natural origin",37, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]], PARAMETER["Longitude of natural origin",-86.1, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]], PARAMETER["Scale factor at natural origin",1.0000383, SCALEUNIT["unity",1], ID["EPSG",8805]], PARAMETER["False easting",10000, LENGTHUNIT["metre",1], ID["EPSG",8806]], PARAMETER["False northing",20000, LENGTHUNIT["metre",1], ID["EPSG",8807]]], CS[Cartesian,2], AXIS["easting (E)",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["northing (N)",north, ORDER[2], LENGTHUNIT["metre",1]]]