/*-------------------------------------------------------------------- * $Id: grdgradient.c,v 1.4.4.2 2002/02/27 17:41:10 pwessel Exp $ * * Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith * See COPYING 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 General Public License as published by * the Free Software Foundation; version 2 of the License. * * 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 General Public License for more details. * * Contact info: gmt.soest.hawaii.edu *--------------------------------------------------------------------*/ /* * grdgradient.c * read a grdfile and compute gradient in azim direction: * * azim = azimuth clockwise from north in degrees. * * gradient = -[(dz/dx)sin(azim) + (dz/dy)cos(azim)]. * * the expression in [] is the correct gradient. We take * -[] in order that data which goes DOWNHILL in the * azim direction will give a positive value; this is * for image shading purposes. * * * Author: W.H.F. Smith * Date: 13 Feb 1991 * Upgraded to v2.0 15-May-1991 Paul Wessel * * Modified: 1 Mar 94 by WHFS to make -M scale change with j latitude * 1 Mar 96 by PW to find gradient direction and magnitude (-S and -D) * 13 Mar 96 by WHFS to add exp trans and user-supplied sigma to -N * option, and add optional second azimuth to -A option. * 11 Sep 97 by PW now may pass average gradient along with sigma in -N * 22 Apr 98 by WHFS to add boundary conditions, switch sense of -S and * -D, and switch -Da to -Dc, for consistency of args. * Version: 3.4.1 */ #include "gmt.h" float *data, *slp; main (int argc, char **argv) { int i, j, ij, k, n, nm, nm2, mx, my; int p[4], n_used = 0; BOOLEAN error = FALSE, map_units = FALSE, normalize = FALSE, atan_trans = FALSE, bad, do_direct_deriv = FALSE; BOOLEAN find_directions = FALSE, do_cartesian = FALSE, do_orientations = FALSE, save_slopes = FALSE, add_ninety = FALSE; BOOLEAN sigma_set = FALSE, offset_set = FALSE, exp_trans = FALSE, two_azims = FALSE; double dx_grid, dy_grid, x_factor, y_factor, dzdx, dzdy, ave_gradient, norm_val = 1.0, sigma = 0.0; double azim, denom, max_gradient = 0.0, min_gradient = 0.0, rpi, m_pr_degree, lat, azim2; double x_factor2, y_factor2, dzdx2, dzdy2, dzds1, dzds2, offset; char *infile = CNULL, *outfile = CNULL, *slopefile = CNULL, format[BUFSIZ]; struct GRD_HEADER header; struct GMT_EDGEINFO edgeinfo; argc = GMT_begin (argc, argv); GMT_boundcond_init (&edgeinfo); for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch (argv[i][1]) { /* Common parameters */ case 'V': case '\0': error += GMT_get_common_args (argv[i], 0, 0, 0, 0); break; /* Supplemental parameters */ case 'A': do_direct_deriv = TRUE; j = sscanf(&argv[i][2], "%lf/%lf", &azim, &azim2); two_azims = (j == 2); break; case 'D': find_directions = TRUE; j = 2; while (argv[i][j]) { switch (argv[i][j]) { case 'C': case 'c': do_cartesian = TRUE; break; case 'O': case 'o': do_orientations = TRUE; break; case 'N': case 'n': add_ninety = TRUE; break; default: fprintf (stderr, "%s: GMT SYNTAX ERROR -S option: Unrecognized modifier\n", GMT_program); error++; break; } j++; } break; case 'G': outfile = &argv[i][2]; break; case 'L': error += GMT_boundcond_parse (&edgeinfo, &argv[i][2]); break; case 'M': map_units = TRUE; break; case 'N': normalize = TRUE; j = 2; if (argv[i][j]) { if (argv[i][j] == 't' || argv[i][j] == 'T') { atan_trans = TRUE; j++; } else if (argv[i][j] == 'e' || argv[i][j] == 'E') { exp_trans = TRUE; j++; } j = sscanf(&argv[i][j], "%lf/%lf/%lf", &norm_val, &sigma, &offset); if (j >= 2) sigma_set = TRUE; if (j == 3) offset_set = TRUE; } break; case 'S': save_slopes = TRUE; slopefile = &argv[i][2]; break; default: error = TRUE; GMT_default_error (argv[i][1]); break; } } else infile = argv[i]; } if (argc == 1 || GMT_quick) { fprintf (stderr,"grdgradient %s - Compute directional gradients from grdfiles\n\n", GMT_VERSION); fprintf (stderr, "usage: grdgradient -G [-A[/]] [-D[a][o][n]] [-L]\n"); fprintf (stderr, "[-M] [-N[t_or_e][[/[/]]]] [-S] [-V]\n\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf (stderr,"\t is name of input grdfile\n"); fprintf (stderr,"\n\tOPTIONS:\n"); fprintf (stderr, "\t-A sets azimuth (0-360 CW from North (+y)) for directional derivatives\n"); fprintf (stderr, "\t -A/ will compute two directions and save the one larger in magnitude.\n"); fprintf (stderr, "\t-D finds the direction of grad z.\n"); fprintf (stderr, "\t Append c to get cartesian angle (0-360 CCW from East (+x)) [Default: azimuth]\n"); fprintf (stderr, "\t Append o to get bidirectional orientations [0-180] rather than directions [0-360]\n"); fprintf (stderr, "\t Append n to add 90 degrees to the values from c or o\n"); fprintf (stderr, "\t-G output file for results from -A or -D\n"); fprintf (stderr, "\t-L sets boundary conditions. can be either\n"); fprintf (stderr, "\t g for geographic boundary conditions\n"); fprintf (stderr, "\t or one or both of\n"); fprintf (stderr, "\t x for periodic boundary conditions on x\n"); fprintf (stderr, "\t y for periodic boundary conditions on y\n"); fprintf (stderr, "\t [Default: Natural conditions]\n"); fprintf (stderr, "\t-M to use map units. In this case, dx,dy of grdfile\n"); fprintf (stderr, "\t will be converted from degrees lon,lat into meters.\n"); fprintf (stderr, "\t Default computes gradient in units of data/grid_distance.\n"); fprintf (stderr, "\t-N will normalize gradients so that max |grad| = [1.0]\n"); fprintf (stderr, "\t -Nt will make atan transform, then scale to [1.0]\n"); fprintf (stderr, "\t -Ne will make exp transform, then scale to [1.0]\n"); fprintf (stderr, "\t -Nt/[/] or -Ne/[/] sets sigma\n"); fprintf (stderr, "\t (and offset) for transform. [sigma, offset estimated from data]\n"); fprintf (stderr, "\t-S output file for |grad z|; requires -D\n"); GMT_explain_option ('V'); exit (EXIT_FAILURE); } if (!(do_direct_deriv || find_directions)) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify -A or -D\n", GMT_program); error++; } if (save_slopes && !slopefile) { fprintf (stderr, "%s: GMT SYNTAX ERROR -S option: Must specify output file\n", GMT_program); error++; } if (!outfile) { fprintf (stderr, "%s: GMT SYNTAX ERROR -G option: Must specify output file\n", GMT_program); error++; } if (!infile) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify input file\n", GMT_program); error++; } if (do_direct_deriv && (azim < 0.0 || azim >= 360.0)) { fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: Use 0-360 degree range\n", GMT_program); error++; } if (two_azims && (azim2 < 0.0 || azim2 >= 360.0)) { fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: Use 0-360 degree range\n", GMT_program); error++; } if (norm_val <= 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -N option: Normalization amplitude must be > 0\n", GMT_program); error++; } if (sigma_set && (sigma <= 0.0) ) { fprintf (stderr, "%s: GMT SYNTAX ERROR -N option: Sigma must be > 0\n", GMT_program); error++; } if (error) exit (EXIT_FAILURE); GMT_put_history (argc, argv); /* Update .gmtcommands */ if (GMT_read_grd_info (infile, &header)) { fprintf (stderr, "%s: Error opening file %s (not found?)\n", infile, GMT_program); exit (EXIT_FAILURE); } GMT_boundcond_param_prep (&header, &edgeinfo); GMT_grd_init (&header, argc, argv, TRUE); nm = header.nx * header.ny; mx = header.nx + 4; my = header.ny + 4; nm2 = mx * my; data = (float *) GMT_memory (VNULL, (size_t)nm2, sizeof (float), GMT_program); GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 2; if (save_slopes) slp = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program); GMT_read_grd (infile, &header, data, header.x_min, header.x_max, header.y_min, header.y_max, GMT_pad, FALSE); /* set boundary conditions: */ GMT_boundcond_set (&header, &edgeinfo, GMT_pad, data); if (map_units) { m_pr_degree = 2.0 * M_PI * gmtdefs.ellipse[gmtdefs.ellipsoid].eq_radius / 360.0; dx_grid = m_pr_degree * header.x_inc * cos (M_PI * (header.y_max + header.y_min) / 360.0); dy_grid = m_pr_degree * header.y_inc; } else { dx_grid = header.x_inc; dy_grid = header.y_inc; } x_factor = -1.0 / (2.0 * dx_grid); y_factor = -1.0 / (2.0 * dy_grid); if (do_direct_deriv) { if (two_azims) { azim2 *= (M_PI / 180.0); x_factor2 = x_factor * sin(azim2); y_factor2 = y_factor * cos(azim2); } azim *= (M_PI / 180.0); x_factor *= sin(azim); y_factor *= cos(azim); } p[0] = 1; p[1] = -1; p[2] = mx; p[3] = -mx; min_gradient = DBL_MAX; max_gradient = -DBL_MAX; ave_gradient = 0.0; for (j = k = 0; j < header.ny; j++) { if (map_units) { lat = (header.node_offset) ? -header.y_inc * (j + 0.5) : -header.y_inc * j; lat += header.y_max; dx_grid = m_pr_degree * header.x_inc * cos (D2R * lat); x_factor = -1.0 / (2.0 * dx_grid); if (do_direct_deriv) { if (two_azims) { x_factor2 = x_factor * sin(azim2); } x_factor *= sin(azim); } } for (i = 0; i < header.nx; i++, k++) { ij = (j + 2) * mx + i + 2; for (n = 0, bad = FALSE; !bad && n < 4; n++) if (GMT_is_fnan (data[ij+p[n]])) bad = TRUE; if (bad) { /* One of corners = NaN, skip */ data[k] = GMT_f_NaN; continue; } dzdx = (data[ij+1] - data[ij-1]) * x_factor; dzdy = (data[ij-mx] - data[ij+mx]) * y_factor; if (two_azims) { dzdx2 = (data[ij+1] - data[ij-1]) * x_factor2; dzdy2 = (data[ij-mx] - data[ij+mx]) * y_factor2; } /* Write output to unused NW corner */ if (do_direct_deriv) { /* Directional derivatives */ if (two_azims) { dzds1 = dzdx + dzdy; dzds2 = dzdx2 + dzdy2; data[k] = (float)((fabs(dzds1) > fabs(dzds2)) ? dzds1 : dzds2); } else { data[k] = (float)(dzdx + dzdy); } ave_gradient += data[k]; min_gradient = MIN (min_gradient, data[k]); max_gradient = MAX (max_gradient, data[k]); } else { azim = (do_cartesian) ? atan2 (-dzdy, -dzdx) * R2D : 90.0 - atan2 (-dzdy, -dzdx) * R2D; if (add_ninety) azim += 90.0; if (azim < 0.0) azim += 360.0; if (azim >= 360.0) azim -= 360.0; if (do_orientations && azim >= 180) azim -= 180.0; data[k] = (float)azim; if (save_slopes) slp[k] = (float)hypot (dzdx, dzdy); } n_used++; } } if (offset_set) ave_gradient = offset; else ave_gradient /= n_used; if (do_direct_deriv) { /* Report some statistics */ if (normalize) { if (atan_trans) { if (sigma_set) { denom = 1.0 / sigma; } else { denom = 0.0; for (k = 0; k < nm; k++) if (!GMT_is_fnan (data[k])) denom += pow(data[k] - ave_gradient, 2.0); denom = sqrt( (n_used - 1) / denom); sigma = 1.0 / denom; } rpi = 2.0 * norm_val / M_PI; for (k = 0; k < nm; k++) if (!GMT_is_fnan (data[k])) data[k] = (float)(rpi * atan((data[k] - ave_gradient)*denom)); header.z_max = rpi * atan((max_gradient - ave_gradient)*denom); header.z_min = rpi * atan((min_gradient - ave_gradient)*denom); } else if (exp_trans) { if (!sigma_set) { sigma = 0.0; for (k = 0; k < nm; k++) if (!GMT_is_fnan (data[k])) sigma += fabs((double)data[k]); sigma = M_SQRT2 * sigma / n_used; } denom = M_SQRT2 / sigma; for (k = 0; k < nm; k++) { if (GMT_is_fnan (data[k])) continue; if (data[k] < ave_gradient) { data[k] = (float)(-norm_val * (1.0 - exp((data[k] - ave_gradient)*denom))); } else { data[k] = (float)(norm_val * (1.0 - exp(-(data[k] - ave_gradient)*denom))); } } header.z_max = norm_val * (1.0 - exp(-(max_gradient - ave_gradient)*denom)); header.z_min = -norm_val * (1.0 - exp((min_gradient - ave_gradient)*denom)); } else { if ( (max_gradient - ave_gradient) > (ave_gradient - min_gradient) ) { denom = norm_val / (max_gradient - ave_gradient); } else { denom = norm_val / (ave_gradient - min_gradient); } for (k = 0; k < nm; k++) if (!GMT_is_fnan (data[k])) data[k] = (float)((data[k] - ave_gradient) * denom); header.z_max = (max_gradient - ave_gradient) * denom; header.z_min = (min_gradient - ave_gradient) * denom; } } } /* Now we write out: */ if (do_direct_deriv) { if (normalize) { strcpy (header.title, "Normalized directional derivative(s)"); } else { strcpy (header.title, "Directional derivative(s)"); } sprintf (format, "\t%s\t%s\t%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format); if (gmtdefs.verbose) { fprintf (stderr, "%s: Min Mean Max sigma intensities:", GMT_program); fprintf (stderr, format, min_gradient, ave_gradient, max_gradient, sigma); } } else strcpy (header.title, "Directions of maximum slopes"); GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 0; /* Because of the shift */ GMT_write_grd (outfile, &header, data, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE); GMT_free ((void *) data); if (save_slopes) { strcpy (header.title, "Magnitude of maximum slopes"); GMT_write_grd (slopefile, &header, slp, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE); GMT_free ((void *)slp); } GMT_end (argc, argv); }