diff --git a/admin/add_windows_cpack.txt b/admin/add_windows_cpack.txt index 2ac499e3bc0..ad246d661c1 100644 --- a/admin/add_windows_cpack.txt +++ b/admin/add_windows_cpack.txt @@ -85,3 +85,13 @@ if (GMTMEX_PATH) DESTINATION ${GMT_BINDIR} COMPONENT Runtime) endif (GMTMEX_PATH) + +install (PROGRAMS + ${VCRUNTIME_PATH}/vcruntime140.dll + ${VCRUNTIME_PATH}/vcruntime140_1.dll + ${VCRUNTIME_PATH}/msvcp140.dll + ${VCRUNTIME_PATH}/msvcp140_1.dll + ${VCRUNTIME_PATH}/msvcp140_2.dll + ${VCRUNTIME_PATH}/msvcp140_codecvt_ids.dll + DESTINATION ${GMT_BINDIR} + COMPONENT Runtime) diff --git a/doc/rst/source/modules-classic.rst b/doc/rst/source/modules-classic.rst index 152145322a5..d4a0a1c2ad4 100644 --- a/doc/rst/source/modules-classic.rst +++ b/doc/rst/source/modules-classic.rst @@ -690,6 +690,10 @@ seis +----------------------------------+--------------------+ | :doc:`/supplements/seis/pssac` | |pssac_purpose| | +----------------------------------+--------------------+ +| :doc:`/supplements/seis/shake` | |shake_purpose| | ++--------------------------------+----------------------+ +| :doc:`/supplements/seis/vs30` | |vs30_purpose| | ++--------------------------------+----------------------+ spotter ------- diff --git a/doc/rst/source/modules.rst b/doc/rst/source/modules.rst index 11ab459dd2f..bfd740f3a8b 100644 --- a/doc/rst/source/modules.rst +++ b/doc/rst/source/modules.rst @@ -160,6 +160,8 @@ All modules are requested via a call to the :doc:`gmt` program. supplements/seis/meca supplements/seis/polar supplements/seis/sac + supplements/seis/shake + supplements/seis/vs30 supplements/spotter/backtracker supplements/spotter/gmtpmodeler supplements/spotter/grdpmodeler @@ -721,6 +723,10 @@ seis +--------------------------------+------------------+ | :doc:`/supplements/seis/sac` | |sac_purpose| | +--------------------------------+------------------+ +| :doc:`/supplements/seis/shake` | |shake_purpose| | ++--------------------------------+------------------+ +| :doc:`/supplements/seis/vs30` | |vs30_purpose| | ++--------------------------------+------------------+ spotter ------- diff --git a/doc/rst/source/reference/supplemental-packages.rst b/doc/rst/source/reference/supplemental-packages.rst index 72c3a719e6d..768ae1dc1fd 100644 --- a/doc/rst/source/reference/supplemental-packages.rst +++ b/doc/rst/source/reference/supplemental-packages.rst @@ -87,10 +87,12 @@ seis: Seismology This package contains the programs :doc:`coupe `, :doc:`meca `, -:doc:`polar `, and -:doc:`sac ` which are used by seismologists +:doc:`polar `, +:doc:`sac `, +:doc:`sac `, and +:doc:`sac `, which are used by seismologists for plotting focal mechanisms (including cross-sections -and polarities) and SAC files. +and polarities), compute Vs30 velocities, intensity maps and SAC files. The coupe, meca, and polar were developed by Kurt Feigl and Genevieve Patau, while Dongdong Tian added sac; the package is now maintained by the GMT team. diff --git a/doc/rst/source/supplements/module_supplements_purpose.rst_ b/doc/rst/source/supplements/module_supplements_purpose.rst_ index 08584edd66d..22de57ec8e5 100644 --- a/doc/rst/source/supplements/module_supplements_purpose.rst_ +++ b/doc/rst/source/supplements/module_supplements_purpose.rst_ @@ -12,8 +12,6 @@ .. |img2grd_purpose| replace:: Extract a subset from an img file in Mercator or Geographic format -.. |terrain_filter_purpose| replace:: Callable function to do image texture - .. |mgd77convert_purpose| replace:: Convert MGD77 data to other formats .. |mgd77header_purpose| replace:: Create MGD77 headers from A77 files @@ -98,11 +96,9 @@ .. |rotsmoother_purpose| replace:: Get mean rotations and covariance matrices from set of finite rotations -.. |grdbarb_purpose| replace:: Plot wind barb field from two component grids - -.. |barb_purpose| replace:: Plot wind barbs in 2-D and 3-D +.. |shake_purpose| replace:: Compute Peak Ground Acceleration/Velocity and Intensity -.. |psbarb_purpose| replace:: Plot wind barbs in 2-D and 3-D +.. |vs30_purpose| replace:: Convert topographic slope to Vs30 velocities .. |x2sys_binlist_purpose| replace:: Create bin index listing from track data files diff --git a/src/seis/CMakeLists.txt b/src/seis/CMakeLists.txt index ca37fcd2f97..644dc756219 100644 --- a/src/seis/CMakeLists.txt +++ b/src/seis/CMakeLists.txt @@ -27,6 +27,6 @@ set (SUPPL_NAME seis) set (SUPPL_HEADERS meca.h meca_symbol.h utilmeca.h seis_defaults.h sacio.h) AUX_SOURCE_DIRECTORY (longopt SUPPL_LONG_OPT_H) -set (SUPPL_PROGS_SRCS psmeca.c pspolar.c pscoupe.c pssac.c ${SUPPL_LONG_OPT_H}) +set (SUPPL_PROGS_SRCS psmeca.c pspolar.c pscoupe.c pssac.c shake.c vs30.c ${SUPPL_LONG_OPT_H}) set (SUPPL_LIB_SRCS ${SUPPL_PROGS_SRCS} utilmeca.c sacio.c) set (SUPPL_EXAMPLE_FILES README.seis) diff --git a/src/seis/longopt/shake_inc.h b/src/seis/longopt/shake_inc.h new file mode 100644 index 00000000000..0f1b7ddd1f5 --- /dev/null +++ b/src/seis/longopt/shake_inc.h @@ -0,0 +1,48 @@ +/*-------------------------------------------------------------------- + * + * Copyright (c) 1991-2023 by the GMT Team (https://www.generic-mapping-tools.org/team.html) + * See LICENSE.TXT file for copying and redistribution conditions. + * + * 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.generic-mapping-tools.org + *--------------------------------------------------------------------*/ + +#ifndef SHAKE_INC_H +#define SHAKE_INC_H + +/* Translation table from long to short module options, directives and modifiers */ + +static struct GMT_KEYWORD_DICTIONARY module_kw[] = { /* Local options for this module */ + /* separator, short_option, long_option, + short_directives, long_directives, + short_modifiers, long_modifiers */ + { 0, 'G', "", + "", "", + "", "" }, + { 0, 'D', "", + "", "", + "", "" }, + { 0, 'L', "", + "", "", + "", "" }, + { 0, 'M', "", + "", "", + "", "" }, + { 0, 'C', "", + "", "", + "", "" }, + { 0, 'F', "", + "", "", + "", "" }, + { 0, '\0', "", "", "", "", ""} /* End of list marked with empty option and strings */ +}; + +#endif /* !SHAKE_INC_H */ diff --git a/src/seis/longopt/vs30_inc.h b/src/seis/longopt/vs30_inc.h new file mode 100644 index 00000000000..b2719ad33d5 --- /dev/null +++ b/src/seis/longopt/vs30_inc.h @@ -0,0 +1,39 @@ +/*-------------------------------------------------------------------- + * + * Copyright (c) 1991-2023 by the GMT Team (https://www.generic-mapping-tools.org/team.html) + * See LICENSE.TXT file for copying and redistribution conditions. + * + * 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.generic-mapping-tools.org + *--------------------------------------------------------------------*/ + +#ifndef VS30_INC_H +#define VS30_INC_H + +/* Translation table from long to short module options, directives and modifiers */ + +static struct GMT_KEYWORD_DICTIONARY module_kw[] = { /* Local options for this module */ + /* separator, short_option, long_option, + short_directives, long_directives, + short_modifiers, long_modifiers */ + { 0, 'G', "", + "", "", + "", "" }, + { 0, 'C', "", + "", "", + "", "" }, + { 0, 'W', "", + "", "", + "", "" }, + { 0, '\0', "", "", "", "", ""} /* End of list marked with empty option and strings */ +}; + +#endif /* !VS30_INC_H */ diff --git a/src/seis/shake.c b/src/seis/shake.c new file mode 100644 index 00000000000..76aff135dde --- /dev/null +++ b/src/seis/shake.c @@ -0,0 +1,506 @@ +/*-------------------------------------------------------------------- + * + * Copyright (c) 1991-2020 by the GMT Team (https://www.generic-mapping-tools.org/team.html) + * See LICENSE.TXT file for copying and redistribution conditions. + * + * 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.generic-mapping-tools.org + *--------------------------------------------------------------------*/ +/* + +/*-------------------------------------------------------------------- + * Compute Peak Ground Acceleration/Velocity and Intensity + * Author: J. Luis, after the FORTRAN version from J M Miranda + * Date: Original from 06-Apr-2010 for Mirone usage. Adapted to GMT in 2020 + * + */ + +#include "gmt_dev.h" + +#define THIS_MODULE_CLASSIC_NAME "shake" +#define THIS_MODULE_MODERN_NAME "shake" +#define THIS_MODULE_LIB "seis" +#define THIS_MODULE_PURPOSE "Compute Peak Ground Acceleration/Velocity and Intensity." +#define THIS_MODULE_KEYS " */ + bool active; + bool do_PGA, do_PGV, do_INT; + int n; /* Number of output grids specified via -G */ + char *file[3]; /* Only first is used for commandline but API may need many */ + } G; + struct SHAKE_L { /* -L[/] */ + bool active; + unsigned int mode; /* 0 = dist to nearest point, 1 = also get the point, 2 = instead get seg#, pt# */ + unsigned int sph; /* 0 = Flat Earth, 1 = spherical [Default], 2 = ellipsoidal */ + char *file; /* Name of file with lines */ + double line[4]; + char unit; + } L; + struct SHAKE_M { + bool active; + double rmw; + } M; +}; + +static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ + struct SHAKE_CTRL *C; + + C = gmt_M_memory (GMT, NULL, 1, struct SHAKE_CTRL); + + /* Initialize values whose defaults are not 0/false/NULL */ + C->F.imeca = 1; + C->L.mode = 1; /* Default returns distance to line */ + C->L.unit = 'k'; /* Default unit is km */ + C->L.sph = GMT_GREATCIRCLE; /* Default is great-circle distances */ + return (C); +} + +static void Free_Ctrl (struct GMT_CTRL *GMT, struct SHAKE_CTRL *C) { /* Deallocate control structure */ + if (!C) return; + gmt_M_free(GMT, C); +} + +static char set_unit_and_mode (char *arg, unsigned int *mode) { + unsigned int k = 0; + *mode = GMT_GREATCIRCLE; /* Default is great circle distances */ + switch (arg[0]) { + case '-': *mode = GMT_FLATEARTH; k = 1; break; + case '+': *mode = GMT_GEODESIC; k = 1; break; + } + return (arg[k]); +} + +static int usage (struct GMTAPI_CTRL *API, int level) { + const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE); + + if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); + GMT_Message (API, GMT_TIME_NONE, "usage: %s -G -L | -Dx0y0/x1/y1 -M [-Ca,v,i] [-F] [%s] [%s] [%s]\n", + name, GMT_Rgeoz_OPT, GMT_V_OPT, GMT_f_OPT); + + if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS); + + GMT_Message (API, GMT_TIME_NONE, "\t The grid with the Vs30 velocities.\n"); + GMT_Message (API, GMT_TIME_NONE, "\t-G Specify file name for output grid file(s).\n"); + GMT_Message (API, GMT_TIME_NONE, "\t If more than one component is set via -C then must contain %%s to format component code.\n"); + GMT_Message (API, GMT_TIME_NONE, "\t-D x0/y0/x1/y1 End points of the fault trace.\n"); + GMT_Message (API, GMT_TIME_NONE, "\t-L Alternatively provide a name of a file with the coordinates of the fault trace.\n"); + GMT_Message (API, GMT_TIME_NONE, "\t-M Select magnitude.\n"); + GMT_Message (API, GMT_TIME_NONE, "\n\tOPTIONS:\n"); + if (API->external) + GMT_Message (API, GMT_TIME_NONE, "\t-C List of comma-separated components to be written as grids. Choose from\n"); + else + GMT_Message (API, GMT_TIME_NONE, "\t-C List of comma-separated components to be written as grids (requires -G). Choose from\n"); + GMT_Message (API, GMT_TIME_NONE, "\t a, v, i. [Default is i(ntensity)].\n"); + GMT_Message (API, GMT_TIME_NONE, "\t-F Select focal mechanism type (e.g. -F1 or -F2 ...).\n"); + GMT_Message (API, GMT_TIME_NONE, "\t 1 unknown [Default].\n"); + GMT_Message (API, GMT_TIME_NONE, "\t 2 strike-slip.\n"); + GMT_Message (API, GMT_TIME_NONE, "\t 3 normal.\n"); + GMT_Message (API, GMT_TIME_NONE, "\t 4 thrust.\n\n"); + GMT_Option (API, "R,V"); + GMT_Option (API, "f,i,:"); + + return (GMT_MODULE_USAGE); +} + + +static int parse (struct GMT_CTRL *GMT, struct SHAKE_CTRL *Ctrl, struct GMT_Z_IO *io, struct GMT_OPTION *options) { + /* This parses the options provided to shake and sets parameters in Ctrl. + * Note Ctrl has already been initialized and non-zero default values set. + * 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, n_errors = 0, n_files = 0, pos = 0; + char txt_a[GMT_LEN256] = {""}, p[GMT_LEN16] = {""}; + struct GMT_OPTION *opt = NULL; + struct GMTAPI_CTRL *API = GMT->parent; + + gmt_M_memset (io, 1, struct GMT_Z_IO); + + for (opt = options; opt; opt = opt->next) { /* Process all the options given */ + + switch (opt->option) { + + case '<': /* Input file (only one is accepted) */ + if (n_files++ > 0) break; + if ((Ctrl->In.active = gmt_check_filearg (GMT, '<', opt->arg, GMT_IN, GMT_IS_GRID)) != 0) + Ctrl->In.file = strdup (opt->arg); + else + n_errors++; + break; + + /* Processes program-specific parameters */ + + case 'C': /* Requires -G and selects which components should be written as grids */ + Ctrl->C.active = true; + while ((gmt_strtok (opt->arg, ",", &pos, p)) && Ctrl->C.n_selected < 3) { + switch (p[0]) { + case 'a': Ctrl->C.selected[0] = Ctrl->G.do_PGA = true; break; + case 'v': Ctrl->C.selected[1] = Ctrl->G.do_PGV = true; break; + case 'i': Ctrl->C.selected[2] = Ctrl->G.do_INT = true; break; + default: + GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Unrecognized field argument %s in -C.!\n", p); + n_errors++; + break; + } + if (!Ctrl->C.selected[1]) Ctrl->C.no_PGV = true; + if (Ctrl->C.selected[2]) Ctrl->C.selected[1] = true; /* Intensity needs velocity too */ + Ctrl->C.n_selected++; + } + if (Ctrl->C.n_selected == 0) { + GMT_Report (GMT->parent, GMT_MSG_NORMAL, "-C requires comma-separated component arguments.\n"); + n_errors++; + } + break; + case 'D': /* x0/y0[x1/y1] */ + Ctrl->D.active = true; + Ctrl->D.line = strdup (opt->arg); + break; + case 'G': /* Output filename */ + if (!GMT->parent->external && Ctrl->G.n) { /* Command line interface */ + GMT_Report (GMT->parent, GMT_MSG_NORMAL, "-G can only be set once!\n"); + n_errors++; + } + else if ((Ctrl->G.active = gmt_check_filearg (GMT, 'G', opt->arg, GMT_OUT, GMT_IS_GRID)) != 0) + Ctrl->G.file[Ctrl->G.n++] = strdup (opt->arg); + else + n_errors++; + + if (!GMT->parent->external) { /* Copy the name into the 3 slots to simplify the grid writing algo */ + Ctrl->G.file[1] = strdup (opt->arg); + Ctrl->G.file[2] = strdup (opt->arg); + } + break; + case 'F': + Ctrl->F.active = true; + Ctrl->F.imeca = atoi (opt->arg); + if (Ctrl->F.imeca < 1 || Ctrl->F.imeca > 4) + GMT_Report (GMT->parent, GMT_MSG_NORMAL, "-F error. 'type' must be in [1 4]\n"); + break; + case 'L': /* -L[+u[+|-]] */ + Ctrl->L.active = true; + Ctrl->L.file = gmt_get_filename(API, opt->arg, "u"); + if (!gmt_check_filearg (GMT, 'L', Ctrl->L.file, GMT_IN, GMT_IS_DATASET)) { + GMT_Report (GMT->parent, GMT_MSG_NORMAL, "-L error. Must provide either an existing file name or line coordinates.\n"); + n_errors++; + } + break; + case 'M': + Ctrl->M.active = true; + Ctrl->M.rmw = atof(opt->arg); + break; + default: /* Report bad options */ + n_errors += gmt_default_error (GMT, opt->option); + break; + } + } + + if (!Ctrl->C.active) { + Ctrl->C.selected[1] = Ctrl->C.selected[2] = Ctrl->G.do_INT = Ctrl->C.no_PGV = true; /* The default */ + Ctrl->C.n_selected = 1; + } + + n_errors += gmt_M_check_condition (GMT, n_files != 1, "Syntax error: Must specify a single grid file\n"); + n_errors += gmt_M_check_condition (GMT, !Ctrl->L.active && !Ctrl->D.active, + "-L or -D option: Must provide a fault file name or coordinates.\n"); + n_errors += gmt_M_check_condition (GMT, !Ctrl->M.active, "-M option: Must specify magnitude.\n"); + + return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); +} + +#define bailout(code) {gmt_M_free_options (mode); return (code);} +#define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);} + +/* --------------------------------------------------------------------------------- */ +EXTERN_MSC int GMT_shake (void *V_API, int mode, void *args) { + uint64_t i, j, ij, k, kk, row, seg; + int error = 0, way, proj_type = 0; /* Geographic */ + char file[GMT_LEN512] = {""}, *code[3] = {"a", "v", "i"}; + + double wesn[4]; + double u = 1, ss = 0, rns = 0, rs = 0, dist = 0.0, xnear = 0.0, ynear = 0.0, xtmp, ytmp; + double rmh_pga, rmh_pgv, r_pga, r_pgv, rfd_pga, rfd_pgv, rfd_pga4nl, k1, k2, k3; + double flin_pga, flin_pgv, pga4nl, bnl_pga, bnl_pgv, fnl_pga, fnl_pgv, fs_pga, fs_pgv; + double rfm_pgv, rfm_pga, c_pga, d_pga, c_pgv, d_pgv, tmp, tmp1, tmp2, tmp3, rlon, rlat; + + struct GMT_DATATABLE *xyline = NULL; + struct GMT_DATASET *Lin = NULL; + struct GMT_GRID *G = NULL; + struct GMT_GRID *Grid[3] = {NULL, NULL, NULL}; + struct GMT_RECORD *Out = NULL; + struct GMT_Z_IO io; + struct GMT_OPTION *opt = NULL; + struct SHAKE_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 (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 ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */ + + /* Parse the command-line arguments */ + + if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, NULL, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */ + if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error); + Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */ + if ((error = parse (GMT, Ctrl, &io, options)) != 0) Return (error); + + /*---------------------------- This is the grd2xyz main code ----------------------------*/ + + gmt_M_memcpy (wesn, GMT->common.R.wesn, 4, double); /* Current -R setting, if any */ + + gmt_set_pad (GMT, 0); /* Change the default pad (ignored if input is a memory grid) */ + if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) /* Get header only */ + Return (API->error); + + if (gmt_M_is_subset (GMT, G->header, wesn)) /* If subset requested make sure wesn matches header spacing */ + gmt_M_err_fail (GMT, gmt_adjust_loose_wesn (GMT, wesn, G->header), ""); + + /* Read data */ + if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn, Ctrl->In.file, G) == NULL) + Return (API->error); + + gmt_M_memcpy (wesn, G->header->wesn, 4, double); /* Need the wesn further down */ + + gmt_set_geographic (GMT, GMT_OUT); + for (k = 0; k < 3; k++) { + if (!Ctrl->C.selected[k]) continue; + /* Create the empty grid and allocate space */ + if ((Grid[k] = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, wesn, + G->header->inc, GMT->common.R.registration, G->header->pad[0], NULL)) == NULL) + Return (API->error); + } + + way = gmt_M_is_geographic (GMT, GMT_IN) ? Ctrl->L.sph: 0; + proj_type = gmt_init_distaz (GMT, Ctrl->L.unit, way, GMT_MAP_DIST); + + if (Ctrl->D.active) { /* Gave -Dx0/y0/x1/y1 */ + int n; + uint64_t par[] = {1,1,2,2}; + double x0, y0, x1, y1; + Lin = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_LINE, GMT_CONTAINER_AND_DATA , par, NULL, NULL, 0, 0, NULL); + n = sscanf (Ctrl->D.line, "%lf/%lf/%lf/%lf", &x0, &y0, &x1, &y1); + if (n != 4 && n != 2) { + GMT_Report (API, GMT_MSG_ERROR, "Option -D: must provide either one x/y or two points x0/y0/x1/y1\n"); + bailout(GMT_PARSE_ERROR); + } + Lin->table[0]->segment[0]->data[GMT_X][0] = x0; Lin->table[0]->segment[0]->data[GMT_Y][0] = y0; + if (n == 4) { + Lin->table[0]->segment[0]->data[GMT_X][1] = x1; Lin->table[0]->segment[0]->data[GMT_Y][1] = y1; + } + else { /* Just repeat the x0/y0 point with a noisy shift */ + Lin->table[0]->segment[0]->data[GMT_X][1] = x0+1e-6; Lin->table[0]->segment[0]->data[GMT_Y][1] = y0+1e-6; + } + } + else { + if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_LINE, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) + Return (API->error); + + if ((Lin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_LINE, GMT_READ_NORMAL, NULL, Ctrl->L.file, NULL)) == NULL) + Return (API->error); + + if (Lin->n_columns < 2) { + GMT_Report (API, GMT_MSG_NORMAL, "Input data have %d column(s) but at least 2 are needed\n", (int)Lin->n_columns); + Return (GMT_DIM_TOO_SMALL); + } + gmt_set_segmentheader (GMT, GMT_OUT, false); /* Since processing of -L file might have turned it on [should be determined below] */ + } + + xyline = Lin->table[0]; /* Can only be one table since we read a single file */ + + if (proj_type == GMT_GEO2CART) { /* Must convert the line points first */ + for (seg = 0; seg < xyline->n_segments; seg++) { + for (row = 0; row < xyline->segment[seg]->n_rows; row++) { + gmt_geo_to_xy (GMT, xyline->segment[seg]->data[GMT_X][row], xyline->segment[seg]->data[GMT_Y][row], &xtmp, &ytmp); + xyline->segment[seg]->data[GMT_X][row] = xtmp; + xyline->segment[seg]->data[GMT_Y][row] = ytmp; + } + } + } + else if (gmt_M_is_geographic (GMT, GMT_IN) && proj_type == GMT_GEOGRAPHIC && !gmt_M_is_spherical (GMT)) { + /* Will need spherical trig so convert to geocentric latitudes if on an ellipsoid */ + for (seg = 0; seg < xyline->n_segments; seg++) { + for (row = 0; row < xyline->segment[seg]->n_rows; row++) { /* Convert to geocentric */ + xyline->segment[seg]->data[GMT_Y][row] = + gmt_lat_swap (GMT, xyline->segment[seg]->data[GMT_Y][row], GMT_LATSWAP_G2O); + } + } + } + + /* ------------------------------------------------------------- */ + if (Ctrl->F.imeca == 2) { + u = rns = rs = 0; ss = 1; + } + else if (Ctrl->F.imeca == 3) { + u = ss = rs = 0; rns = 1; + } + else if (Ctrl->F.imeca == 4) { + u = ss = rns = 0; rs = 1; + } + + /* - ------------------------------------------------------------ + ! calcula o termo de magnitude scaling para o pga e pgv + ! nao depende da distancia mas apenas da magnitude + ! ------------------------------------------------------------- */ + rmh_pga = 6.75; + rmh_pgv = 8.50; + rfm_pga = -0.53804*u - 0.50350*ss - 0.75472*rns - 0.50970*rs; + if (Ctrl->M.rmw <= rmh_pga) + rfm_pga += (0.28805*(Ctrl->M.rmw-rmh_pga) - 0.10164*(Ctrl->M.rmw-rmh_pga)*(Ctrl->M.rmw-rmh_pga)); + + rfm_pgv = 5.00121*u + 5.04727*ss + 4.63188*rns + 5.08210*rs; + if (Ctrl->M.rmw <= rmh_pgv) + rfm_pgv += (0.18322*(Ctrl->M.rmw-rmh_pgv) - 0.12736*(Ctrl->M.rmw-rmh_pgv)*(Ctrl->M.rmw-rmh_pgv)); + + /* - ------------------------------------------------------------ + ! iterate the positions matrix pga (i,j), pgv (i,j) e pint (i,j) + ! to calculate the scaling distance term + ! ------------------------------------------------------------- */ + /* The k2 & k3 coeff bellow is to simplify the computations of origianl equations such: + c_pga = +bnl_pga * (3.*0.40546510 - 1.0986123)/(1.0986123 * 1.0986123); */ + k2 = (3.*0.40546510 - 1.0986123) / (1.0986123 * 1.0986123); + k3 = (2.*0.40546510 - 1.0986123) / (1.0986123 * 1.0986123 * 1.0986123); + k1 = log(0.6); + + for (j = k = 0; j < G->header->n_rows; j++) { + rlat = G->header->wesn[YHI] - j * G->header->inc[GMT_Y]; + for (i = 0; i < G->header->n_columns; i++, k++) { + rlon = G->header->wesn[XLO] + i * G->header->inc[GMT_X]; + ij = gmt_M_ijp(G->header, j, i); + //ij = k; + + (void)gmt_near_lines (GMT, rlon, rlat, xyline, Ctrl->L.mode, &dist, &xnear, &ynear); + + r_pga = sqrt (dist*dist + 1.35*1.35); + r_pgv = sqrt (dist*dist + 2.54*2.54); + rfd_pga = (-0.66050+0.11970*(Ctrl->M.rmw-4.5))*log(r_pga/1.0) - 0.01151*(r_pga-1.0); + rfd_pgv = (-0.87370+0.10060*(Ctrl->M.rmw-4.5))*log(r_pgv/1.0) - 0.00334*(r_pgv-1.0); + rfd_pga4nl = (-0.66050+0.11970*(Ctrl->M.rmw-4.5))*log(r_pga/5.0) - 0.01151*(r_pga-5.0); + + /* ------------------------------------------------------------- + ! calcula o efeito de sitio + ! ------------------------------------------------------------- */ + tmp = log(G->data[ij] / 760); + flin_pga = -0.360 * tmp; + flin_pgv = -0.600 * tmp; + + pga4nl = exp(rfm_pga + rfd_pga4nl); + + if (G->data[ij] <= 180.) { + bnl_pga = -0.640; + bnl_pgv = -0.600; + } + else if (G->data[ij] <= 300 && G->data[ij] > 180.) { + tmp = log(G->data[ij]/300.) / log(180./300.); + bnl_pga = (-0.64+0.14) * tmp - 0.14; + bnl_pgv = (-0.60+0.06) * tmp - 0.14; + } + else if (G->data[ij] <= 760 && G->data[ij] > 300.) { + tmp = log(G->data[ij]/760.) / log(300./760.); + bnl_pga = -0.14 * tmp; + bnl_pgv = -0.06 * tmp; + } + else { + bnl_pga = bnl_pgv = 0.; + } + + /* ---- verificar se as condicoes s�o sempre em pga !!!!!!!!!!!!!!!!!!!!! */ + c_pga = +bnl_pga * k2; d_pga = -bnl_pga * k3; + c_pgv = +bnl_pgv * k2; d_pgv = -bnl_pgv * k3; + + if (pga4nl <= 0.03) { + fnl_pga = bnl_pga * k1; + fnl_pgv = bnl_pgv * k1; + } + else if (pga4nl > 0.03 && pga4nl <= 0.09) { + tmp3 = log(pga4nl/0.03); tmp1 = tmp3 * tmp3; tmp2 = tmp1 * tmp3; + fnl_pga = bnl_pga * k1 + c_pga * tmp1 + d_pga * tmp2; + fnl_pgv = bnl_pgv * k1 + c_pgv * tmp1 + d_pgv * tmp2; + } + else if (pga4nl > 0.09) { + tmp = log (pga4nl/0.1); + fnl_pga = bnl_pga * tmp; + fnl_pgv = bnl_pgv * tmp; + } + + fs_pga = flin_pga + fnl_pga; + fs_pgv = flin_pgv + fnl_pgv; + + if (Ctrl->G.do_PGA) + Grid[0]->data[ij] = (float)(exp (rfm_pga + rfd_pga + fs_pga) * 980.); + + if (Ctrl->G.do_PGV || Ctrl->G.do_INT) + Grid[1]->data[ij] = (float)exp (rfm_pgv + rfd_pgv + fs_pgv); + + if (Ctrl->G.do_INT) { + tmp = log10(Grid[1]->data[ij]); + Grid[2]->data[ij] = (float)(4.398 + 1.916 * tmp + 0.280 * tmp*tmp); + } + } + } + + /* Now write the one to three grids */ + if (Ctrl->C.no_PGV) Ctrl->C.selected[1] = false; /* Don't save PGV, its usage was only temporary */ + for (k = kk = 0; k < 3; k++) { + if (!Ctrl->C.selected[k]) continue; + if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Grid[k])) + Return (API->error); + + if (!API->external) kk = k; /* On command line we pick item k from an array of 3 items */ + if (strstr (Ctrl->G.file[kk], "%s")) + sprintf (file, Ctrl->G.file[kk], code[k]); + else + strncpy (file, Ctrl->G.file[kk], GMT_LEN512-1); + + if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, file, Grid[k]) != GMT_NOERROR) { + Return (API->error); + } + kk++; /* For the external interface we take them in the order given as there is not an array of 3 */ + } + + if (Ctrl->L.active) + GMT_End_IO (API, GMT_IN, 0); /* Disables further data input */ + + if (GMT_Destroy_Data (GMT->parent, &G) != GMT_NOERROR) + GMT_Report (API, GMT_MSG_NORMAL, "Failed to free G\n"); + + Return (GMT_NOERROR); +} diff --git a/src/seis/vs30.c b/src/seis/vs30.c new file mode 100644 index 00000000000..f0eaff114bb --- /dev/null +++ b/src/seis/vs30.c @@ -0,0 +1,453 @@ +/* + * grad2vs30: convert topographic slope to Vs30 + * + * Takes three input files, each formatted GMT grd files + * and point for point co-registered with each other; + * gradient_file contains the topographic slope expressed as + * a unitless ratio (e.g., meters per meter), landmask_file + * is 0 for water and 1 (one) for land; craton_file is a + * weight ranging from 1 (one) on stable shields (craton) and + * 0 in active tectonic regions -- values in between will + * be computed as the weighted average of the craton and + * tectonic models. + * The optional numerical argument "water" is the value + * that water-covered areas will be set to; the default is 600. + * The output file name is specified with "output_file" (required). + * + * LICENSE + * (https://github.com/usgs/earthquake-global_vs30/blob/master/LICENSE.md) + * + * Unless otherwise noted, This project is in the public domain in the United States because it + * contains materials that originally came from the United States Geological Survey, an agency + * of the United States Department of Interior. For more information, see the official USGS + * copyright policy at https://www2.usgs.gov/visual-id/credit_usgs.html#copyright + * + * Additionally, we waive copyright and related rights in the work worldwide through the CC0 1.0 + * Universal public domain dedication. + * + * The getpar package (src/getpar.c, src/libget.h, src/getpar.3) is included with the express + * permission of the author, Robert W. Clayton. + * + * J. Luis + * re-writen after https://github.com/usgs/earthquake-global_vs30/blob/master/src/grad2vs30.c + * Version: 6 API + */ + +//shake vs30.grd -Glixo.grd -Lline.dat+uk -Ci -Vl -Rvs30.grd + +#include "gmt_dev.h" + +#define THIS_MODULE_CLASSIC_NAME "vs30" +#define THIS_MODULE_MODERN_NAME "vs30" +#define THIS_MODULE_LIB "seis" +#define THIS_MODULE_PURPOSE "Compute VS30" +#define THIS_MODULE_KEYS " | fname[+g] */ + bool active; + bool is_grid; /* To signal that a craton grid was provided */ + char *file; + float val; /* a [0 1] val where 0 (craton) and 1 (active) */ + } C; + struct VS30_G { + bool active; + char *file; + } G; + struct VS30_W { + bool active; + float water; + } W; +}; + +static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ + struct VS30_CTRL *C; + + C = gmt_M_memory (GMT, NULL, 1, struct VS30_CTRL); + + /* Initialize values whose defaults are not 0/false/NULL */ + C->W.water = 600; /* */ + return (C); +} + +static void Free_Ctrl (struct GMT_CTRL *GMT, struct VS30_CTRL *C) { /* Deallocate control structure */ + if (!C) return; + gmt_M_free(GMT, C); +} + +/* + * Slope to Vs30 uses Wald & Allen (2007) for craton, and Allen & Wald (2009) for active tectonic. + * + * Split up the tables just in case future work makes them have different numbers of rows + * Columns are: vs30_min vs30_max slope_min slope_max + */ +const unsigned int rows_active = 6; +double active_table[6][4] = + {{180, 240, 3.0e-4, 3.5e-3}, + {240, 300, 3.5e-3, 0.01}, + {300, 360, 0.01, 0.018}, + {360, 490, 0.018, 0.05}, + {490, 620, 0.05, 0.10}, + {620, 760, 0.10, 0.14}}; + +const unsigned int rows_craton = 6; +double craton_table[6][4] = + {{180, 240, 2.0e-5, 2.0e-3}, + {240, 300, 2.0e-3, 4.0e-3}, + {300, 360, 4.0e-3, 7.2e-3}, + {360, 490, 7.2e-3, 0.013}, + {490, 620, 0.013, 0.018}, + {620, 760, 0.018, 0.025}}; + + +static int usage (struct GMTAPI_CTRL *API, int level) { + const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE); + + if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); + GMT_Message (API, GMT_TIME_NONE, "usage: %s -G [-C|fname[+g]] [-W] [%s] [%s]\n", + name, GMT_Rgeoz_OPT, GMT_V_OPT); + + if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS); + + GMT_Message (API, GMT_TIME_NONE, "\t The input grid name of the topography.\n"); + GMT_Message (API, GMT_TIME_NONE, "\t-G File name for output grid with the Vs30 velocities.\n"); + GMT_Message (API, GMT_TIME_NONE, "\t-C Argument can be one of three:\n"); + GMT_Message (API, GMT_TIME_NONE, "\t - A value between 0 and 1, where 0 means a stable Craton and 1 an Active region.\n"); + GMT_Message (API, GMT_TIME_NONE, "\t - The name of a multi-segment file with the 'cratons' polygons. In this case the\n"); + GMT_Message (API, GMT_TIME_NONE, "\t polygons will be feed to grdmask to compute a cratons/active tectonic mask.\n"); + GMT_Message (API, GMT_TIME_NONE, "\t - The name of a grid with the cratons/active tectonic regions.\n"); + GMT_Message (API, GMT_TIME_NONE, "\n\tOPTIONS:\n"); + GMT_Message (API, GMT_TIME_NONE, "\t-W water_vel sets the Vs30 value used in areas designated as water in the\n"); + GMT_Message (API, GMT_TIME_NONE, "\t landmask [default=600]\n\n"); + GMT_Option (API, "R,V"); + GMT_Option (API, "bi,i,r,:"); + + return (GMT_MODULE_USAGE); +} + +static int parse (struct GMT_CTRL *GMT, struct VS30_CTRL *Ctrl, struct GMT_Z_IO *io, struct GMT_OPTION *options) { + /* This parses the options provided to shake and sets parameters in Ctrl. + * Note Ctrl has already been initialized and non-zero default values set. + * 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, pos = 0; + char txt_a[GMT_LEN256] = {""}, p[GMT_LEN16] = {""}, *pch; + struct GMT_OPTION *opt = NULL; + struct GMTAPI_CTRL *API = GMT->parent; + + gmt_M_memset (io, 1, struct GMT_Z_IO); + + for (opt = options; opt; opt = opt->next) { /* Process all the options given */ + + switch (opt->option) { + + case '<': /* Input file (only one is accepted) */ + if (n_files++ > 0) break; + if ((Ctrl->In.active = gmt_check_filearg (GMT, '<', opt->arg, GMT_IN, GMT_IS_GRID)) != 0) + Ctrl->In.file = strdup (opt->arg); + else + n_errors++; + break; + + /* Processes program-specific parameters */ + + case 'C': /* Polygon file with the Cratons limits */ + Ctrl->C.active = true; + Ctrl->C.file = gmt_get_filename (API, opt->arg, ""); + if ((pch = strstr(opt->arg, "+g")) != NULL) + Ctrl->C.is_grid = true; /* Useless if Ctrl->C.file = NULL is set below */ + + if (!gmt_check_filearg (GMT, 'C', Ctrl->C.file, GMT_IN, GMT_IS_DATASET)) { + Ctrl->C.file = NULL; + Ctrl->C.val = atof (opt->arg); + if (Ctrl->C.val < 0 || Ctrl->C.val > 1) { + GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Errorn in -C. Must provide either a file name or a value in the [0 1] interval.\n", p); + n_errors++; + } + } + break; + case 'G': /* Output file */ + if ((Ctrl->G.active = gmt_check_filearg (GMT, 'G', opt->arg, GMT_OUT, GMT_IS_GRID))) + Ctrl->G.file = strdup (opt->arg); + else + n_errors++; + break; + case 'W': + Ctrl->W.water = (float)atof (opt->arg); + break; + default: /* Report bad options */ + n_errors += gmt_default_error (GMT, opt->option); + break; + } + } + + n_errors += gmt_M_check_condition (GMT, !Ctrl->G.file, "Syntax error -G option: Must specify output grid file\n"); + n_errors += gmt_M_check_condition (GMT, !Ctrl->C.active, "Syntax error -C option: Must specify a value or a file name.\n"); + + return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); +} + +static int check_grid_compat (struct GMTAPI_CTRL *API, struct GMT_GRID *A, struct GMT_GRID *B) { + /* Check that grids A & B are compatible. If they are, return 0. Otherwise: + return 1 if increments differ more than 0.2 % + return 1 if any of the corners differ more than 1/5 of the grid spacing + return 3 if registrations are not equal + */ + if (fabs((A->header->inc[GMT_X] - B->header->inc[GMT_X]) / A->header->inc[GMT_X]) > 0.002 || + fabs((A->header->inc[GMT_Y] - B->header->inc[GMT_Y]) / A->header->inc[GMT_Y]) > 0.002) + return 1; + + if (fabs((A->header->wesn[XLO] - B->header->wesn[XLO]) / A->header->inc[GMT_X]) > 0.2 || + fabs((A->header->wesn[XHI] - B->header->wesn[XHI]) / A->header->inc[GMT_X]) > 0.2 || + fabs((A->header->wesn[YLO] - B->header->wesn[YLO]) / A->header->inc[GMT_Y]) > 0.2 || + fabs((A->header->wesn[YHI] - B->header->wesn[YHI]) / A->header->inc[GMT_Y]) > 0.2) + return 2; + + if (A->header->registration != B->header->registration) + return 3; + + return 0; +} + +/* + * Function interpVs30 interpolates (or extrapolates) vs30 from a row of one of the tables + * (tables have already been log()'ed so the exp() returns vs30 in linear units) + * (I was going to pre-compute the differences (tt[1]-tt[0] and * tt[3]-tt[2]) and make this + * function a #define for speed, but the execution time is utterly dwarfed by the read/write + * times so there was no point.) + */ +inline double interpVs30(double *tt, double lg) { + return exp(tt[0] + (tt[1] - tt[0]) * (lg - tt[2]) / (tt[3] - tt[2])); +} + +#define bailout(code) {gmt_M_free_options (mode); return (code);} +#define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);} + +/* --------------------------------------------------------------------------------- */ +EXTERN_MSC int GMT_vs30 (void *V_API, int mode, void *args) { + unsigned int row, col, j, nr, k; + uint64_t ij; + int error = 0; + char cmd[GMT_LEN256] = {""}, data_grd[GMT_LEN16] = {""}; + char crat_grd[GMT_LEN16] = {""}, mask_grd[GMT_LEN16] = {""}, grad_grd[GMT_LEN16] = {""}; + float crat, lg; + double (*table)[4], tvs[2], vv, wesn[4]; + + struct GMT_GRID *G = NULL, *Ggrad = NULL, *Gcrat = NULL, *Gland = NULL, *Gout = NULL; + struct GMT_Z_IO io; + struct GMT_OPTION *opt = NULL; + struct VS30_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 (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 ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */ + + /* Parse the command-line arguments */ + + if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, NULL, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */ + if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error); + Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */ + if ((error = parse (GMT, Ctrl, &io, options)) != 0) Return (error); + + /*---------------------------- This is the grd2xyz main code ----------------------------*/ + + /* Initialize the input objects and open the files */ + GMT_Report (API, GMT_MSG_VERBOSE, "Reading input files...\n"); + + gmt_M_memcpy (wesn, GMT->common.R.wesn, 4, double); /* Current -R setting, if any */ + + if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) /* Get header only */ + Return (API->error); + + if (gmt_M_is_subset (GMT, G->header, wesn)) /* If subset requested make sure wesn matches header spacing */ + gmt_M_err_fail (GMT, gmt_adjust_loose_wesn (GMT, wesn, G->header), ""); + + /* Read data */ + if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn, Ctrl->In.file, G) == NULL) + Return (API->error); + + gmt_M_memcpy (wesn, G->header->wesn, 4, double); /* Need the wesn further down */ + + /* ------------------------------------------------------------------------------- */ + GMT_Report (API, GMT_MSG_VERBOSE, "Compute landmask grid\n"); + /* Create a virtual file to hold the landmask grid */ + if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT, NULL, mask_grd)) + Return (API->error); + + /* Prepare the grdlandmask arguments */ + sprintf (cmd, "-G%s -I%f/%f -Df -R%.16g/%.16g/%.16g/%.16g --GMT_HISTORY=false ", + mask_grd, G->header->inc[GMT_X], G->header->inc[GMT_Y], wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI]); + + GMT_Report (API, GMT_MSG_LONG_VERBOSE, "Calling grdlandmask with args %s\n", cmd); + if (GMT_Call_Module (API, "grdlandmask", GMT_MODULE_CMD, cmd)) + Return (API->error); + + /* Obtain the data from the virtual file */ + if ((Gland = GMT_Read_VirtualFile (API, mask_grd)) == NULL) + Return (API->error); + /* ------------------------------------------------------------------------------- */ + + /* ------------------------------------------------------------------------------- */ + GMT_Report (API, GMT_MSG_VERBOSE, "Derive slope grid from input data grid\n"); + /* Create a virtual file to hold the gradient grid */ + if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT, NULL, grad_grd)) + Return (API->error); + + /* Prepare the grdgradient arguments */ + sprintf (cmd, "-S%s -n+bg -fg -D -R%.16g/%.16g/%.16g/%.16g --GMT_HISTORY=false ", + grad_grd, wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI]); + + strcat (cmd, Ctrl->In.file); + GMT_Report (API, GMT_MSG_LONG_VERBOSE, "Calling grdgradient with args %s\n", cmd); + if (GMT_Call_Module (API, "grdgradient", GMT_MODULE_CMD, cmd)) + Return (API->error); + + /* Obtain the data from the virtual file */ + if ((Ggrad = GMT_Read_VirtualFile (API, grad_grd)) == NULL) + Return (API->error); + /* ------------------------------------------------------------------------------- */ + + if (Ctrl->C.file && !Ctrl->C.is_grid) { /* Read the Cratons multi-seg polygons and compute a mask */ + GMT_Report (API, GMT_MSG_VERBOSE, "Derive cratons grid from input file %s\n", Ctrl->C.file); + /* Create a virtual file to hold the gradient grid */ + if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT, NULL, crat_grd)) + Return (API->error); + + /* Prepare the grdmask arguments */ + sprintf (cmd, "-G%s -I%f/%f -R%.16g/%.16g/%.16g/%.16g -N1/0/0 --GMT_HISTORY=false ", + crat_grd, G->header->inc[GMT_X], G->header->inc[GMT_Y], wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI]); + + strcat (cmd, Ctrl->C.file); + GMT_Report (API, GMT_MSG_LONG_VERBOSE, "Calling grdmask with args %s\n", cmd); + if (GMT_Call_Module (API, "grdmask", GMT_MODULE_CMD, cmd)) + Return (API->error); + + /* Obtain the data from the virtual file */ + if ((Gcrat = GMT_Read_VirtualFile (API, crat_grd)) == NULL) + Return (API->error); + } + else if (Ctrl->C.file && Ctrl->C.is_grid) { /* Read a cratons grid */ + if ((Gcrat = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, + NULL, Ctrl->C.file, NULL)) == NULL) { /* Get header only */ + Return (API->error); + } + + if (check_grid_compat (API, Gcrat, G)) { + GMT_Report (API, GMT_MSG_NORMAL, "Cratons and topography grids are not compatible.\n", cmd); + Return (API->error); + } + + /* Read data */ + if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn, Ctrl->C.file, Gcrat) == NULL) { + Return (API->error); + } + } + + /* The output file has the same dimensions as Ggrad * so write the header, + * prep the output object, then open the output for writing. + */ + if ((Gout = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, wesn, G->header->inc, + GMT->common.R.registration, 2, NULL)) == NULL) + Return (API->error); + + /* * We're doing log-log interpolation, so log() everything in the tables first, for efficiency */ + for (row = 0; row < rows_active; row++) + for (col = 0; col < 4; col++) + active_table[row][col] = log(active_table[row][col]); + + for (row = 0; row < rows_craton; row++) + for (col = 0; col < 4; col++) + craton_table[row][col] = log(craton_table[row][col]); + + for (row = 0; row < Gout->header->n_rows; row++) { + for (col = 0; col < Gout->header->n_columns; col++) { + ij = gmt_M_ijp(G->header, row, col); + if (Gland->data[ij] == 0) { /* Set areas covered by water to the water value */ + Gout->data[ij] = Ctrl->W.water; + continue; + } + + crat = Ctrl->C.file != NULL ? Gcrat->data[ij] : Ctrl->C.val; + + /* This is the slope value to be converted to vs30 */ + lg = log(Ggrad->data[ij]); + + for (k = 0; k < 2; k++) { /* Get the Vs30 for both craton and active */ + /* k == 0 => craton, k == 1 => active */ + if (k == 0) { + table = craton_table; + nr = rows_craton; + } + else { + table = active_table; + nr = rows_active; + } + + /* + * Handle slopes lower than the minimum in the table by extrapolation capped by the minimum vs30 + * (this isn't necessary when the table contains the minimum vs30, but it's cheap insurance if + * we change the table or minimum -- ditto for the max, below) + */ + if (lg <= table[0][2]) { + vv = interpVs30(table[0], lg); + if (vv < vs30_min) vv = vs30_min; + tvs[k] = vv; + continue; + } + /* + * Handle slopes greater than the maximum in the table by extrapolation capped by the maximum vs30 + */ + if (lg >= table[nr-1][3]) { + vv = interpVs30(table[nr-1],lg); + if (vv > vs30_max) vv = vs30_max; + tvs[k] = vv; + continue; + } + /* All other slopes should be handled within the tables */ + for (j = 0; j < nr; j++) { + if (lg <= table[j][3]) { + tvs[k] = interpVs30(table[j], lg); + break; + } + } + } + /* Do a weighted average of craton and active vs30 */ + Gout->data[ij] = (float)(crat * tvs[0] + (1.0 - crat) * tvs[1]); + } + + //if (row % 100 == 0) + //GMT_Report (API, GMT_MSG_VERBOSE, "Done with %ld of %ld elements\r", row * Gout->header->n_columns, Gout->header->nm); + } + + GMT_Report (API, GMT_MSG_VERBOSE, "Writing output file...\n"); + if (GMT_Write_Data(API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Gout) != 0) + Return (API->error); + + GMT_Report (API, GMT_MSG_VERBOSE, "Done.\n"); + + if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) /* Disables further data output */ + Return (API->error); + + Return (GMT_NOERROR); +}