src/findentrances.cc

   1 /* findentrances.cc
   2  * Simple converter from survex .3d files to a list of entrances in GPX format
   3  *
   4  * Copyright (C) 2012 Olaf Kähler
   5  * Copyright (C) 2012,2013,2015 Olly Betts
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 2 of the License, or
  10  * (at your option) any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20  */
  21
  22 /*
  23  * This program parses the survex file, creates a list of entrances and does
  24  * the coordinate transformations from almost arbitrary formats into WGS84
  25  * using the PROJ.4 library. The output is then written as a GPX file ready
  26  * for use in your favourite GPS software. Everything else is kept as simple
  27  * and minimalistic as possible.
  28  *
  29  * Usage:
  30  *   findentrances [-d <+proj +datum +string>] <input.3d>
  31  *
  32  * Example for data given in BMN M31 (Totes Gebirge, Austria):
  33  *   findentrances -d '+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' cucc_austria.3d > ent.gpx
  34  *
  35  * Example for data given in british grid SD (Yorkshire):
  36  *   findentrances -d '+proj=tmerc +lat_0=49d +lon_0=-2d +k=0.999601 +x_0=100000 +y_0=-500000 +ellps=airy +towgs84=375,-111,431,0,0,0,0' yorkshire/all.3d > ent.gpx
  37  *
  38  * Example for data given as proper british grid reference:
  39  *   findentrances -d '+proj=tmerc +lat_0=49d +lon_0=-2d +k=0.999601 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=375,-111,431,0,0,0,0' all.3d > ent.gpx
  40  *
  41  */
  42
  43 #include <stdio.h>
  44 #include <string>
  45 #include <vector>
  46 #include <algorithm>
  47 #include <math.h>
  48
  49 #include <proj_api.h>
  50
  51 #include "message.h"
  52 #include "cmdline.h"
  53 #include "img_hosted.h"
  54
  55 using namespace std;
  56
  57 #define WGS84_DATUM_STRING "+proj=longlat +ellps=WGS84 +datum=WGS84"
  58
  59 struct Point {
  60     string label;
  61     double x, y, z;
  62 };
  63
  64 static void
  65 read_survey(const char *filename, vector<Point> & points, const char ** p_proj)
  66 {
  67     img *survey = img_open_survey(filename, NULL);
  68     if (!survey) {
  69         fatalerror(img_error2msg(img_error()), filename);
  70     }
  71
  72     if (survey->cs) {
  73         if (*p_proj) {
  74             fprintf(stderr, "3d file specifies the coordinate system, ignoring datum specified on command line\n");
  75         }
  76         *p_proj = osstrdup(survey->cs);
  77     } else if (!*p_proj) {
  78         cmdline_syntax();
  79         fprintf(stderr, "3d file does not specify the coordinate system - you need to specify the datum on the command line\n");
  80         exit(1);
  81     }
  82
  83     int result;
  84     do {
  85         img_point pt;
  86         result = img_read_item(survey, &pt);
  87
  88         if (result == img_LABEL) {
  89             if (survey->flags & img_SFLAG_ENTRANCE) {
  90                 Point newPt;
  91                 newPt.label.assign(survey->label);
  92                 newPt.x = pt.x;
  93                 newPt.y = pt.y;
  94                 newPt.z = pt.z;
  95                 points.push_back(newPt);
  96             }
  97         } else if (result == img_BAD) {
  98             img_close(survey);
  99             fatalerror(img_error2msg(img_error()), filename);
 100         }
 101     } while (result != img_STOP);
 102
 103     img_close(survey);
 104 }
 105
 106 static void
 107 convert_coordinates(vector<Point> & points, const char *inputDatum, const char *outputDatum)
 108 {
 109     projPJ pj_input, pj_output;
 110     if (!(pj_input = pj_init_plus(inputDatum))) {
 111         fatalerror(/*Failed to initialise input coordinate system “%s”*/287, inputDatum);
 112     }
 113     if (!(pj_output = pj_init_plus(outputDatum))) {
 114         fatalerror(/*Failed to initialise output coordinate system “%s”*/288, outputDatum);
 115     }
 116
 117     for (size_t i=0; i<points.size(); ++i) {
 118         pj_transform(pj_input, pj_output, 1, 1, &(points[i].x), &(points[i].y), &(points[i].z));
 119     }
 120 }
 121
 122 struct SortPointsByLabel {
 123     bool operator()(const Point & a, const Point & b)
 124     { return a.label < b.label; }
 125 } SortPointsByLabel;
 126
 127 static void
 128 sort_points(vector<Point> & points)
 129 {
 130     sort(points.begin(), points.end(), SortPointsByLabel);
 131 }
 132
 133 static void
 134 write_gpx(const vector<Point> & points, FILE *file)
 135 {
 136     fprintf(file, "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n<gpx version=\"1.0\" creator=\"survex - findentrances\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" xmlns=\"http://www.topografix.com/GPX/1/0\" xsi:schemaLocation=\"http://www.topografix.com/GPX/1/0 http://www.topografix.com/GPX/1/0/gpx.xsd\">\n");
 137     for (size_t i=0; i<points.size(); ++i) {
 138         const Point & pt = points[i];
 139         // %.8f is at worst just over 1mm.
 140         fprintf(file, "<wpt lon=\"%.8f\" lat=\"%.8f\"><ele>%.2f</ele><name>%s</name></wpt>\n", pt.x*180.0/M_PI, pt.y*180.0/M_PI, pt.z, pt.label.c_str());
 141     }
 142
 143     fprintf(file, "</gpx>\n");
 144 }
 145
 146 static const struct option long_opts[] = {
 147     {"datum", required_argument, 0, 'd'},
 148     {"help", no_argument, 0, HLP_HELP},
 149     {"version", no_argument, 0, HLP_VERSION},
 150     {0, 0, 0, 0}
 151 };
 152
 153 static const char *short_opts = "d:";
 154
 155 static struct help_msg help[] = {
 156 /*                              <-- */
 157    /* TRANSLATORS: The PROJ library is used to do coordinate transformations
 158     * (https://trac.osgeo.org/proj/) - the user passes a string to tell PROJ
 159     * what the input datum is. */
 160    {HLP_ENCODELONG(0),        /*input datum as string to pass to PROJ*/389, 0},
 161    {0, 0, 0}
 162 };
 163
 164 int main(int argc, char **argv)
 165 {
 166     msg_init(argv);
 167
 168     const char *datum_string = NULL;
 169     cmdline_set_syntax_message(/*-d PROJ_DATUM 3D_FILE*/388, 0, NULL);
 170     cmdline_init(argc, argv, short_opts, long_opts, NULL, help, 1, 1);
 171     while (1) {
 172         int opt = cmdline_getopt();
 173         if (opt == EOF) break;
 174         else if (opt == 'd') datum_string = optarg;
 175     }
 176
 177     const char *survey_filename = argv[optind];
 178
 179     vector<Point> points;
 180     read_survey(survey_filename, points, &datum_string);
 181     convert_coordinates(points, datum_string, WGS84_DATUM_STRING);
 182     sort_points(points);
 183     write_gpx(points, stdout);
 184 }