/* $Id: img2grd.c 12857 2014-02-06 00:26:26Z pwessel $ * * Copyright (c) 1991-2014 by P. Wessel, W. H. F. Smith, R. Scharroo, J. Luis and F. Wobbe * See LICENSE.TXT file for copying and redistribution conditions. * * img2grd.c * * img2grd reads an "img" format file (used by Sandwell and Smith * in their estimates of gravity and depth from satellite altimetry), * extracts a Region [with optional N by N averaging], and writes out * the data [or Track control] as a pixel-registered GMT "grd" file, * preserving the Mercator projection (gmtdefaults PROJ_ELLIPSOID = Sphere) * inherent in the img file. * * img file is read from dir GMT_IMGDIR if this is defined as an * environment variable; else img file name is opened directly. * * The coordinates in the img file are determined by the latitude and * longitude span of the img file and the width of an img pixel in * minutes of longitude. These may be changed from their default values * using the lower-case [-W] [-D/] [-I] * but must be set to match the coverage of the input img file. * * The user-supplied -R/// will be adjusted by rounding * outward so that the actual range spanned by the file falls exactly on * the edges of the input pixels. If the user uses the option to * average the input pixels into N by N square block averages, then the * -R will be adjusted so that the range spans the block averages. * * The coordinates in the output file are in spherical mercator projected * units ranging from 0,0 in the lower left corner of the output to * xmax, ymax in the upper right corner, where xmax,ymax are the width * and height of the (spherical) Mercator map with -Rww/ee/ss/nn and -Jm1. * Here ww/ee/ss/nn are the outward-rounded range values. * * Further details on the coordinate system can be obtained from the * comments in gmt_imgsubs.c and gmt_imgsubs.h * * This is a complete rebuild (for v3.1) of the old program by this name. * New functionality added here is the averaging option [-N] and the * options to define -I, -W, -y, which permit handling very early versions * of these files and anticipate higher-resolutions to come down the road. * Also added is the option to look for the img file in GMT_IMGDIR * * Author: Walter H F Smith * Date: 8 October, 1998 * **WHFS 27 Nov 1998 fixed bug so that -T0 gives correct results ** PW 18 Oct 1999 Use WORDS_BIGENDIAN macro set by configure. * PW 12 Apr 2006 Replaced -x -y with -W -D * PW 28 Nov 2006 Added -C for setting origin to main img origin (0,0) * */ #define THIS_MODULE_NAME "img2grd" #define THIS_MODULE_LIB "img" #define THIS_MODULE_PURPOSE "Extract a subset from an img file in Mercator or Geographic format" #include "gmt_imgsubs.h" #include "gmt_dev.h" #include "common_byteswap.h" #define GMT_PROG_OPTIONS "-VRn" GMT_OPT("m") struct IMG2GRD_CTRL { struct In { /* Input file name */ bool active; char *file; /* Input file name */ } In; struct C { /* -C */ bool active; } C; struct D { /* -D[/] */ bool active; double min, max; } D; struct E { /* -E */ bool active; } E; struct G { /* -G */ bool active; char *file; } G; struct I { /* -I */ bool active; double value; } I; struct M { /* -M */ bool active; } M; struct N { /* -N */ bool active; int value; } N; struct S { /* -S */ bool active; unsigned int mode; double value; } S; struct T { /* -T */ bool active; int value; } T; struct W { /* -W */ bool active; double value; } W; }; void *New_img2grd_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ struct IMG2GRD_CTRL *C; C = GMT_memory (GMT, NULL, 1, struct IMG2GRD_CTRL); /* Initialize values whose defaults are not 0/false/NULL */ C->D.min = GMT_IMG_MINLAT; C->D.max = GMT_IMG_MAXLAT; C->N.value = 1; /* No averaging */ C->T.value = 1; /* Default img type */ C->I.value = GMT_IMG_MPIXEL; return (C); } void Free_img2grd_Ctrl (struct GMT_CTRL *GMT, struct IMG2GRD_CTRL *C) { /* Deallocate control structure */ if (!C) return; if (C->In.file) free (C->In.file); if (C->G.file) free (C->G.file); GMT_free (GMT, C); } int GMT_img2grd_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: img2grd %s -G -T [-C]\n", GMT_Rgeo_OPT); GMT_Message (API, GMT_TIME_NONE, "\t[-D[/]] [-E] [-I] [-M] [-N] [-S[]] [%s]\n\t[-W] [%s]\n\n", GMT_V_OPT, GMT_n_OPT); if (level == GMT_SYNOPSIS) return (EXIT_FAILURE); GMT_Message (API, GMT_TIME_NONE, "\t gives name of img file.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-G Set filename for the output grid file.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-R Specify the region in decimal degrees or degrees:minutes.\n"); GMT_Message (API, GMT_TIME_NONE, "\n\tOPTIONS:\n"); GMT_Message (API, GMT_TIME_NONE, "\t-C Refer Mercator coordinates to img source origin and requires -M\n"); GMT_Message (API, GMT_TIME_NONE, "\t [Default sets lower left to 0,0].\n"); GMT_Message (API, GMT_TIME_NONE, "\t-D Set input img file bottom and top latitudes [%.3f/%.3f].\n", GMT_IMG_MINLAT, GMT_IMG_MAXLAT); GMT_Message (API, GMT_TIME_NONE, "\t If no latitudes are given it is taken to mean %.3f/%.3f.\n", GMT_IMG_MINLAT_80, GMT_IMG_MAXLAT_80); GMT_Message (API, GMT_TIME_NONE, "\t Without -D we automatically determine the extent from the file size.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-E Resample geographic grid to the specified -R. Cannot be used with -M .\n"); GMT_Message (API, GMT_TIME_NONE, "\t (Default gives the exact -R of the Mercator grid).\n"); GMT_Message (API, GMT_TIME_NONE, "\t-I Set input img pixels to be minutes of longitude wide [2.0].\n"); GMT_Message (API, GMT_TIME_NONE, "\t Without -I we automatically determine the pixel size from the file size.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-M Write a Mercator grid [Default writes a geographic grid].\n"); GMT_Message (API, GMT_TIME_NONE, "\t-N Output averages of input in navg by navg squares [no averaging].\n"); GMT_Message (API, GMT_TIME_NONE, "\t-S Multiply img integer values by before output [1].\n"); GMT_Message (API, GMT_TIME_NONE, "\t To set scale based on information encoded in filename, just give -S.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-T Select the img type format:\n"); GMT_Message (API, GMT_TIME_NONE, "\t -T0 for obsolete img files w/ no constraint code, gets data.\n"); GMT_Message (API, GMT_TIME_NONE, "\t -T1 for new img file w/ constraints coded, gets data at all points [Default].\n"); GMT_Message (API, GMT_TIME_NONE, "\t -T2 for new img file w/ constraints coded, gets data only at constrained points, NaN elsewhere.\n"); GMT_Message (API, GMT_TIME_NONE, "\t -T3 for new img file w/ constraints coded, gets 1 at constraints, 0 elsewhere.\n"); GMT_Option (API, "V"); GMT_Message (API, GMT_TIME_NONE, "\t-W Input img file runs from 0 to longitude [360.0].\n"); GMT_Option (API, "n,."); return (EXIT_FAILURE); } int GMT_img2grd_parse (struct GMT_CTRL *GMT, struct IMG2GRD_CTRL *Ctrl, struct GMT_OPTION *options) { /* This parses the options provided to grdcut 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; struct GMTAPI_CTRL *API = GMT->parent; for (opt = options; opt; opt = opt->next) { /* Process all the options given */ switch (opt->option) { /* Common parameters */ case '<': /* Input files */ if (n_files++ == 0 && GMT_check_filearg (GMT, '<', opt->arg, GMT_IN)) Ctrl->In.file = strdup (opt->arg); else n_errors++; break; /* Processes program-specific parameters */ case 'C': Ctrl->C.active = true; break; case 'D': Ctrl->D.active = true; if (opt->arg[0] && (sscanf (opt->arg, "%lf/%lf", &Ctrl->D.min, &Ctrl->D.max)) != 2) { n_errors++; GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -D: Failed to decode /.\n"); } else { Ctrl->D.min = GMT_IMG_MINLAT_80; Ctrl->D.max = GMT_IMG_MAXLAT_80; } break; case 'E': Ctrl->E.active = true; 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 'm': if (GMT_compat_check (GMT, 4)) /* Warn and fall through */ GMT_Report (API, GMT_MSG_COMPAT, "Warning: -m is deprecated; use -I instead.\n"); else { n_errors += GMT_default_error (GMT, opt->option); break; } case 'I': Ctrl->I.active = true; if ((sscanf (opt->arg, "%lf", &Ctrl->I.value)) != 1) { n_errors++; GMT_Report (API, GMT_MSG_NORMAL, "Syntax error: -I requires a positive value.\n"); } break; case 'M': Ctrl->M.active = true; break; case 'N': Ctrl->N.active = true; if ((sscanf (opt->arg, "%d", &Ctrl->N.value)) != 1 || Ctrl->N.value < 1) { n_errors++; GMT_Report (API, GMT_MSG_NORMAL, "Syntax error: -N requires an integer > 1.\n"); } break; case 'S': Ctrl->S.active = true; if (sscanf (opt->arg, "%lf", &Ctrl->S.value) != 1) Ctrl->S.mode = 1; break; case 'T': Ctrl->T.active = true; if ((sscanf (opt->arg, "%d", &Ctrl->T.value)) != 1) { n_errors++; GMT_Report (API, GMT_MSG_NORMAL, "Syntax error: -T requires an output type 0-3.\n"); } break; case 'W': Ctrl->W.active = true; if ((sscanf (opt->arg, "%lf", &Ctrl->W.value)) != 1) { n_errors++; GMT_Report (API, GMT_MSG_NORMAL, "Syntax error: -W requires a longitude >= 360.0.\n"); } break; default: /* Report bad options */ n_errors += GMT_default_error (GMT, opt->option); break; } } n_errors += GMT_check_condition (GMT, Ctrl->In.file == NULL, "Syntax error: Must specify input imgfile name.\n"); n_errors += GMT_check_condition (GMT, n_files > 1, "Syntax error: More than one world image file name given.\n"); n_errors += GMT_check_condition (GMT, !GMT->common.R.active, "Syntax error: Must specify -R option.\n"); n_errors += GMT_check_condition (GMT, GMT->common.R.active && (GMT->common.R.wesn[XLO] >= GMT->common.R.wesn[XHI] || GMT->common.R.wesn[YLO] >= GMT->common.R.wesn[YHI]), "Syntax error:Must specify -R with west < east and south < north.\n"); n_errors += GMT_check_condition (GMT, !Ctrl->G.active || Ctrl->G.file == NULL, "Syntax error: Must specify output grid file name with -G.\n"); n_errors += GMT_check_condition (GMT, Ctrl->D.active && (Ctrl->D.min <= -90 || Ctrl->D.max >= 90.0 || Ctrl->D.max <= Ctrl->D.min), "Syntax error: Min/max latitudes are invalid.\n"); n_errors += GMT_check_condition (GMT, Ctrl->T.value < 0 || Ctrl->T.value > 3, "Syntax error: Must specify output type in the range 0-3 with -T.\n"); n_errors += GMT_check_condition (GMT, Ctrl->W.active && Ctrl->W.value < 360.0, "Syntax error: Requires a maximum longitude >= 360.0 with -W.\n"); n_errors += GMT_check_condition (GMT, Ctrl->I.active && Ctrl->I.value <= 0.0, "Syntax error: Requires a positive value with -I.\n"); n_errors += GMT_check_condition (GMT, Ctrl->C.active && !Ctrl->M.active, "Syntax error: -C requires -M.\n"); n_errors += GMT_check_condition (GMT, Ctrl->E.active && Ctrl->M.active, "Syntax error: -E cannot be used with -M.\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_OK); } #define bailout(code) {GMT_Free_Options (mode); return (code);} #define Return(code) {Free_img2grd_Ctrl (GMT, Ctrl); GMT_end_module (GMT, GMT_cpy); bailout (code);} int GMT_img2grd (void *V_API, int mode, void *args) { int error = 0; unsigned int navgsq, navg; /* navg by navg pixels are averaged if navg > 1; else if navg == 1 do nothing */ unsigned int nx, ny, iout, jout, jinstart, jinstop, k, kk, ion, jj, iin, jin2, ii, kstart, *ix = NULL; int in_ID, out_ID = GMT_NOTSET, jin, iinstart, iinstop; uint64_t ij; int16_t tempint; size_t n_expected; /* Expected items per row */ double west, east, south, north, wesn[4], toplat, botlat, dx; double south2, north2, rnavgsq, csum, dsum, left, bottom, inc[2]; int16_t *row = NULL; uint16_t *u2 = NULL; char infile[GMT_BUFSIZ] = {""}, cmd[GMT_BUFSIZ] = {""}, s_in_ID[GMT_STR16] = {""}, s_out_ID[GMT_LEN256] = {""}; char z_units[GMT_GRID_UNIT_LEN80] = {""}; FILE *fp = NULL; struct GMT_IMG_COORD imgcoord; struct GMT_IMG_RANGE imgrange = { GMT_IMG_MAXLON, GMT_IMG_MINLAT, GMT_IMG_MAXLAT, GMT_IMG_MPIXEL }; struct GMT_GRID *Merc = NULL; struct IMG2GRD_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_img2grd_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_img2grd_usage (API, GMT_USAGE)); /* Return the usage message */ if (options->option == GMT_OPT_SYNOPSIS) bailout (GMT_img2grd_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 (GMT_Parse_Common (API, GMT_PROG_OPTIONS, options)) Return (API->error); Ctrl = New_img2grd_Ctrl (GMT); /* Allocate and initialize a new control structure */ if ((error = GMT_img2grd_parse (GMT, Ctrl, options))) Return (error); /*---------------------------- This is the img2grd main code ----------------------------*/ GMT_set_pad (GMT, 2U); /* Ensure space for BCs in case an API passed pad == 0 */ /* Set up default settings if not specified */ if (Ctrl->W.active) imgrange.maxlon = Ctrl->W.value; if (Ctrl->I.active) imgrange.mpixel = Ctrl->I.value; if (Ctrl->D.active) { imgrange.minlat = Ctrl->D.min; imgrange.maxlat = Ctrl->D.max; } if (!GMT_getdatapath (GMT, Ctrl->In.file, infile, R_OK)) { GMT_Report (API, GMT_MSG_NORMAL, "img file %s not found\n", Ctrl->In.file); Return (GMT_GRDIO_FILE_NOT_FOUND); } if (! (Ctrl->I.active && Ctrl->D.active)) { int min; double lat = 0.0; struct stat buf; if (stat (infile, &buf)) return (GMT_GRDIO_STAT_FAILED); /* Inquiry about file failed somehow */ switch (buf.st_size) { /* Known sizes are 1 or 2 min at lat_max = ~72 or 80 */ case GMT_IMG_NLON_1M*GMT_IMG_NLAT_1M_80*GMT_IMG_ITEMSIZE: if (lat == 0.0) lat = GMT_IMG_MAXLAT_80; case GMT_IMG_NLON_1M*GMT_IMG_NLAT_1M_72*GMT_IMG_ITEMSIZE: if (lat == 0.0) lat = GMT_IMG_MAXLAT_72; min = 1; break; case GMT_IMG_NLON_2M*GMT_IMG_NLAT_2M_80*GMT_IMG_ITEMSIZE: if (lat == 0.0) lat = GMT_IMG_MAXLAT_80; case GMT_IMG_NLON_2M*GMT_IMG_NLAT_2M_72*GMT_IMG_ITEMSIZE: if (lat == 0.0) lat = GMT_IMG_MAXLAT_72; min = 2; break; default: if (lat == 0.0) return (GMT_GRDIO_BAD_IMG_LAT); min = (buf.st_size > GMT_IMG_NLON_2M*GMT_IMG_NLAT_2M_80*GMT_IMG_ITEMSIZE) ? 1 : 2; if (!Ctrl->I.active) GMT_Report (API, GMT_MSG_NORMAL, "img file %s has unusual size - grid increment defaults to %d min\n", infile, min); break; } if (!Ctrl->D.active) {imgrange.minlat = -lat; imgrange.maxlat = +lat;} if (!Ctrl->I.active) imgrange.mpixel = Ctrl->I.value = (double) min; GMT_Report (API, GMT_MSG_LONG_VERBOSE, "img file %s determined to have an increment of %d min and a latitude bound at +/- %g\n", infile, min, lat); } strcpy (z_units, "meters, mGal, Eotvos, micro-radians or Myr, depending on img file and -S."); if (Ctrl->S.mode) { /* Guess the scaling */ if (strstr (infile, "age")) { Ctrl->S.value = 0.01; strcpy (z_units, "Myr"); GMT_Report (API, GMT_MSG_VERBOSE, "img file %s determined to be crustal ages with scale 0.01.\n", infile); } if (strstr (infile, "topo")) { Ctrl->S.value = 1.0; strcpy (z_units, "meter"); GMT_Report (API, GMT_MSG_VERBOSE, "img file %s determined to be bathymetry with scale 1.\n", infile); } else if (strstr (infile, "grav")) { Ctrl->S.value = 0.1; strcpy (z_units, "mGal"); GMT_Report (API, GMT_MSG_VERBOSE, "img file %s determined to be free-air anomalies with scale 0.1.\n", infile); } else if (strstr (infile, "curv") || strstr (infile, "vgg")) { Ctrl->S.value = (Ctrl->I.value == 2.0) ? 0.05 : 0.02; strcpy (z_units, "Eotvos"); GMT_Report (API, GMT_MSG_VERBOSE, "img file %s determined to be VGG anomalies with scale %g.\n", infile, Ctrl->S.value); } else { GMT_Report (API, GMT_MSG_VERBOSE, "Unable to determine data type for img file %s; scale not used.\n", infile); Ctrl->S.value = 1.0; } } if (Ctrl->T.value == 3) strcpy (z_units, "T/F, one or more constraints fell in this pixel."); if (GMT_img_setup_coord (GMT, &imgrange, &imgcoord) ) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error: Error in img coordinate specification [-I -W or -D].\n"); Return (GMT_RUNTIME_ERROR); } else if (Ctrl->N.value && (imgcoord.nx360%Ctrl->N.value != 0 || imgcoord.nyrow%Ctrl->N.value != 0) ) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error: Bad choice of navg in -N. Must divide %d and %d\n", imgcoord.nx360, imgcoord.nyrow); Return (GMT_RUNTIME_ERROR); } if ((fp = fopen (infile, "rb")) == NULL) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error: Cannot open %s for binary read.\n", infile); Return (GMT_RUNTIME_ERROR); } GMT_set_geographic (GMT, GMT_IN); GMT_set_cartesian (GMT, GMT_OUT); /* Since output is no longer lon/lat */ /* Keep original -R for later */ west = wesn[XLO] = GMT->common.R.wesn[XLO]; east = wesn[XHI] = GMT->common.R.wesn[XHI]; south = wesn[YLO] = GMT->common.R.wesn[YLO]; north = wesn[YHI] = GMT->common.R.wesn[YHI]; navg = Ctrl->N.value; /* Expected edges of input image based on coordinate initialization (might not exactly match user spec) */ toplat = GMT_img_ypix_to_lat (0.0, &imgcoord); botlat = GMT_img_ypix_to_lat ((double)imgcoord.nyrow, &imgcoord); dx = 1.0 / ((double)imgcoord.nx360 / 360.0); inc[GMT_X] = inc[GMT_Y] = dx * navg; GMT_Report (API, GMT_MSG_VERBOSE, "Expects %s to be %d by %d pixels spanning 0/%5.1f/%.8g/%.8g.\n", infile, imgcoord.nxcol, imgcoord.nyrow, dx*imgcoord.nxcol, botlat, toplat); if (toplat < wesn[YHI]) { GMT_Report (API, GMT_MSG_VERBOSE, "Warning: Your top latitude (%.12g) lies outside top latitude of input (%.12g) - now truncated.\n", wesn[YHI], toplat); wesn[YHI] = toplat - GMT_CONV_LIMIT; /* To ensure proper round-off in calculating ny */ } if (botlat > wesn[YLO]) { GMT_Report (API, GMT_MSG_VERBOSE, "Warning: Your bottom latitude (%.12g) lies outside bottom latitude of input (%.12g) - now truncated.\n", wesn[YLO], botlat); wesn[YLO] = botlat + GMT_CONV_LIMIT; /* To ensure proper round-off in calculating ny */ } /* Re-adjust user's -R so that it falls on pixel coordinate boundaries */ jinstart = navg * (urint (floor (GMT_img_lat_to_ypix (wesn[YHI], &imgcoord) / navg))); jinstop = navg * (urint (ceil (GMT_img_lat_to_ypix (wesn[YLO], &imgcoord) / navg))); /* jinstart <= jinputrow < jinstop */ ny = (int)((jinstop - jinstart) / navg); north2 = wesn[YHI] = GMT_img_ypix_to_lat ((double)jinstart, &imgcoord); south2 = wesn[YLO] = GMT_img_ypix_to_lat ((double)jinstop, &imgcoord); iinstart = navg * (urint (floor (wesn[XLO]/inc[GMT_X]))); iinstop = navg * (urint (ceil (wesn[XHI]/inc[GMT_X]))); /* iinstart <= ipixelcol < iinstop, but modulo all with imgcoord.nx360 */ /* Reset left and right edges of user area */ wesn[XLO] = iinstart * dx; wesn[XHI] = iinstop * dx; nx = (int)((iinstop - iinstart) / navg); GMT_Report (API, GMT_MSG_VERBOSE, "To fit [averaged] input, your %s is adjusted to -R%.12g/%.12g/%.12g/%.12g.\n", Ctrl->In.file, wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI]); GMT_Report (API, GMT_MSG_VERBOSE, "The output grid size will be %d by %d pixels.\n", nx, ny); /* Set iinstart so that it is non-negative, for use to index pixels. */ while (iinstart < 0) iinstart += imgcoord.nx360; /* Set navgsq, rnavgsq, for the averaging */ navgsq = navg * navg; rnavgsq = 1.0 / navgsq; /* Set up header with Mercatorized dimensions assuming -Jm1 */ if (Ctrl->C.active) { int equator; equator = irint (GMT_img_lat_to_ypix (0.0, &imgcoord)); wesn[XLO] = iinstart * inc[GMT_X]; wesn[XHI] = wesn[XLO] + nx * inc[GMT_X]; wesn[YHI] = (imgcoord.nyrow - (int)jinstart - equator) * inc[GMT_Y]; wesn[YLO] = wesn[YHI] - ny * inc[GMT_Y]; left = bottom = 0.0; if (wesn[XHI] > 360.0) { wesn[XHI] -= 360.0; wesn[XLO] -= 360.0; } else if (west < 0.0 && east < 0.0 && west >= -180.0 && wesn[XLO] > 0.0) { /* Gave reasonable negative region, honor that for -C */ wesn[XHI] -= 360.0; wesn[XLO] -= 360.0; } } else { wesn[XLO] = 0.0; wesn[XHI] = nx * inc[GMT_X]; wesn[YLO] = 0.0; wesn[YHI] = ny * inc[GMT_Y]; left = wesn[XLO]; bottom = wesn[YLO]; } GMT_Report (API, GMT_MSG_DEBUG, "Allocate Grid container for Mercator data\n"); if ((Merc = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_GRID_ALL, NULL, wesn, inc, \ GMT_GRID_PIXEL_REG, GMT_NOTSET, Ctrl->G.file)) == NULL) Return (API->error); if (Ctrl->M.active) { sprintf (Merc->header->x_units, "Spherical Mercator projected Longitude, -Jm1, length from %.12g", left); sprintf (Merc->header->y_units, "Spherical Mercator projected Latitude, -Jm1, length from %.12g", bottom); sprintf (Merc->header->remark, "Spherical Mercator Projected with -Jm1 -R%.12g/%.12g/%.12g/%.12g", wesn[XLO], wesn[XHI], south2, north2); } else { sprintf (Merc->header->x_units, "longitude [degrees_east]"); sprintf (Merc->header->y_units, "latitude [degrees_north]"); } strncpy (Merc->header->z_units, z_units, GMT_GRID_UNIT_LEN80); strcpy (Merc->header->title, "Data from Altimetry"); Merc->header->z_min = DBL_MAX; Merc->header->z_max = -DBL_MAX; /* Now malloc some space for integer pixel index, and int16_t data buffer. */ row = GMT_memory (GMT, NULL, navg * imgcoord.nxcol, int16_t); ix = GMT_memory (GMT, NULL, navgsq * Merc->header->nx, unsigned int); /* Load ix with the index to the correct column, for each output desired. This helps for Greenwich, also faster averaging of the file, etc. Note for averaging each n by n block is looped in turn. */ if (navg > 1) { for (iout = k = 0; iout < Merc->header->nx; iout++) { ion = iout * navg; for (jin2 = 0; jin2 < navg; jin2++) { jj = jin2 * imgcoord.nxcol; for (iin = 0; iin < navg; iin++) { ii = (iin + iinstart + ion) % imgcoord.nx360; ix[k] = ii + jj; k++; } } } } else { for (iout = 0; iout < Merc->header->nx; iout++) ix[iout] = (iout + iinstart) % imgcoord.nx360; } /* Now before beginning data loop, fseek if needed. */ if (jinstart > 0 && abs(jinstart) < imgcoord.nyrow) fseek (fp, 2LL * imgcoord.nxcol * jinstart, SEEK_SET); /* Now loop over output points, reading and handling data as needed */ n_expected = navg * imgcoord.nxcol; for (jout = 0; jout < Merc->header->ny; jout++) { jin = jinstart + navg * jout; ij = GMT_IJP (Merc->header, jout, 0); /* Left-most index of this row */ if (jin < 0 || jin >= imgcoord.nyrow) { /* Outside latitude range; set row to NaNs */ for (iout = 0; iout < Merc->header->nx; iout++, ij++) Merc->data[ij] = GMT->session.f_NaN; continue; } if ((fread (row, sizeof (int16_t), n_expected, fp) ) != n_expected) { GMT_Report (API, GMT_MSG_NORMAL, "Error: Read failure at jin = %d.\n", jin); exit (EXIT_FAILURE); } #ifndef WORDS_BIGENDIAN u2 = (uint16_t *)row; for (iout = 0; iout < navg * imgcoord.nxcol; iout++) u2[iout] = bswap16 (u2[iout]); #endif for (iout = 0, kstart = 0; iout < Merc->header->nx; iout++, ij++, kstart += navgsq) { if (navg) { csum = dsum = 0.0; for (k = 0, kk = kstart; k < navgsq; k++, kk++) { tempint = row[ix[kk]]; if (Ctrl->T.value && abs (tempint) % 2 != 0) { csum += 1.0; tempint--; } dsum += (double) tempint; } csum *= rnavgsq; dsum *= rnavgsq; } else { tempint = row[ix[iout]]; if (Ctrl->T.value && abs (tempint) %2 != 0) { csum = 1.0; tempint--; } else csum = 0.0; dsum = (double) tempint; } if (Ctrl->S.active) dsum *= Ctrl->S.value; switch (Ctrl->T.value) { case 0: case 1: Merc->data[ij] = (float) dsum; break; case 2: Merc->data[ij] = (float)((csum >= 0.5) ? dsum : GMT->session.f_NaN); break; case 3: Merc->data[ij] = (float)csum; break; } if (Ctrl->T.value != 2 || csum >= 0.5) { if (Merc->header->z_min > Merc->data[ij]) Merc->header->z_min = Merc->data[ij]; if (Merc->header->z_max < Merc->data[ij]) Merc->header->z_max = Merc->data[ij]; } } } fclose (fp); GMT_free (GMT, ix); GMT_free (GMT, row); /* We now have the Mercator grid in Grid. */ GMT_Report (API, GMT_MSG_VERBOSE, "Created %d by %d Mercatorized grid file. Min, Max values are %.8g %.8g\n", Merc->header->nx, Merc->header->ny, Merc->header->z_min, Merc->header->z_max); if (Ctrl->M.active) { /* Write out the Mercator grid and return, no projection needed */ GMT_set_pad (GMT, API->pad); /* Reset to session default pad before output */ if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Merc)) Return (API->error); if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_ALL, NULL, Ctrl->G.file, Merc) != GMT_OK) { Return (API->error); } Return (GMT_OK); } /* Here we need to reproject the data, and Merc becomes a temporary grid */ GMT_Report (API, GMT_MSG_VERBOSE, "Undo the implicit spherical Mercator -Jm1i projection.\n"); /* Preparing source and destination for GMT_grdproject */ /* a. Register the Mercator grid to be the source read by GMT_grdproject by passing a pointer */ GMT_Report (API, GMT_MSG_DEBUG, "Register Mercator Grid container as grdproject input\n"); if ((in_ID = GMT_Register_IO (API, GMT_IS_GRID, GMT_IS_REFERENCE, GMT_IS_SURFACE, GMT_IN, NULL, Merc)) == GMT_NOTSET) { Return (API->error); } if (GMT_Encode_ID (API, s_in_ID, in_ID) != GMT_OK) { Return (API->error); /* Make filename with embedded object ID */ } /* b. If -E: Register a grid struct Geo to be the destination allocated and written to by GMT_grdproject, else write to -G */ if (Ctrl->E.active) { /* Since we will resample again, register a memory location for the result */ GMT_Report (API, GMT_MSG_DEBUG, "Register memory Grid container as grdproject output\n"); if ((out_ID = GMT_Register_IO (API, GMT_IS_GRID, GMT_IS_REFERENCE, GMT_IS_SURFACE, GMT_OUT, NULL, NULL)) == GMT_NOTSET) { Return (API->error); } if (GMT_Encode_ID (API, s_out_ID, out_ID) != GMT_OK) { Return (API->error); /* Make filename with embedded object ID */ } } else /* The output here is the final result */ strncpy (s_out_ID, Ctrl->G.file, GMT_LEN256); sprintf (cmd, "-R%g/%g/%g/%g -Jm1i -I %s -G%s --PROJ_ELLIPSOID=Sphere --PROJ_LENGTH_UNIT=inch", west, east, south2, north2, s_in_ID, s_out_ID); GMT_Report (API, GMT_MSG_DEBUG, "Calling grdproject %s.\n", cmd); if (GMT_Call_Module (API, "grdproject", GMT_MODULE_CMD, cmd)!= GMT_OK) { /* Inverse project the grid or fail */ Return (API->error); } if (GMT_Destroy_Data (API, &Merc) != GMT_OK) { Return (API->error); } if (Ctrl->E.active) { /* Resample again using the given -R and the dx/dy in even minutes */ /* Preparing source and destination for GMT_grdsample */ /* a. Register the Geographic grid returned by GMT_grdproject to be the source read by GMT_grdsample by passing a pointer */ struct GMT_GRID *Geo = NULL; GMT_Report (API, GMT_MSG_DEBUG, "Retrieve Geo Grid container as grdproject output\n"); if ((Geo = GMT_Retrieve_Data (API, out_ID)) == NULL) { Return (API->error); } strcpy (Geo->header->title, "Data from Altimetry"); strncpy (Geo->header->z_units, z_units, GMT_GRID_UNIT_LEN80); sprintf (Geo->header->x_units, "longitude [degrees_east]"); sprintf (Geo->header->y_units, "latitude [degrees_north]"); GMT_Report (API, GMT_MSG_DEBUG, "Register Geo Grid container as grdsample input\n"); if ((in_ID = GMT_Register_IO (API, GMT_IS_GRID, GMT_IS_REFERENCE, GMT_IS_SURFACE, GMT_IN, NULL, Geo)) == GMT_NOTSET) { Return (API->error); } if (GMT_Encode_ID (API, s_in_ID, in_ID) != GMT_OK) { /* Make filename with embedded object ID */ Return (API->error); } sprintf (cmd, "-R%g/%g/%g/%g -I%gm %s -G%s -fg", west, east, south, north, Ctrl->I.value, s_in_ID, Ctrl->G.file); if (GMT_Call_Module (API, "grdsample", GMT_MODULE_CMD, cmd) != GMT_OK) { /* Resample the grid or fail */ Return (API->error); } } Return (GMT_OK); }