home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / mapproject.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-18  |  9.8 KB  |  289 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)mapproject.c    2.17  18 Jun 1995
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  * mapproject reads a pair of coordinates [+ optional data fields] from
  9.  * standard input or file(s) and transforms the coordinates according to the
  10.  * map projection selected. See man gmt-system for projections currently supported.
  11.  *
  12.  * The default is to expect longitude, latitude, [and optional datavalues],
  13.  * and return x, y, [ and optional datavalues], but if the -I option is used,
  14.  * the reverse is true.  Specifying -C means that the origin of projected coordinates
  15.  * should be set to origin of projection  [Default origin is lower left corner of "map"].
  16.  * If your data is lat lon instead of lon lat [Default], use -: to toggle x/y -> y/x.
  17.  * Note that only unprojected values are affected by the -: switch.  True x,y values are
  18.  * always printed out as x,y.
  19.  *
  20.  *
  21.  * Author:    Paul Wessel
  22.  * Date:    1-MAR-1990
  23.  * Version:    2.0, based on old version 1.1
  24.  *
  25.  */
  26.  
  27. #include "gmt.h"
  28.  
  29. main (argc, argv)
  30. int argc;
  31. char **argv; {
  32.     int i, fno, n = 0, n_read = 0, n_cols = 0, n_files = 0, x, y, n_args;
  33.     
  34.     BOOLEAN error = FALSE, inverse = FALSE, suppress = FALSE, one_to_one = FALSE;
  35.     BOOLEAN map_center = FALSE, nofile = TRUE, done = FALSE, first = TRUE, km = FALSE;
  36.     BOOLEAN multi_segments = FALSE;
  37.     
  38.     double west = 0.0, east = 0.0, south = 0.0, north = 0.0;
  39.     double xmin, xmax, ymin, ymax, xy_in[2], xy_out[2];
  40.     double x_in_min, x_in_max, y_in_min, y_in_max;
  41.     double x_out_min, x_out_max, y_out_min, y_out_max, u_scale;
  42.     
  43.     char data[512], line[512], format1[512], format2[512], EOL_flag = '>';
  44.     
  45.     FILE *fp = NULL;
  46.     
  47.     argc = gmt_begin (argc, argv);
  48.     
  49.     for (i = 1; i < argc; i++) {
  50.         if (argv[i][0] == '-') {
  51.             switch (argv[i][1]) {
  52.         
  53.                 /* Common parameters */
  54.             
  55.                 case 'H':
  56.                 case 'J':
  57.                 case 'R':
  58.                 case 'V':
  59.                 case ':':
  60.                 case '\0':
  61.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  62.                     break;
  63.                 
  64.                 /* Supplemental parameters */
  65.                 
  66.                 case 'C':
  67.                     map_center = TRUE;
  68.                     break;
  69.                 case 'F':
  70.                     one_to_one = TRUE;
  71.                     if (argv[i][2] == 'k' || argv[i][2] == 'K') km = TRUE;
  72.                     if (argv[i][2] == 'm' || argv[i][2] == 'M') km = 2;
  73.                     break;
  74.                 case 'I':
  75.                     inverse = TRUE;
  76.                     break;
  77.                 case 'M':               /* Multiple line segments */
  78.                     multi_segments = TRUE;
  79.                     if (argv[i][2]) EOL_flag = argv[i][2];
  80.                     break;
  81.                 case 'S':
  82.                     suppress = TRUE;
  83.                     break;
  84.                 default:
  85.                     error = TRUE;
  86.                     gmt_default_error (argv[i][1]);
  87.                     break;
  88.             }
  89.         }
  90.         else
  91.             n_files++;
  92.     }
  93.     
  94.     if (argc == 1 || gmt_quick) {
  95.         fprintf (stderr, "mapproject %s - Forward and Inverse map transformation of 2-D coordinates\n\n", GMT_VERSION);
  96.         fprintf (stderr, "usage: mapproject <infiles> -J<parameters> -R<west/east/south/north>\n");
  97.         fprintf (stderr, "    [-C] [-F[k|m]] [-H] [-I] [-M[<flag>]] [-S] [-V] [-:]\n\n");
  98.         
  99.         if (gmt_quick) exit (-1);
  100.         
  101.         fprintf (stderr, "    infiles (in ASCII) has 2 or more columns.  If no file(s) is given, standard input is read.\n");
  102.         explain_option ('J');
  103.         explain_option ('R');
  104.         fprintf (stderr, "\n\tOPTIONS:\n");
  105.         fprintf (stderr, "    -C returns x/y relative to projection center [Default is relative to lower left corner]\n");
  106.         fprintf (stderr, "    -F force projected values to be in actual meters [Default uses the given plot scale]\n");
  107.         fprintf (stderr, "       Append k to use km as the unit or m for miles\n");
  108.         explain_option ('H');
  109.         fprintf (stderr, "    -I means Inverse, i.e., get lon/lat from x/y input. [Default is lon/lat -> x/y]\n");
  110.         fprintf (stderr, "    -M Input files each consist of multiple segments separated by a record\n");
  111.         fprintf (stderr, "       whose first character is <flag> ['>']\n");
  112.         fprintf (stderr, "    -S means Suppress points outside region\n");
  113.         explain_option ('V');
  114.         explain_option (':');
  115.         explain_option ('.');
  116.         exit (-1);
  117.     }
  118.     
  119.     if (!project_info.region_supplied) {
  120.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
  121.         error++;
  122.     }
  123.  
  124.     if (error) exit (-1);
  125.  
  126.     map_setup (west, east, south, north);
  127.     
  128.     if (gmtdefs.verbose) {
  129.         xmin = (map_center) ? project_info.xmin - project_info.x0 : project_info.xmin;
  130.         xmax = (map_center) ? project_info.xmax - project_info.x0 : project_info.xmax;
  131.         ymin = (map_center) ? project_info.ymin - project_info.y0 : project_info.ymin;
  132.         ymax = (map_center) ? project_info.ymax - project_info.y0 : project_info.ymax;
  133.         if (one_to_one) {
  134.             xmin /= project_info.x_scale;
  135.             xmax /= project_info.x_scale;
  136.             ymin /= project_info.y_scale;
  137.             ymax /= project_info.y_scale;
  138.         }
  139.         fprintf (stderr, "mapproject:  Transform (%lg/%lg/%lg/%lg ",
  140.             project_info.w, project_info.e, project_info.s, project_info.n);
  141.         (inverse) ? fprintf (stderr, " <- ") : fprintf (stderr, " -> ");
  142.         fprintf (stderr, "%lg/%lg/%lg/%lg\n", xmin, xmax, ymin, ymax);
  143.     }
  144.         
  145.     /* Now we are ready to take on some input values */
  146.     
  147.     x = (gmtdefs.xy_toggle) ? 1 : 0;    y = 1 - x;        /* Set up which columns have x and y */
  148.     
  149.     sprintf (format1, "%s\t%s\t%%s\n\0", gmtdefs.d_format, gmtdefs.d_format);
  150.     sprintf (format2, "%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format);
  151.     
  152.     if (n_files > 0)
  153.         nofile = FALSE;
  154.     else
  155.         n_files = 1;
  156.     n_args = (argc > 1) ? argc : 2;
  157.     
  158.     x_in_min = y_in_min = x_out_min = y_out_min = MAXDOUBLE;
  159.     x_in_max = y_in_max = x_out_max = y_out_max = -MAXDOUBLE;
  160.  
  161.     if (km == 2)
  162.         u_scale = (inverse) ? 1609.334 : 0.0006213711922;
  163.     else
  164.         u_scale = (inverse) ? 1000.0 : 0.001;
  165.     
  166.     for (fno = 1; !done && fno < n_args; fno++) {    /* Loop over input files, if any */
  167.         if (!nofile && argv[fno][0] == '-') continue;
  168.         
  169.         if (nofile) {    /* Just read standard input */
  170.             fp = stdin;
  171.             done = TRUE;
  172.             if (gmtdefs.verbose) fprintf (stderr, "mapproject: Reading from standard input\n");
  173.         }
  174.         else if ((fp = fopen (argv[fno], "r")) == NULL) {
  175.             fprintf (stderr, "mapproject: Cannot open file %s\n", argv[fno]);
  176.             continue;
  177.         }
  178.  
  179.         if (!nofile && gmtdefs.verbose) fprintf (stderr, "mapproject: Working on file %s\n", argv[fno]);
  180.         
  181.         if (gmtdefs.io_header) {
  182.             for (i = 0; i < gmtdefs.n_header_recs; i++) {
  183.                 fgets (data, 512, fp);
  184.                 if (first) printf ("%s", data);
  185.             }
  186.             first = FALSE;
  187.         }
  188.         if (multi_segments) {
  189.             fgets (line, 512, fp);
  190.             printf ("%s", line);
  191.         }
  192.     
  193.         if (inverse) {        /* Do inverse transformation */
  194.             while (fgets (line, 512, fp)) {    /* Read x,y */
  195.                 if (multi_segments && line[0] == EOL_flag) {
  196.                     printf ("%s", line);
  197.                     continue;
  198.                 }
  199.                 
  200.                 n_cols = sscanf (line, "%lf %lf%[^\n]\n", &xy_in[0], &xy_in[1], data);
  201.                 if (one_to_one) {    /* Convert from 1:1 scale */
  202.                     if (km) {
  203.                         xy_in[0] *= u_scale;
  204.                         xy_in[1] *= u_scale;
  205.                     }
  206.                     xy_in[0] *= project_info.x_scale;
  207.                     xy_in[1] *= project_info.y_scale;
  208.                 }
  209.                 if (map_center) {    /* Then correct so lower left corner is (0,0) */
  210.                     xy_in[0] += project_info.x0;
  211.                     xy_in[1] += project_info.y0;
  212.                 }
  213.                 xy_to_geo (&xy_out[0], &xy_out[1], xy_in[0], xy_in[1]);
  214.                 n_read++;
  215.                 if (suppress && map_outside (xy_out[0], xy_out[1])) continue;
  216.                 if (gmtdefs.verbose) {
  217.                     x_in_min = MIN (x_in_min, xy_in[0]);
  218.                     x_in_max = MAX (x_in_max, xy_in[0]);
  219.                     y_in_min = MIN (y_in_min, xy_in[1]);
  220.                     y_in_max = MAX (y_in_max, xy_in[1]);
  221.                     x_out_min = MIN (x_out_min, xy_out[0]);
  222.                     x_out_max = MAX (x_out_max, xy_out[0]);
  223.                     y_out_min = MIN (y_out_min, xy_out[1]);
  224.                     y_out_max = MAX (y_out_max, xy_out[1]);
  225.                 }
  226.                 if (n_cols == 2)
  227.                     printf (format2, xy_out[x], xy_out[y]);
  228.                 else
  229.                     printf (format1, xy_out[x], xy_out[y], data);
  230.                 n++;
  231.             }
  232.         }
  233.         else {        /* Do forward transformation */
  234.             while (fgets (line, 512, fp)) { /* Read lon,lat OR lat,lon */
  235.                 if (multi_segments && line[0] == EOL_flag) {
  236.                     printf ("%s", line);
  237.                     continue;
  238.                 }
  239.                 n_cols = sscanf (line, "%lf %lf%[^\n]", &xy_in[x], &xy_in[y], data);
  240.                 n_read++;
  241.                 if (suppress && map_outside (xy_in[0], xy_in[1])) continue;
  242.                 geo_to_xy (xy_in[0], xy_in[1], &xy_out[0], &xy_out[1]);
  243.                 if (map_center) {    /* Change origin from lower left to projection center */
  244.                     xy_out[0] -= project_info.x0;
  245.                     xy_out[1] -= project_info.y0;
  246.                 }
  247.                 if (one_to_one) {    /* Convert to 1:1 scale */
  248.                     xy_out[0] /= project_info.x_scale;
  249.                     xy_out[1] /= project_info.y_scale;
  250.                     if (km) {
  251.                         xy_out[0] *= u_scale;
  252.                         xy_out[1] *= u_scale;
  253.                     }
  254.                 }
  255.                 if (gmtdefs.verbose) {
  256.                     x_in_min = MIN (x_in_min, xy_in[0]);
  257.                     x_in_max = MAX (x_in_max, xy_in[0]);
  258.                     y_in_min = MIN (y_in_min, xy_in[1]);
  259.                     y_in_max = MAX (y_in_max, xy_in[1]);
  260.                     x_out_min = MIN (x_out_min, xy_out[0]);
  261.                     x_out_max = MAX (x_out_max, xy_out[0]);
  262.                     y_out_min = MIN (y_out_min, xy_out[1]);
  263.                     y_out_max = MAX (y_out_max, xy_out[1]);
  264.                 }
  265.                 if (n_cols == 2)
  266.                     printf (format2, xy_out[0], xy_out[1]);
  267.                 else
  268.                     printf (format1, xy_out[0], xy_out[1], data);
  269.                 n++;
  270.             }
  271.         }
  272.         if (fp != stdin) fclose (fp);
  273.     }
  274.     
  275.     if (gmtdefs.verbose && n_read > 0) {
  276.         fprintf (stderr, "mapproject: Input extreme values:  Xmin: %lg Xmax: %lg Ymin: %lg Ymax %lg\n",
  277.             x_in_min, x_in_max, y_in_min, y_in_max);
  278.         fprintf (stderr, "mapproject: Output extreme values: Xmin: %lg Xmax: %lg Ymin: %lg Ymax %lg\n",
  279.             x_out_min, x_out_max, y_out_min, y_out_max);
  280.         if (inverse)
  281.             fprintf (stderr, "mapproject: Mapped %d x-y pairs to lon-lat\n", n);
  282.         else
  283.             fprintf (stderr, "mapproject: Mapped %d lon-lat pairs to x-y\n", n);
  284.         if (suppress && n != n_read) fprintf (stderr, "mapproject: %d fell outside region\n", n_read - n);
  285.     }
  286.     
  287.     gmt_end (argc, argv);
  288. }
  289.