/*-------------------------------------------------------------------- * $Id: grdpmodeler.c 12822 2014-01-31 23:39:56Z remko $ * * Copyright (c) 1999-2014 by P. Wessel * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by * the Free Software Foundation; version 3 or any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * Contact info: www.soest.hawaii.edu/pwessel *--------------------------------------------------------------------*/ /* * grdpmodeler will read an age grid file and a plate motion model and * evaluates chosen attributes such as speed, direction, lon/lat origin * of the grid location. * * Author: Paul Wessel * Date: 1-NOV-2010 * Ver: 5 */ #define THIS_MODULE_NAME "grdpmodeler" #define THIS_MODULE_LIB "spotter" #define THIS_MODULE_PURPOSE "Evaluate a plate model on a geographic grid" #include "spotter.h" #define GMT_PROG_OPTIONS "-:RVbhi" #define N_PM_ITEMS 8 #define PM_RATE 0 #define PM_AZIM 1 #define PM_OMEGA 2 #define PM_DLON 3 #define PM_LON 4 #define PM_DLAT 5 #define PM_LAT 6 #define PM_DIST 7 struct GRDROTATER_CTRL { /* All control options for this program (except common args) */ /* active is true if the option has been activated */ struct In { bool active; char *file; } In; struct E { /* -Erotfile */ bool active; char *file; } E; struct F { /* -Fpolfile */ bool active; char *file; } F; struct G { /* -Goutfile */ bool active; char *file; } G; struct I { /* -Idx[/dy] */ bool active; double inc[2]; } I; struct S { /* -Sa|d|r|w|x|y|X|Y */ bool active; unsigned int mode; } S; struct T { /* -T */ bool active; double value; } T; }; void *New_grdpmodeler_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ struct GRDROTATER_CTRL *C; C = GMT_memory (GMT, NULL, 1, struct GRDROTATER_CTRL); /* Initialize values whose defaults are not 0/false/NULL */ return (C); } void Free_grdpmodeler_Ctrl (struct GMT_CTRL *GMT, struct GRDROTATER_CTRL *C) { /* Deallocate control structure */ if (!C) return; if (C->In.file) free (C->In.file); if (C->E.file) free (C->E.file); if (C->F.file) free (C->F.file); if (C->G.file) free (C->G.file); GMT_free (GMT, C); } int GMT_grdpmodeler_usage (struct GMTAPI_CTRL *API, int level) { GMT_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_NAME, THIS_MODULE_PURPOSE); if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); GMT_Message (API, GMT_TIME_NONE, "usage: grdpmodeler -E -G [-F]\n"); GMT_Message (API, GMT_TIME_NONE, "\t[%s] [%s] [-Sa|d|r|w|x|y|X|Y]\n\t[-T] [%s] [%s]\n\t[%s] [%s]\n\t[%s]\n\n", GMT_Id_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_b_OPT, GMT_h_OPT, GMT_i_OPT, GMT_r_OPT); if (level == GMT_SYNOPSIS) return (EXIT_FAILURE); GMT_Message (API, GMT_TIME_NONE, "\t is a gridded data file in geographic coordinates with crustal ages.\n"); spotter_rot_usage (API, 'E'); GMT_Message (API, GMT_TIME_NONE, "\t-G Set output filename with the model predictions.\n"); GMT_Message (API, GMT_TIME_NONE, "\n\tOPTIONS:\n"); GMT_Message (API, GMT_TIME_NONE, "\t-F Specify a multi-segment closed polygon file that describes the area\n"); GMT_Message (API, GMT_TIME_NONE, "\t of the grid to work on [Default works on the entire grid].\n"); GMT_Option (API, "Rg"); GMT_Message (API, GMT_TIME_NONE, "\t-S Select a model prediction as a function of crustal age. Choose among:\n"); GMT_Message (API, GMT_TIME_NONE, "\t a : Plate spreading azimuth.\n"); GMT_Message (API, GMT_TIME_NONE, "\t d : Distance to origin of crust in km.\n"); GMT_Message (API, GMT_TIME_NONE, "\t r : Plate motion rate in mm/yr or km/Myr.\n"); GMT_Message (API, GMT_TIME_NONE, "\t w : Rotation rate in degrees/Myr.\n"); GMT_Message (API, GMT_TIME_NONE, "\t x : Change in longitude since formation.\n"); GMT_Message (API, GMT_TIME_NONE, "\t y : Change in latitude since formation.\n"); GMT_Message (API, GMT_TIME_NONE, "\t X : Longitude at origin of crust.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Y : Latitude at origin of crust.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-T Set fixed time of reconstruction to override age grid.\n"); GMT_Option (API, "V,bi2,h,i,r,."); return (EXIT_FAILURE); } int GMT_grdpmodeler_parse (struct GMT_CTRL *GMT, struct GRDROTATER_CTRL *Ctrl, struct GMT_OPTION *options) { /* This parses the options provided to grdpmodeler and sets parameters in CTRL. * Any GMT common options will override values set previously by other commands. * It also replaces any file names specified as input or output with the data ID * returned when registering these sources/destinations with the API. */ unsigned int n_errors = 0, n_files = 0; struct GMT_OPTION *opt = NULL; for (opt = options; opt; opt = opt->next) { switch (opt->option) { case '<': /* Input files */ if (n_files++ > 0) break; if ((Ctrl->In.active = GMT_check_filearg (GMT, '<', opt->arg, GMT_IN))) Ctrl->In.file = strdup (opt->arg); else n_errors++; break; /* Supplemental parameters */ case 'E': /* File with stage poles */ if ((Ctrl->E.active = GMT_check_filearg (GMT, 'E', opt->arg, GMT_IN))) Ctrl->E.file = strdup (opt->arg); else n_errors++; break; case 'F': if ((Ctrl->F.active = GMT_check_filearg (GMT, 'F', opt->arg, GMT_IN))) Ctrl->F.file = strdup (opt->arg); else n_errors++; break; case 'G': if ((Ctrl->G.active = GMT_check_filearg (GMT, 'G', opt->arg, GMT_OUT))) Ctrl->G.file = strdup (opt->arg); else n_errors++; break; case 'I': Ctrl->I.active = true; if (GMT_getinc (GMT, opt->arg, Ctrl->I.inc)) { GMT_inc_syntax (GMT, 'I', 1); n_errors++; } break; case 'S': Ctrl->S.active = true; switch (opt->arg[0]) { case 'a': /* Plate spreading azimuth */ Ctrl->S.mode = PM_AZIM; break; case 'd': /* Distance from point to origin at ridge */ Ctrl->S.mode = PM_DIST; break; case 'r': /* Plate spreading rate */ Ctrl->S.mode = PM_RATE; break; case 'w': /* Plate rotation omega */ Ctrl->S.mode = PM_OMEGA; break; case 'x': /* Change in longitude since origin */ Ctrl->S.mode = PM_DLON; break; case 'y': /* Change in latitude since origin */ Ctrl->S.mode = PM_DLAT; break; case 'X': /* Plate longitude at crust origin */ Ctrl->S.mode = PM_LON; break; case 'Y': /* Plate latitude at crust origin */ Ctrl->S.mode = PM_LAT; break; default: n_errors++; break; } break; case 'T': Ctrl->T.active = true; Ctrl->T.value = atof (opt->arg); break; default: /* Report bad options */ n_errors += GMT_default_error (GMT, opt->option); break; } } if (!Ctrl->In.file) { /* Must have -R -I [-r] */ n_errors += GMT_check_condition (GMT, !GMT->common.R.active && !Ctrl->I.active, "Syntax error: Must specify input file or -R -I [-r]\n"); } else { /* Must not have -I -r */ n_errors += GMT_check_condition (GMT, Ctrl->I.active || GMT->common.r.active, "Syntax error: Cannot specify input file AND -R -r\n"); } n_errors += GMT_check_condition (GMT, !Ctrl->G.file, "Syntax error -G: Must specify output file\n"); n_errors += GMT_check_condition (GMT, !Ctrl->E.active, "Syntax error: Must specify -E\n"); n_errors += GMT_check_condition (GMT, !Ctrl->S.active, "Syntax error: Must specify -S\n"); n_errors += GMT_check_condition (GMT, Ctrl->T.value < 0.0, "Syntax error -T: Must specify positive age.\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_OK); } #define bailout(code) {GMT_Free_Options (mode); return (code);} #define Return(code) {Free_grdpmodeler_Ctrl (GMT, Ctrl); GMT_end_module (GMT, GMT_cpy); bailout (code);} int GMT_grdpmodeler (void *V_API, int mode, void *args) { unsigned int col, row, inside, stage, n_stages, registration; int retval, error = 0; uint64_t node, seg; double lon, lat, d, value = 0.0, t_max = 0.0, age, wesn[4], inc[2], *grd_x = NULL, *grd_y = NULL, *grd_yc = NULL; char *quantity[N_PM_ITEMS] = { "velocity", "azimuth", "rotation rate", "longitude displacement", \ "longitude", "latitude displacement", "latitude", "distance displacement"}; struct GMT_DATASET *D = NULL; struct GMT_DATATABLE *pol = NULL; struct EULER *p = NULL; /* Pointer to array of stage poles */ struct GMT_OPTION *ptr = NULL; struct GMT_GRID *G_age = NULL, *G_mod = NULL; struct GRDROTATER_CTRL *Ctrl = NULL; struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL; struct GMT_OPTION *options = NULL; struct GMTAPI_CTRL *API = GMT_get_API_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */ /*----------------------- Standard module initialization and parsing ----------------------*/ if (API == NULL) return (GMT_NOT_A_SESSION); if (mode == GMT_MODULE_PURPOSE) return (GMT_grdpmodeler_usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */ options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */ if (!options || options->option == GMT_OPT_USAGE) bailout (GMT_grdpmodeler_usage (API, GMT_USAGE)); /* Return the usage message */ if (options->option == GMT_OPT_SYNOPSIS) bailout (GMT_grdpmodeler_usage (API, GMT_SYNOPSIS)); /* Return the synopsis */ /* Parse the command-line arguments */ GMT = GMT_begin_module (API, THIS_MODULE_LIB, THIS_MODULE_NAME, &GMT_cpy); /* Save current state */ if ((ptr = GMT_Find_Option (API, 'f', options)) == NULL) GMT_parse_common_options (GMT, "f", 'f', "g"); /* Did not set -f, implicitly set -fg */ if (GMT_Parse_Common (API, GMT_PROG_OPTIONS, options)) Return (API->error); Ctrl = New_grdpmodeler_Ctrl (GMT); /* Allocate and initialize a new control structure */ if ((error = GMT_grdpmodeler_parse (GMT, Ctrl, options))) Return (error); /*---------------------------- This is the grdpmodeler main code ----------------------------*/ /* Check limits and get data file */ if (Ctrl->In.file) { /* Gave an age grid */ if ((G_age = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_HEADER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) { /* Get header only */ Return (API->error); } GMT_memcpy (wesn, (GMT->common.R.active ? GMT->common.R.wesn : G_age->header->wesn), 4, double); if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_DATA_ONLY, wesn, Ctrl->In.file, G_age) == NULL) { Return (API->error); /* Get header only */ } GMT_memcpy (inc, G_age->header->inc, 2, double); /* Use same increment for output grid */ registration = G_age->header->registration; GMT_memcpy (GMT->common.R.wesn, G_age->header->wesn, 4, double); GMT->common.R.active = true; } else { GMT_memcpy (inc, Ctrl->I.inc, 2, double); registration = GMT->common.r.registration; } if (Ctrl->F.active) { /* Read the user's clip polygon file */ if ((D = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_READ_NORMAL, NULL, Ctrl->F.file, NULL)) == NULL) { Return (API->error); } pol = D->table[0]; /* Since it is a single file */ GMT_Report (API, GMT_MSG_VERBOSE, "Restrict evalution to within polygons in file %s\n", Ctrl->F.file); } n_stages = spotter_init (GMT, Ctrl->E.file, &p, false, false, 0, &t_max); for (stage = 0; stage < n_stages; stage++) { if (p[stage].omega < 0.0) { /* Ensure all stages have positive rotation angles */ p[stage].omega = -p[stage].omega; p[stage].lat = -p[stage].lat; p[stage].lon += 180.0; if (p[stage].lon > 360.0) p[stage].lon -= 360.0; } } if (Ctrl->T.active && Ctrl->T.value > t_max) { GMT_Report (API, GMT_MSG_NORMAL, "Requested a fixed reconstruction time outside range of rotation table\n"); GMT_free (GMT, p); Return (EXIT_FAILURE); } if ((G_mod = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_GRID_ALL, NULL, NULL, inc, \ registration, GMT_NOTSET, Ctrl->G.file)) == NULL) Return (API->error); grd_x = GMT_memory (GMT, NULL, G_mod->header->nx, double); grd_y = GMT_memory (GMT, NULL, G_mod->header->ny, double); grd_yc = GMT_memory (GMT, NULL, G_mod->header->ny, double); /* Precalculate node coordinates in both degrees and radians */ for (row = 0; row < G_mod->header->ny; row++) { grd_y[row] = GMT_grd_row_to_y (GMT, row, G_mod->header); grd_yc[row] = GMT_lat_swap (GMT, grd_y[row], GMT_LATSWAP_G2O); } for (col = 0; col < G_mod->header->nx; col++) grd_x[col] = GMT_grd_col_to_x (GMT, col, G_mod->header); /* Loop over all nodes in the new rotated grid and find those inside the reconstructed polygon */ GMT_Report (API, GMT_MSG_VERBOSE, "Evalute model prediction grid\n"); GMT_init_distaz (GMT, (Ctrl->S.mode == PM_DIST) ? 'k' : 'd', GMT_GREATCIRCLE, GMT_MAP_DIST); /* Great circle distances in degrees, or km if -Sd */ if (Ctrl->S.mode == PM_DLON) GMT->current.io.geo.range = GMT_IS_M180_TO_P180_RANGE; /* Need +- around 0 here */ GMT_grd_loop (GMT, G_mod, row, col, node) { G_mod->data[node] = GMT->session.f_NaN; if (Ctrl->F.active) { for (seg = inside = 0; seg < pol->n_segments && !inside; seg++) { /* Use degrees since function expects it */ if (GMT_polygon_is_hole (pol->segment[seg])) continue; /* Holes are handled within GMT_inonout */ inside = (GMT_inonout (GMT, grd_x[col], grd_y[row], pol->segment[seg]) > 0); } if (!inside) continue; /* Outside the polygon(s) */ } /* Here we are inside; get the coordinates and rotate back to original grid coordinates */ age = (Ctrl->T.active) ? Ctrl->T.value : G_age->data[node]; if (GMT_is_dnan (age)) continue; /* No crustal age */ if ((retval = spotter_stage (GMT, age, p, n_stages)) < 0) continue; /* Outside valid stage rotation range */ stage = retval; switch (Ctrl->S.mode) { case PM_RATE: /* Compute plate motion speed at this point in time/space */ d = GMT_distance (GMT, grd_x[col], grd_yc[row], p[stage].lon, p[stage].lat); value = sind (d) * p[stage].omega * GMT->current.proj.DIST_KM_PR_DEG; /* km/Myr or mm/yr */ break; case PM_AZIM: /* Compute plate motion direction at this point in time/space */ value = GMT_az_backaz (GMT, grd_x[col], grd_yc[row], p[stage].lon, p[stage].lat, false) - 90.0; GMT_lon_range_adjust (GMT->current.io.geo.range, &value); break; case PM_OMEGA: /* Compute plate rotation rate omega */ value = p[stage].omega; /* degree/Myr */ break; case PM_DLAT: /* Compute latitude where this point was formed in the model */ lon = grd_x[col] * D2R; lat = grd_yc[row] * D2R; (void)spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL); value = grd_y[row] - GMT_lat_swap (GMT, lat * R2D, GMT_LATSWAP_O2G); /* Convert back to geodetic */ break; case PM_LAT: /* Compute latitude where this point was formed in the model */ lon = grd_x[col] * D2R; lat = grd_yc[row] * D2R; (void)spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL); value = GMT_lat_swap (GMT, lat * R2D, GMT_LATSWAP_O2G); /* Convert back to geodetic */ break; case PM_DLON: /* Compute latitude where this point was formed in the model */ lon = grd_x[col] * D2R; lat = grd_yc[row] * D2R; (void)spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL); value = grd_x[col] - lon * R2D; if (fabs (value) > 180.0) value = copysign (360.0 - fabs (value), -value); break; case PM_LON: /* Compute latitude where this point was formed in the model */ lon = grd_x[col] * D2R; lat = grd_yc[row] * D2R; (void)spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL); value = lon * R2D; break; case PM_DIST: /* Compute distance between node and point of origin at ridge */ lon = grd_x[col] * D2R; lat = grd_yc[row] * D2R; (void)spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL); value = GMT_distance (GMT, grd_x[col], grd_yc[row], lon * R2D, lat * R2D); break; } G_mod->data[node] = (float)value; } /* Now write model prediction grid */ GMT_Report (API, GMT_MSG_VERBOSE, "Write predicted grid\n"); strcpy (G_mod->header->x_units, "degrees_east"); strcpy (G_mod->header->y_units, "degrees_north"); switch (Ctrl->S.mode) { case PM_RATE: /* Compute plate motion speed at this point in time/space */ strcpy (G_mod->header->z_units, "mm/yr"); break; case PM_AZIM: /* Compute plate motion direction at this point in time/space */ strcpy (G_mod->header->z_units, "degree"); break; case PM_OMEGA: /* Compute plate rotation rate omega */ strcpy (G_mod->header->z_units, "degree/Myr"); break; case PM_DLAT: /* Difference in latitude relative to where this point was formed in the model */ case PM_LAT: /* Latitude where this point was formed in the model */ strcpy (G_mod->header->z_units, "degrees_north"); break; case PM_DLON: /* Difference in longitude relative to where this point was formed in the model */ case PM_LON: /* Longitude where this point was formed in the model */ strcpy (G_mod->header->z_units, "degrees_east"); break; case PM_DIST: /* Distance to origin in km */ strcpy (G_mod->header->z_units, "km"); break; } sprintf (G_mod->header->remark, "Plate Model predictions of %s for model %s", quantity[Ctrl->S.mode], Ctrl->E.file); if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, G_mod)) Return (API->error); if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_ALL, NULL, Ctrl->G.file, G_mod) != GMT_OK) { Return (API->error); } GMT_free (GMT, grd_x); GMT_free (GMT, grd_y); GMT_free (GMT, grd_yc); GMT_free (GMT, p); GMT_Report (API, GMT_MSG_VERBOSE, "Done!\n"); Return (GMT_OK); }