Ticket #150: 0001-cavern-Allow-cs-custom-to-read-the-CRS-from-a-file.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

@@ -501 +501 @@
        ; 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 custom @mammoth.prj
+
 Description
    ``*cs`` allows the coordinate systems used for fixed points and for
    processed survey data to be specified.
@@ -529 +534 @@
      resource for finding the EPSG code you want.  For example, ``EPSG:4167``
      is NZGD2000.  Supported since Survex 1.2.15.
 
-   * ``CUSTOM`` followed by a PROJ string (like in the example above).
+   * ``CUSTOM`` followed by a PROJ string, or by ``@`` and the name of a file
+     containing a description of the coordinate system - the last two examples
+     above show each of these in turn.
+
+     Such a file may contain a PROJ string, but the point of it is that it can
+     also hold 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.  This 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.
+
+     It's worth using WKT when the datum you need is one PROJ has no
+     ``+datum=`` keyword for, which is the case for most datums in current use
+     (such as NAD83(2011), ETRS89 and GDA2020).  A PROJ string can only express
+     a datum as an ellipsoid plus a static shift, so PROJ quietly drops such a
+     datum when it converts a coordinate system to a PROJ string.
+
+     The filename is relative to the directory containing the file the ``*cs``
+     command is in, and an extension of ``.prj`` is assumed if the name as
+     given doesn't exist.  Use double quotes around the name if it contains
+     spaces - the ``@`` may go either side of the opening quote, so both
+     ``@"My Cave.prj"`` and ``"@My Cave.prj"`` work.
+
+     The lines of the file are joined together with a single space, since the
+     coordinate system is stored in the ``.3d`` file as a single line.  This
+     only affects whitespace which the formats treat as insignificant, because
+     neither WKT nor PROJJSON allows a newline inside a quoted name.
+
+     Reading the coordinate system from a file was added in Survex 1.4.23.
 
    * ``ESRI:`` followed by a positive integer code.  ESRI codes are used by
      ArcGIS to specify coordinate systems (in a similar way to EPSG codes)

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"
@@ -2486 +2487 @@
      {NULL,       CS_NONE}
 };
 
+/* Read a coordinate system from the file FNM, as specified by
+ * `*cs custom @FILENAME`.  FP is the position of FILENAME in the current file,
+ * which any diagnostic is reported against.
+ *
+ * WKT and PROJJSON are made up of double quoted strings and are usually
+ * written over several lines, which makes them awkward to write in a .svx
+ * file.  Keeping such a description in its own file also means a .prj file
+ * such as those which accompany ESRI shapefiles can be used directly.
+ *
+ * Returns the coordinate system 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;
+   }
+
+   int c = GETC(fh);
+   if (c == 0xef) {
+      /* Skip a UTF-8 "BOM" if there is one - PROJ rejects a description which
+       * starts with one. */
+      if (GETC(fh) == 0xbb && GETC(fh) == 0xbf) {
+         c = GETC(fh);
+      } else {
+         rewind(fh);
+         c = GETC(fh);
+      }
+   }
+
+   /* We store the coordinate system in the .3d file as part of a
+    * newline-terminated line, so it can't contain a newline.  Join the lines
+    * with a single space, dropping blanks at the start and end of each line.
+    * Neither WKT nor PROJJSON allows a newline inside a quoted name, so only
+    * insignificant whitespace is affected.
+    */
+   /* The s_clear() calls give each string a buffer, which for an empty file it
+    * would otherwise still lack by the time we use it - s_steal() writes the
+    * terminating zero byte to one, and s_appends() reads from one.
+    */
+   string cs = S_INIT;
+   s_clear(&cs);
+   string blanks = S_INIT;
+   s_clear(&blanks);
+   bool line_break = false;
+   for ( ; 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)
 {
@@ -2545 +2632 @@
        switch (cs) {
          case CS_NONE:
            break;
-         case CS_CUSTOM:
+         case CS_CUSTOM: {
            ok_for_output = MAYBE;
+           skipblanks();
+           /* `@FILENAME` reads the coordinate system from a file.  If FILENAME
+            * is quoted then the `@` may be written either side of the opening
+            * quote.
+            */
+           bool from_file = (ch == '@');
+           if (from_file) nextch();
            get_pos(&fp);
            string str = S_INIT;
            read_string(&str);
-           proj_str = s_steal(&str);
+           const char *p = s_str(&str);
+           if (!from_file && *p == '@') {
+              from_file = true;
+              ++p;
+           }
+           if (from_file) {
+              proj_str = read_cs_from_file(p, &fp);
+              s_free(&str);
+              if (!proj_str) {
+                 skipline();
+                 return;
+              }
+           } else {
+              proj_str = s_steal(&str);
+           }
            cs_sub = 0;
            break;
+         }
          case CS_EPSG: case CS_ESRI:
            ok_for_output = MAYBE;
            if (ch == ':' && isdigit(nextch())) {
@@ -2766 +2875 @@
          /* 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
@@ -2779 +2889 @@
                                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

@@ -177 +177 @@
 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\
 omitfixaroundsolve.out omitfixaroundsolve.svx\
 repeatreading.svx repeatreading.out repeatreading.pos\
 mixedeols.out mixedeols.svx\

tests/cavern.tst

@@ -120 +120 @@
  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 custom @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 CUSTOM @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 custom @FILENAME assumes a .prj extension if the name as given
+; doesn't exist, and that a UTF-8 "BOM" at the start of the file is skipped.
+; The coordinate system is quoted here too, with the `@` inside the quotes.
+*cs EPSG:6318
+*cs out CUSTOM "@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:4:17: error: Couldn't open file "nosuchfile.prj"
+ *cs out CUSTOM @nosuchfile.prj
+                 ^~~~~~~~~~~~~~
+./csprjbad.svx:6: 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 2 error(s) - no output files produced.

new file tests/csprjbad.svx

@@ +1 @@
+; pos=fail error=2
+; Test the error from *cs custom @FILENAME naming a file which isn't there.
+*cs EPSG:6318
+*cs out CUSTOM @nosuchfile.prj
+
+*fix entrance -86.0905 37.1966 200
+
+entrance a 30.0 090 -05