From 70cef2cd4fe5af3700f64413798c5ac826d888c2 Mon Sep 17 00:00:00 2001 From: Philip Schuchardt Date: Tue, 11 Aug 2026 19:55:12 -0500 Subject: [PATCH] cavern: Add *cs file to read the coordinate system from a file `*cs file FILENAME` reads the coordinate system description from a separate file, so it can be a PROJ string, WKT or PROJJSON. WKT and PROJJSON are made up of double quoted strings written over several lines, which is awkward to write in a .svx file, and they can express datums a PROJ string can't, such as NAD83(2011), ETRS89 and GDA2020. It also means a .prj file such as those which accompany ESRI shapefiles can be used directly, so a survey can be georeferenced to match GIS data you already have. The filename is resolved relative to the directory of the file the *cs is in, a .prj extension is tried if the name as given isn't found, and the lines of the file are joined with a space since the coordinate system is stored in the .3d file as a single line. This addresses the review of the patch on https://trac.survex.com/ticket/150: `file` is a new coordinate system keyword rather than an `@` prefix on `custom`, so `custom` still always means a literal PROJ/WKT string, and there's no UTF-8 BOM skipping - PROJ rejects a BOM and so does GDAL's .prj reader, so cavern now does too. --- doc/3dformat.htm | 9 ++++- doc/datafile.rst | 29 +++++++++++++++ src/commands.c | 88 +++++++++++++++++++++++++++++++++++++++++++++- src/filelist.h | 1 + tests/Makefile.am | 3 ++ tests/cavern.tst | 3 +- tests/csprj.out | 24 +++++++++++++ tests/csprj.pos | 4 +++ tests/csprj.prj | 33 +++++++++++++++++ tests/csprj.svx | 11 ++++++ tests/csprj2.out | 24 +++++++++++++ tests/csprj2.pos | 4 +++ tests/csprj2.prj | 1 + tests/csprj2.svx | 11 ++++++ tests/csprjbad.out | 36 +++++++++++++++++++ tests/csprjbad.svx | 13 +++++++ tests/csprjbom.prj | 1 + 17 files changed, 292 insertions(+), 3 deletions(-) create mode 100644 tests/csprj.out create mode 100644 tests/csprj.pos create mode 100644 tests/csprj.prj create mode 100644 tests/csprj.svx create mode 100644 tests/csprj2.out create mode 100644 tests/csprj2.pos create mode 100644 tests/csprj2.prj create mode 100644 tests/csprj2.svx create mode 100644 tests/csprjbad.out create mode 100644 tests/csprjbad.svx create mode 100644 tests/csprjbom.prj diff --git a/doc/3dformat.htm b/doc/3dformat.htm index b6e8ab42..17546c91 100644 --- a/doc/3dformat.htm +++ b/doc/3dformat.htm @@ -60,7 +60,14 @@ There's no length limit on this string. which can be passed to PROJ. For a coordinate system with an EPSG code EPSG: 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 ESRI: followed by the code number. +with ESRI: 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 +datum= 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).
  • Survey hierarchy separator character. Survey station names form a hierarchy, and this character separates levels in the hierarchy. E.g. 161.entrance.6. Defaults to . if this item is not diff --git a/doc/datafile.rst b/doc/datafile.rst index 86848568..76421e9e 100644 --- a/doc/datafile.rst +++ b/doc/datafile.rst @@ -523,6 +523,11 @@ Example ; 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,6 +565,30 @@ Description * ``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. diff --git a/src/commands.c b/src/commands.c index 1f247a1d..2e5e4fc8 100644 --- a/src/commands.c +++ b/src/commands.c @@ -33,6 +33,7 @@ #include "datain.h" #include "date.h" #include "debug.h" +#include "filelist.h" #include "filename.h" #include "message.h" #include "netbits.h" @@ -2466,6 +2467,7 @@ typedef enum { CS_EPSG, CS_ESRI, CS_EUR79Z30, + CS_FILE, CS_IJTSK, CS_IJTSK03, CS_JTSK, @@ -2483,6 +2485,7 @@ static const sztok cs_tab[] = { {"EPSG", CS_EPSG}, /* EPSG: */ {"ESRI", CS_ESRI}, /* ESRI: */ {"EUR79Z30", CS_EUR79Z30}, + {"FILE", CS_FILE}, /* FILE */ {"IJTSK", CS_IJTSK}, {"IJTSK03", CS_IJTSK03}, {"JTSK", CS_JTSK}, @@ -2497,6 +2500,71 @@ static const sztok cs_tab[] = { {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,6 +2632,20 @@ cmd_cs(void) 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,7 +2717,8 @@ cmd_cs(void) 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,6 +2729,7 @@ cmd_cs(void) case CS_NONE: break; case CS_CUSTOM: + case CS_FILE: /* proj_str already set */ break; case CS_EPSG: @@ -2777,6 +2861,7 @@ cmd_cs(void) /* 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,6 +2875,7 @@ cmd_cs(void) proj_context_errno_string(PJ_DEFAULT_CTX, proj_context_errno(PJ_DEFAULT_CTX))); skipline(); + free(proj_str); return; } proj_destroy(pj); diff --git a/src/filelist.h b/src/filelist.h index 215cbbff..694d5d78 100644 --- a/src/filelist.h +++ b/src/filelist.h @@ -24,3 +24,4 @@ #define EXT_SVX_MSG "msg" #define EXT_INI "ini" #define EXT_LOG "log" +#define EXT_PRJ "prj" diff --git a/tests/Makefile.am b/tests/Makefile.am index 51bd3a6a..e41e8f04 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -178,6 +178,9 @@ csbad.altout csbad.out csbad.svx\ 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\ diff --git a/tests/cavern.tst b/tests/cavern.tst index 51e4c564..662b4027 100755 --- a/tests/cavern.tst +++ b/tests/cavern.tst @@ -121,7 +121,8 @@ TESTS_= 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\ diff --git a/tests/csprj.out b/tests/csprj.out new file mode 100644 index 00000000..c4cb8d94 --- /dev/null +++ b/tests/csprj.out @@ -0,0 +1,24 @@ + +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) diff --git a/tests/csprj.pos b/tests/csprj.pos new file mode 100644 index 00000000..b2af0f79 --- /dev/null +++ b/tests/csprj.pos @@ -0,0 +1,4 @@ +( Easting, Northing, Altitude ) +( 917.71, 1109.86, 197.39 ) a +( 935.37, 1092.19, 198.26 ) b +( 887.82, 1109.86, 200.00 ) entrance diff --git a/tests/csprj.prj b/tests/csprj.prj new file mode 100644 index 00000000..bb40c2e4 --- /dev/null +++ b/tests/csprj.prj @@ -0,0 +1,33 @@ +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]]] diff --git a/tests/csprj.svx b/tests/csprj.svx new file mode 100644 index 00000000..f099422d --- /dev/null +++ b/tests/csprj.svx @@ -0,0 +1,11 @@ +; 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 diff --git a/tests/csprj2.out b/tests/csprj2.out new file mode 100644 index 00000000..cac0d37f --- /dev/null +++ b/tests/csprj2.out @@ -0,0 +1,24 @@ + +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) diff --git a/tests/csprj2.pos b/tests/csprj2.pos new file mode 100644 index 00000000..7fb569af --- /dev/null +++ b/tests/csprj2.pos @@ -0,0 +1,4 @@ +( Easting, Northing, Altitude ) +( 10873.35, 41819.44, 197.39 ) a +( 10891.01, 41801.78, 198.26 ) b +( 10843.46, 41819.44, 200.00 ) entrance diff --git a/tests/csprj2.prj b/tests/csprj2.prj new file mode 100644 index 00000000..d89bdac5 --- /dev/null +++ b/tests/csprj2.prj @@ -0,0 +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]]] diff --git a/tests/csprj2.svx b/tests/csprj2.svx new file mode 100644 index 00000000..6785fb35 --- /dev/null +++ b/tests/csprj2.svx @@ -0,0 +1,11 @@ +; 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 diff --git a/tests/csprjbad.out b/tests/csprjbad.out new file mode 100644 index 00000000..12f4acae --- /dev/null +++ b/tests/csprjbad.out @@ -0,0 +1,36 @@ +./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. diff --git a/tests/csprjbad.svx b/tests/csprjbad.svx new file mode 100644 index 00000000..990df42e --- /dev/null +++ b/tests/csprjbad.svx @@ -0,0 +1,13 @@ +; 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 diff --git a/tests/csprjbom.prj b/tests/csprjbom.prj new file mode 100644 index 00000000..1a1c8d4a --- /dev/null +++ b/tests/csprjbom.prj @@ -0,0 +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]]] -- 2.50.1 (Apple Git-155)