diff --git a/doc/manual.sgml b/doc/manual.sgml index 0cb34cae..321bb875 100644 --- a/doc/manual.sgml +++ b/doc/manual.sgml @@ -4250,7 +4250,7 @@ can use a *equate with *case preserve on: - Survex understands basic MAK file features. Known current limitations and + Survex understands most MAK file features. Known current limitations and assumptions: @@ -4308,10 +4308,15 @@ the MAK version you can change the A16 equate to: *equate file1.A16 file2.A16 file3.A16 + + 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. + 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) diff --git a/src/commands.c b/src/commands.c index 98802b0b..fdafb81d 100644 --- a/src/commands.c +++ b/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); diff --git a/src/commands.h b/src/commands.h index d0e514a3..66e69d5a 100644 --- a/src/commands.h +++ b/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); diff --git a/src/datain.c b/src/datain.c index eae7dc6c..5f7c4705 100644 --- a/src/datain.c +++ b/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; diff --git a/tests/Makefile.am b/tests/Makefile.am index fb1b5a81..4c53f3b7 100644 --- a/tests/Makefile.am +++ b/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\ diff --git a/tests/backread.dat b/tests/backread.dat index b327c844..24270424 100644 --- a/tests/backread.dat +++ b/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 diff --git a/tests/backread.dump b/tests/backread.dump index 25e470ad..5d3a542b 100644 --- a/tests/backread.dump +++ b/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 diff --git a/tests/badinc5.out b/tests/badinc5.out index ed8691c0..dba6e111 100644 --- a/tests/badinc5.out +++ b/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. diff --git a/tests/cavern.tst b/tests/cavern.tst index 4d76a6d1..113d864d 100755 --- a/tests/cavern.tst +++ b/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. diff --git a/tests/fixfeet.pos b/tests/fixfeet.pos index 90d5c444..6dd280da 100644 --- a/tests/fixfeet.pos +++ b/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 diff --git a/tests/utm.dump b/tests/utm.dump index 00638785..c1466ae7 100644 --- a/tests/utm.dump +++ b/tests/utm.dump @@ -5,11 +5,15 @@ CS EPSG:32760 VERSION 8 SEPARATOR '.' -- -LEG 1313799.91 5427953.18 30.00 1313799.91 5427954.18 30.00 [] STYLE=NORMAL 1986.10.13 -LEG 1313799.91 5427953.18 30.00 1313799.91 5427954.18 30.00 [] STYLE=NORMAL 1986.10.13 -LEG 1313799.91 5427954.18 30.00 1313800.50 5427955.33 30.52 [] STYLE=DIVING -NODE 1313800.50 5427955.33 30.52 [D1] UNDERGROUND -NODE 1313799.91 5427954.18 30.00 [C2] UNDERGROUND -NODE 1313799.91 5427954.18 30.00 [C3] UNDERGROUND +LEG 1313799.91 5427953.18 30.00 1313800.39 5427954.06 30.00 [] STYLE=NORMAL 1986.10.13 +LEG 1313799.91 5427953.18 30.00 1313800.90 5427953.07 30.02 [] STYLE=NORMAL +LEG 1313800.90 5427953.07 30.02 1313800.79 5427952.08 30.00 [] STYLE=NORMAL +LEG 1313799.91 5427953.18 30.00 1313800.39 5427954.06 30.00 [] STYLE=NORMAL 1986.10.13 +LEG 1313800.39 5427954.06 30.00 1313801.45 5427954.78 30.52 [] STYLE=DIVING 1999.12.01 +NODE 1313801.45 5427954.78 30.52 [D1] UNDERGROUND +NODE 1313800.39 5427954.06 30.00 [C2] UNDERGROUND +NODE 1313800.79 5427952.08 30.00 [B2] UNDERGROUND +NODE 1313800.90 5427953.07 30.02 [B1] UNDERGROUND +NODE 1313800.39 5427954.06 30.00 [C3] UNDERGROUND NODE 1313799.91 5427953.18 30.00 [C1] UNDERGROUND FIXED STOP diff --git a/tests/utm.mak b/tests/utm.mak index ea2667e0..7be5959d 100644 --- a/tests/utm.mak +++ b/tests/utm.mak @@ -1,3 +1,4 @@ +@1314000.0,5428000.0,50.0,-60,0.0; $-60; &WGS 1984; *fix bh reference 174.7767 -41.2784 30 diff --git a/tests/utm.out b/tests/utm.out new file mode 100644 index 00000000..b24199eb --- /dev/null +++ b/tests/utm.out @@ -0,0 +1,32 @@ +In file included from ./utm.mak:5: +./backread.dat:22: warning: No survey date specified - using 0 for magnetic declination + C1 B1 3.281 90 1 1 1 1 1 270 -1 +./utm.mak:1: info: Declination: 22.1dg @ 1986-10-13 / 22.4dg @ 1999-12-01, grid convergence: -6.4dg + @1314000.0,5428000.0,50.0,-60,0.0; + +Removing trailing traverses... + +Concatenating traverses... + +Simplifying network... + +Calculating network... + +Calculating traverses... + +Calculating trailing traverses... + +Calculating statistics... + +Survey contains 6 survey stations, joined by 5 legs. +There are 0 loops. +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 30.52m to C2 at 30.00m) +North-South range = 2.71m (from D1 at 5427954.78m to B2 at 5427952.08m) +East-West range = 1.54m (from D1 at 1313801.45m to C1 at 1313799.91m) + 3 1-nodes. + 2 2-nodes. + 1 3-node. +There were 1 warning(s).