/*-------------------------------------------------------------------- * $Id: gmttrack.c,v 1.1.1.1.4.2 2002/02/27 17:41:10 pwessel Exp $ * * Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith * See README file for copying and redistribution conditions. *--------------------------------------------------------------------*/ /* * gmttrack will read *.gmt-files and and write PostScript code * that will create a navigation plot on a basemap using the * projection specified by the user. WESN must be specified on the * command line along with other options for scaling, anotation etc. * * Author: Paul Wessel * Date: 10-MAR-1988 * Revised: 22-NOV-1988 Color Version * 13-JUN-1991 v2.0 * 10-MAY-1995 v3.0b * 16-MAR-1999 v3.2, no longer -F * 30-JUN-1999 v3.3.1, Forgot to place INCH in ps_plotinit * No mention of -W in explanation. * 16-FEB-2000 v3.3.4, Only write to PS file via pslib * */ #include "gmt.h" #include "gmt_mgg.h" #define DST 0.25 #define ANSIZE 0.125 #define MPRDEG 111194.9 double heading (int rec, int n_records); struct GMTMGG_REC *record; int *distance; struct ANOT { int anot_int_dist; int anot_int_time; int tick_int_dist; int tick_int_time; } info; struct GMT_PEN pen; int get_anotinfo (char *args, struct ANOT *info); int gmttrack_outside (int lon, int lat, int w, int e, int s, int n); void anot_legname (double x, double y, int rec, int nrecs, char *text, int size, int w, int e, int s, int n); main (int argc, char **argv) { double west = 0.0, east = 0.0, south = 0.0, north = 0.0, lat, lon, angle; double last_lon, last_lat, anotsize = ANSIZE, factor = 0.001, dlon; double dx, dy, x, y, x0, y0, dist, start_dist = 0., stop_dist = 1.0E100; int i, iw, ie, is, in, anot_dist, tick_dist, anot_time, tick_time; int nfiles = 0, n_records, leg_year, rec, pen_type = 0, k, white[2], black[3]; int error = FALSE, first, last = FALSE, anot_name = FALSE, both = FALSE; int anot_tick = FALSE, draw_tick = FALSE, start_time = 0, stop_time = 2000000000; int this_julian, last_julian, year, month, day, hour, minute, second; int mon1, day1 = 0, year1, hour1, min1, mon2, day2 = 0, year2, hour2, min2; char agency[10], gmtfile[BUFSIZ], label[50], comment[BUFSIZ]; FILE *fp; struct GMTMGG_TIME *gmt; GMT_init_pen (&pen, GMT_PENWIDTH); for (i = 0; i < 3; i++) black[i] = 0, white[i] = 255; argc = GMT_begin (argc, argv); for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch(argv[i][1]) { case 'B': case 'P': case 'O': case 'K': case 'R': case 'J': case 'U': case 'V': case 'X': case 'x': case 'Y': case 'y': case 'c': case '\0': error += GMT_get_common_args (argv[i], &west, &east, &south, &north); break; case 'W': GMT_getpen (&argv[i][2], &pen); break; case 'A': anot_name = TRUE; if (argv[i][2]) anotsize = atof (&argv[i][2]) * GMT_u2u[GMT_PT][GMT_INCH]; break; case 'M': error += get_anotinfo (&argv[i][2], &info); break; case 'D': /* Assign start/stop times for sub-section */ if (argv[i][2] == 'a') { /* Start date */ sscanf(&argv[i][3], "%d/%d/%d/%d:%d", &mon1, &day1, &year1, &hour1, &min1); } else if (argv[i][2] == 'b') { /* Stop date */ sscanf(&argv[i][3], "%d/%d/%d/%d:%d", &mon2, &day2, &year2, &hour2, &min2); } else error = TRUE; break; case 'S': /* Assign start/stop position for sub-section */ if (argv[i][2] == 'a') /* Start position */ start_dist = atof(&argv[i][3]) * 1000.0; /* km -> m */ else if (argv[i][2] == 'b') /* Stop position */ stop_dist = atof(&argv[i][3]) * 1000.0; /* km -> m */ else error = TRUE; break; } } else nfiles++; } if (start_dist > stop_dist || start_time > stop_time) error = TRUE; if ((day1 > 0 && start_dist > 0.) || (day2 > 0 && stop_dist < 1.0e100)) error = TRUE; if (nfiles == 0) error = TRUE; if (GMT_check_rgb (pen.rgb)) error = TRUE; if (error) { fprintf(stderr,"usage: gmttrack leg(s) -R/// -J [-A[size]] [-B]\n"); fprintf(stderr," [-Da] [-Db] [-K] [-M] [-O] [-P] [-Sa]\n"); fprintf(stderr," [-Sb] [-U[dx/dy/][label]] [-V] [-W] [-X] [-Y] [-c\n\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf (stderr, " Give one or more names of .gmt files\n"); GMT_explain_option ('J'); GMT_explain_option ('R'); fprintf(stderr," OPTIONS:\n\n"); fprintf (stderr, " -A will Annotate legs when they enter the grid. is optional character size in points [9]\n"); GMT_explain_option ('B'); fprintf(stderr," -Da only plots from date (given as mm/dd/yyyy/hh:mm)\n"); fprintf(stderr," -Db only plots up to date (given as mm/dd/yyyy/hh:mm)\n"); GMT_explain_option ('K'); fprintf (stderr, " -M to put time/distance Marks on the track. E.g. a500ka24ht6h means (a)nnotate\n"); fprintf (stderr, " every 500 km (k) and 24 h(ours), with (t)ickmarks every 500 km and 6 (h)ours\n"); GMT_explain_option ('O'); GMT_explain_option ('P'); fprintf(stderr," -Sa only plots from dist (in km)\n"); fprintf(stderr," -Sb only plots up to dist (in km)\n"); GMT_explain_option ('U'); GMT_explain_option ('V'); fprintf (stderr, " -W sets track pen attributes [width = %lgp, color = (%d/%d/%d), texture = solid line].\n", pen.width, pen.rgb[0], pen.rgb[1], pen.rgb[2]); GMT_explain_option ('X'); GMT_explain_option ('c'); GMT_explain_option ('.'); exit (EXIT_FAILURE); } GMT_put_history (argc, argv); /* Update .gmtdefaults */ gmtmggpath_init(); if (east < west) west -= 360.0; GMT_map_setup (west, east, south, north); ps_plotinit (NULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin, gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies, gmtdefs.dpi, GMT_INCH, gmtdefs.paper_width, gmtdefs.page_rgb, GMT_epsinfo (argv[0])); GMT_echo_command (argc, argv); if (gmtdefs.unix_time) GMT_timestamp (argc, argv); ps_setpaint (gmtdefs.basemap_frame_rgb); GMT_setpen (&pen); distance = (int *) GMT_memory (VNULL, (size_t)1, sizeof (int), "gmttrack"); record = (struct GMTMGG_REC *) GMT_memory (VNULL, (size_t)1, sizeof (struct GMTMGG_REC), "gmttrack"); ps_setfont (7); /* Times Italic */ GMT_map_clip_on (GMT_no_rgb, 3); iw = irint (west * 1.0e6); ie = irint (east * 1.0e6); is = irint (south * 1.0e6); in = irint (north * 1.0e6); both = (info.anot_int_time && info.anot_int_dist); ps_setpaint (pen.rgb); for (i = 1; i < argc; i++) { /* Loop over arguments and open/plot the gmtfiles */ if (argv[i][0] == '-') continue; if (gmtmggpath_func(gmtfile,argv[i])) { fprintf(stderr,"gmttrack : Cannot find leg %s\n", argv[i]); continue; } if ((fp = fopen(gmtfile, "rb")) == NULL) { fprintf(stderr, "gmttrack: Could not open file %s\n", gmtfile); continue; } /* Read first record of file containing start-year, n_records and info */ if (fread((void *)(&leg_year), (size_t)4, (size_t)1, fp) != 1) { fprintf(stderr,"gmttrack: Error while reading first year\n"); exit(EXIT_FAILURE); } if (fread((void *)(&n_records), (size_t)4, (size_t)1, fp) != 1) { fprintf(stderr,"gmttrack: Error while reading no of records\n"); exit(EXIT_FAILURE); } if (fread((void *)agency, (size_t)10, (size_t)1, fp) != 1) { fprintf(stderr,"gmttrack: Error while reading info-header\n"); exit(EXIT_FAILURE); } distance = (int *) GMT_memory ((void *)distance, (size_t)n_records, sizeof (int), "gmttrack"); record = (struct GMTMGG_REC *) GMT_memory ((void *)record, (size_t)n_records, sizeof (struct GMTMGG_REC), "gmttrack"); gmt = gmtmgg_init(leg_year); /* Initialize gmt_structure */ last_julian = -1; /* Decode date to time in sec if needed */ if (day1 > 0) gmtmgg_time (&start_time, year1, mon1, day1, hour1, min1, 0, gmt); if (day2 > 0) gmtmgg_time (&stop_time, year2, mon2, day2, hour2, min2, 0, gmt); /* Read cruise data */ for (rec = k = 0; rec < n_records; rec++) { if (fread((void *)(&record[k]), (size_t)18, (size_t)1, fp) != 1) { fprintf(stderr,"gmttrack: Error reading data record no %d\n",rec); exit(EXIT_FAILURE); } lon = record[k].lon * MDEG2DEG; lat = record[k].lat * MDEG2DEG; if (rec == 0) { dist = 0; last_lon = lon; last_lat = lat; } else { dlon = fabs (lon - last_lon); if (dlon > 180.0) dlon = 360.0 - dlon; dx = dlon * cos (0.5 * D2R * (lat+last_lat)); dy = (lat - last_lat); last_lon = lon; last_lat = lat; dist += hypot (dx, dy) * MPRDEG; /* Dist in meters */ } distance[k] = irint (dist); if (dist < start_dist || dist > stop_dist) continue; if (record[k].time < start_time || record[k].time > stop_time) continue; k++; } fclose (fp); n_records = k; if (gmtdefs.verbose) fprintf (stderr, "gmttrack: Plotting %s\n", agency); sprintf (comment, "Tracking %s", agency); ps_comment (comment); /* First draw the track line, clip segments outside the area */ first = TRUE; last = FALSE; for (rec = 0; rec < n_records; rec++) { if (gmttrack_outside (record[rec].lon, record[rec].lat, iw, ie, is, in)) { if (last && rec < (n_records-1)) { lon = record[rec].lon * 1.0e-6; lat = record[rec].lat * 1.0e-6; GMT_geo_to_xy (lon, lat, &x, &y); ps_plot (x, y, 2); last = FALSE; } first = TRUE; continue; } lon = record[rec].lon * 1.0e-6; lat = record[rec].lat * 1.0e-6; GMT_geo_to_xy (lon, lat, &x, &y); if (first) { if (rec > 0) { lon = record[rec].lon * 1.0e-6; lat = record[rec].lat * 1.0e-6; GMT_geo_to_xy (lon, lat, &x0, &y0); ps_plot (x0, y0, 3); } else ps_plot (x, y, 3); if (anot_name) anot_legname (x, y, rec, n_records, argv[i], irint (GMT_u2u[GMT_INCH][GMT_PT] * 1.25 * anotsize), iw ,ie, is, in); first = FALSE; if (info.anot_int_dist > 0) anot_dist = (distance[rec] / info.anot_int_dist + 1) * info.anot_int_dist; if (info.tick_int_dist > 0) tick_dist = (distance[rec] / info.tick_int_dist + 1) * info.tick_int_dist; if (info.anot_int_time > 0) anot_time = (record[rec].time / info.anot_int_time + 1) * info.anot_int_time; if (info.tick_int_time > 0) tick_time = (record[rec].time / info.tick_int_time + 1) * info.tick_int_time; } ps_plot (x, y, 2); last = TRUE; /* See if we need to anotate/tick the trackline for time/km marks */ if (info.anot_int_time && (record[rec].time >= anot_time)) { anot_time += info.anot_int_time; anot_tick = 1; } if (info.anot_int_dist && (distance[rec] >= anot_dist)) { anot_dist += info.anot_int_dist; anot_tick += 2; } if (info.tick_int_time && (record[rec].time >= tick_time)) { tick_time += info.tick_int_time; draw_tick = 1; } if (info.tick_int_dist && (distance[rec] >= tick_dist)) { tick_dist += info.tick_int_dist; draw_tick += 2; } if (anot_tick) { angle = heading (rec, n_records); if (angle < 0.0) angle += 90.0; else angle -= 90.0; if (anot_tick & 1) { /* Time mark */ this_julian = gmtmgg_date (anot_time - info.anot_int_time,&year,&month,&day,&hour,&minute,&second,gmt); ps_command ("S"); if (this_julian != last_julian) { sprintf (label, " %.2d:%.2d %d/%d", hour, minute, month, day); ps_circle (x, y, 0.2*anotsize, black, TRUE); ps_setfont (7); } else { sprintf (label, " %.2d:%.2d", hour, minute); ps_circle (x, y, 0.2*anotsize, white, TRUE); } ps_setpaint (gmtdefs.basemap_frame_rgb); ps_text (x, y, (int)rint (GMT_u2u[GMT_INCH][GMT_PT] * anotsize), label, angle, 5, 0); if (this_julian != last_julian) ps_setfont (6); last_julian = this_julian; } if (anot_tick & 2) { /* Distance mark */ ps_command ("S"); sprintf (label, "%d km ", (int)((anot_dist - info.anot_int_dist) * factor)); ps_square (x, y, 0.4*anotsize, gmtdefs.basemap_frame_rgb, FALSE); ps_setpaint (gmtdefs.basemap_frame_rgb); ps_text (x, y, irint (GMT_u2u[GMT_INCH][GMT_PT] * anotsize), label, angle, 7, 0); } } if (both && !(anot_tick & 1) && (draw_tick & 1)) { ps_command ("S"); ps_circle (x, y, 0.2*anotsize, gmtdefs.basemap_frame_rgb, FALSE); } if (both && !(anot_tick & 2) && (draw_tick & 2)) { ps_command ("S"); ps_square (x, y, 0.4*anotsize, gmtdefs.basemap_frame_rgb, FALSE); } if (draw_tick) { ps_command ("S"); ps_setpaint (gmtdefs.basemap_frame_rgb); ps_cross (x, y, 0.4*anotsize); } if (anot_tick || draw_tick) { ps_setpaint (pen.rgb); ps_plot (x, y, 3); anot_tick = draw_tick = FALSE; } } } GMT_map_clip_off(); if (pen_type > 0) ps_setdash (NULL, 0); if (frame_info.plot) { ps_setpaint (gmtdefs.basemap_frame_rgb); GMT_map_basemap (); } ps_plotend (gmtdefs.last_page); GMT_free ((void *)record); GMT_free ((void *)distance); GMT_end (argc, argv); } double heading (int rec, int n_records) { int i1, i2; double angle, lon, lat, x1, x2, y1, y2; i1 = rec - 25; if (i1 < 0) i1 = 0; i2 = i1 + 25; if (i2 > (n_records-1)) i2 = n_records - 1; lon = record[i2].lon * 1.0e-6; lat = record[i2].lat * 1.0e-6; GMT_geo_to_xy (lon, lat, &x2, &y2); lon = record[i1].lon * 1.0e-6; lat = record[i1].lat * 1.0e-6; GMT_geo_to_xy (lon, lat, &x1, &y1); angle = atan2 (y2 - y1, x2 - x1) * R2D; if (angle > 90.) angle -= 180; else if (angle < -90.0) angle += 180.0; return (angle); } int get_anotinfo (char *args, struct ANOT *info) { int i1, i2, error = FALSE; int flag1, flag2, type; double value; info->anot_int_dist = info->tick_int_dist = 0; info->anot_int_time = info->tick_int_time = 0; i1 = 0; while (args[i1]) { flag1 = 'a'; if (isalpha (args[i1])) { flag1 = args[i1]; if (isupper(flag1)) flag1 = tolower (flag1); i1++; } i2 = i1; while (args[i2] && strchr ("athkm", (int)args[i2]) == NULL) i2++; value = atof (&args[i1]); flag2 = args[i2]; if (isupper(flag2)) flag2 = tolower (flag2); if (flag2 == 'h') { /* Hours */ value *= 3600; type = 't'; } else if (flag2 == 'k') { /* kilometers */ value *= 1000; type = 'd'; } else if (flag2 == 'm') { /* Minutes */ value *= 60; type = 't'; } else /* Default is seconds */ type = 't'; i2++; if (flag1 == 'a') { /* Anotation interval */ if (type == 'd') /* Distance */ info->anot_int_dist = (int)value; else info->anot_int_time = (int)value; } else { /* Tickmark interval */ if (type == 'd') /* Distance */ info->tick_int_dist = (int)value; else info->tick_int_time = (int)value; } i1 = i2; } if (info->anot_int_dist <= 0 && info->tick_int_dist <= 0 && info->anot_int_time <= 0 && info->tick_int_time <= 0) error = TRUE; if (info->anot_int_dist <= 0) info->anot_int_dist = info->tick_int_dist; else if (info->tick_int_dist <= 0) info->tick_int_dist = info->anot_int_dist; if (info->anot_int_time <= 0) info->anot_int_time = info->tick_int_time; else if (info->tick_int_time <= 0) info->tick_int_time = info->anot_int_time; return (error); } int gmttrack_outside (int lon, int lat, int w, int e, int s, int n) { if (lat < s || lat > n) return (TRUE); lon -= 360000000; while (lon < w) lon += 360000000; if (lon > e) return (TRUE); return (FALSE); } void anot_legname (double x, double y, int rec, int nrecs, char *text, int size, int w, int e, int s, int n) { int just; double angle; if (rec == 0) { angle = 0.0; just = 6; } else { angle = heading (rec, nrecs); if (record[rec-1].lat < s) just = (angle >= 0.0) ? 1 : 3; else if (record[rec-1].lat > n) just = (angle >= 0.0) ? 11 : 9; else if (record[rec-1].lon < w) just = (angle >= 0.0) ? 9 : 1; else just = (angle >= 0.0) ? 3 : 11; } ps_setfont (5); ps_setpaint (gmtdefs.basemap_frame_rgb); ps_text (x, y, size, text, angle, just, 0); ps_setpaint (pen.rgb); ps_plot (x, y, 3); ps_setfont (7); }