1 /*
2 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
3 * See LICENSE.TXT file for copying and redistribution conditions.
4 *
5 * img2grd.c
6 *
7 * img2grd reads an "img" format file (used by Sandwell and Smith
8 * in their estimates of gravity and depth from satellite altimetry),
9 * extracts a Region [with optional N by N averaging], and writes out
10 * the data [or Track control] as a pixel-registered GMT "grd" file,
11 * preserving the Mercator projection (gmtdefaults PROJ_ELLIPSOID = Sphere)
12 * inherent in the img file.
13 *
14 * img file is read from dir GMT_IMGDIR if this is defined as an
15 * environment variable; else img file name is opened directly.
16 *
17 * The coordinates in the img file are determined by the latitude and
18 * longitude span of the img file and the width of an img pixel in
19 * minutes of longitude. These may be changed from their default values
20 * using the lower-case [-W<maxlon>] [-D<minlat>/<maxlat>] [-I<minutes>]
21 * but must be set to match the coverage of the input img file.
22 *
23 * The user-supplied -R<w>/<e>/<s>/<n> will be adjusted by rounding
24 * outward so that the actual range spanned by the file falls exactly on
25 * the edges of the input pixels. If the user uses the option to
26 * average the input pixels into N by N square block averages, then the
27 * -R will be adjusted so that the range spans the block averages.
28 *
29 * The coordinates in the output file are in spherical mercator projected
30 * units ranging from 0,0 in the lower left corner of the output to
31 * xmax, ymax in the upper right corner, where xmax,ymax are the width
32 * and height of the (spherical) Mercator map with -Rww/ee/ss/nn and -Jm1.
33 * Here ww/ee/ss/nn are the outward-rounded range values.
34 *
35 * Further details on the coordinate system can be obtained from the
36 * comments in gmt_imgsubs.c and gmt_imgsubs.h
37 *
38 * This is a complete rebuild (for v3.1) of the old program by this name.
39 * New functionality added here is the averaging option [-N] and the
40 * options to define -I, -W, -y, which permit handling very early versions
41 * of these files and anticipate higher-resolutions to come down the road.
42 * Also added is the option to look for the img file in GMT_IMGDIR
43 *
44 * Author: Walter H F Smith
45 * Date: 8 October, 1998
46 *
47 **WHFS 27 Nov 1998 fixed bug so that -T0 gives correct results
48 ** PW 18 Oct 1999 Use WORDS_BIGENDIAN macro set by configure.
49 * PW 12 Apr 2006 Replaced -x -y with -W -D
50 * PW 28 Nov 2006 Added -C for setting origin to main img origin (0,0)
51 * PW 13 Apr 2018 Get ready for grids with < 1min spacing
52 * PW 12 May 2019 Make -C the default, add backwards option -F for old default
53 *
54 */
55
56 #include "gmt_dev.h"
57 #include "gmt_common_byteswap.h"
58
59 #define THIS_MODULE_CLASSIC_NAME "img2grd"
60 #define THIS_MODULE_MODERN_NAME "img2grd"
61 #define THIS_MODULE_LIB "img"
62 #define THIS_MODULE_PURPOSE "Extract a subset from an img file in Mercator or Geographic format"
63 #define THIS_MODULE_KEYS "<G{,GG}"
64 #define THIS_MODULE_NEEDS "R"
65 #define THIS_MODULE_OPTIONS "-VRn" GMT_OPT("m")
66
67 /* The following values are used to initialize the default values
68 controlling the range covered by the img file and the size
69 of a pixel in minutes of longitude. The values shown here
70 are for the 2-minute files currently in use (world_grav 7.2
71 and topo_polish 6.2). */
72
73 #define GMT_IMG_MAXLON 360.0
74 #define GMT_IMG_MINLAT -72.0059773539
75 #define GMT_IMG_MAXLAT +72.0059773539
76 #define GMT_IMG_MINLAT_80 -80.7380086280
77 #define GMT_IMG_MAXLAT_80 +80.7380086280
78
79 #define GMT_IMG_MPIXEL 2.0
80
81 /* The following structure contains info corresponding to
82 * the above values, which may be altered by the user at
83 * run time via argv[][] switches: */
84
85 struct GMT_IMG_RANGE {
86 double maxlon;
87 double minlat;
88 double maxlat;
89 double mpixel;
90 };
91
92 /* The following structure contains info used to set up the
93 * coordinate transformations. These are to be determined
94 * based on the values in GMT_IMG_RANGE after argv[][] has
95 * been parsed. Two different structures will be used, one
96 * representing the input file, and the other the output,
97 * to simplify the computation of coordinates for N by N
98 * averages in the output: */
99
100 struct GMT_IMG_COORD {
101 double radius; /* # of pixels in 1 radian of longitude */
102 int nx360; /* # of pixels in 360 degrees of longtd */
103 int nxcol; /* # of columns in input img file */
104 int nyrow; /* # of rows in input img file */
105 int nytop; /* # of rows from Equator to top edge */
106 };
107
108 struct IMG2GRD_CTRL {
109 struct IMG2GRD_In { /* Input file name */
110 bool active;
111 char *file; /* Input file name */
112 } In;
113 struct IMG2GRD_D { /* -D[<minlat>/<maxlat>] */
114 bool active;
115 double min, max;
116 } D;
117 struct IMG2GRD_E { /* -E */
118 bool active;
119 } E;
120 struct IMG2GRD_F { /* -F */
121 bool active;
122 } F;
123 struct IMG2GRD_G { /* -G<output_grdfile> */
124 bool active;
125 char *file;
126 } G;
127 struct IMG2GRD_I { /* -I<minutes> */
128 bool active;
129 double value;
130 } I;
131 struct IMG2GRD_M { /* -M */
132 bool active;
133 } M;
134 struct IMG2GRD_N { /* -N<ave> */
135 bool active;
136 int value;
137 } N;
138 struct IMG2GRD_S { /* -S<scale> */
139 bool active;
140 unsigned int mode;
141 double value;
142 } S;
143 struct IMG2GRD_T { /* -T<type> */
144 bool active;
145 int value;
146 } T;
147 struct IMG2GRD_W { /* -W<maxlon> */
148 bool active;
149 double value;
150 } W;
151 };
152
New_Ctrl(struct GMT_CTRL * GMT)153 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
154 struct IMG2GRD_CTRL *C;
155
156 C = gmt_M_memory (GMT, NULL, 1, struct IMG2GRD_CTRL);
157
158 /* Initialize values whose defaults are not 0/false/NULL */
159 C->D.min = GMT_IMG_MINLAT;
160 C->D.max = GMT_IMG_MAXLAT;
161 C->N.value = 1; /* No averaging */
162 C->T.value = 1; /* Default img type */
163 C->I.value = GMT_IMG_MPIXEL;
164
165 return (C);
166 }
167
Free_Ctrl(struct GMT_CTRL * GMT,struct IMG2GRD_CTRL * C)168 static void Free_Ctrl (struct GMT_CTRL *GMT, struct IMG2GRD_CTRL *C) { /* Deallocate control structure */
169 if (!C) return;
170 gmt_M_str_free (C->In.file);
171 gmt_M_str_free (C->G.file);
172 gmt_M_free (GMT, C);
173 }
174
usage(struct GMTAPI_CTRL * API,int level)175 static int usage (struct GMTAPI_CTRL *API, int level) {
176 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
177 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
178 GMT_Usage (API, 0, "usage: %s <world_image_filename> -G%s %s [-D[<minlat>/<maxlat>]] [-E] [-F] "
179 "[-I<min>[m|s]] [-M] [-N<navg>] [-S[<scale>]] [-T<type>] [%s] [-W<maxlon>] [%s] [%s]\n", name, GMT_OUTGRID, GMT_Rgeo_OPT, GMT_V_OPT, GMT_n_OPT, GMT_PAR_OPT);
180
181 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
182
183 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
184 GMT_Usage (API, 1, "\n<world_image_filename> gives name of img file.");
185 gmt_outgrid_syntax (API, 'G', "Set name of the output grid file");
186 GMT_Usage (API, 1, "\n%s", GMT_Rgeo_OPT);
187 GMT_Usage (API, -2, "Specify the region in decimal degrees or degrees:minutes.");
188 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
189 GMT_Usage (API, 1, "\n-D[<minlat>/<maxlat>]");
190 GMT_Usage (API, -2, "Set input img file bottom and top latitudes [%.3f/%.3f]. "
191 "If no latitudes are given it is taken to mean %.3f/%.3f. "
192 "Without -D we automatically determine the extent from the file size.",
193 GMT_IMG_MINLAT, GMT_IMG_MAXLAT, GMT_IMG_MINLAT_80, GMT_IMG_MAXLAT_80);
194 GMT_Usage (API, 1, "\n-E Resample geographic grid to the specified -R. Cannot be used with -M "
195 "[Default gives the exact -R of the Mercator grid].");
196 GMT_Usage (API, 1, "\n-F Translate Mercator coordinates so lower left is (0,0); requires -M.");
197 GMT_Usage (API, 1, "\n-I<min>[m|s]");
198 GMT_Usage (API, -2, "Set input img pixels to be <min> minutes of longitude wide [2.0]. "
199 "Without -I we automatically determine the pixel size from the file size.");
200 GMT_Usage (API, 1, "\n-M Write a Mercator grid [Default writes a geographic grid].");
201 GMT_Usage (API, 1, "\n-N<navg>");
202 GMT_Usage (API, -2, "Output averages of input in <navg> by <navg> squares [no averaging].");
203 GMT_Usage (API, 1, "\n-S[<scale>]");
204 GMT_Usage (API, -2, "Multiply img integer values by <scale> before output [1]. "
205 "To set scale based on information encoded in filename, just give -S.");
206 GMT_Usage (API, 1, "\n-T<type>");
207 GMT_Usage (API, -2, "Select the img type format:");
208 GMT_Usage (API, 3, "0: Obsolete img files w/ no constraint code, gets data.");
209 GMT_Usage (API, 3, "1: New img file w/ constraints coded, gets data at all points [Default].");
210 GMT_Usage (API, 3, "2: New img file w/ constraints coded, gets data only at constrained points, NaN elsewhere.");
211 GMT_Usage (API, 3, "3: New img file w/ constraints coded, gets 1 at constraints, 0 elsewhere.");
212 GMT_Option (API, "V");
213 GMT_Usage (API, 1, "\n-W<maxlon>");
214 GMT_Usage (API, -2, "Input img file runs from 0 to <maxlon> longitude [360.0].");
215 GMT_Option (API, "n,.");
216
217 return (GMT_MODULE_USAGE);
218 }
219
parse(struct GMT_CTRL * GMT,struct IMG2GRD_CTRL * Ctrl,struct GMT_OPTION * options)220 static int parse (struct GMT_CTRL *GMT, struct IMG2GRD_CTRL *Ctrl, struct GMT_OPTION *options) {
221
222 /* This parses the options provided to grdcut and sets parameters in CTRL.
223 * Any GMT common options will override values set previously by other commands.
224 * It also replaces any file names specified as input or output with the data ID
225 * returned when registering these sources/destinations with the API.
226 */
227
228 unsigned int n_errors = 0, n_files = 0;
229 bool sec = false, min = false;
230 size_t L = 0;
231 struct GMT_OPTION *opt = NULL;
232 struct GMTAPI_CTRL *API = GMT->parent;
233
234 for (opt = options; opt; opt = opt->next) { /* Process all the options given */
235
236 switch (opt->option) {
237 /* Common parameters */
238
239 case '<': /* Input files */
240 if (n_files++ > 0) break;
241 Ctrl->In.active = true;
242 if (opt->arg[0]) Ctrl->In.file = strdup (opt->arg);
243 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file))) n_errors++;
244 break;
245
246 /* Processes program-specific parameters */
247
248 case 'F': /* Backwards helper for old non-C behavior */
249 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
250 Ctrl->F.active = true;
251 break;
252 case 'D':
253 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
254 Ctrl->D.active = true;
255 if (opt->arg[0] && (sscanf (opt->arg, "%lf/%lf", &Ctrl->D.min, &Ctrl->D.max)) != 2) {
256 n_errors++;
257 GMT_Report (API, GMT_MSG_ERROR, "Option -D: Failed to decode <minlat>/<maxlat>.\n");
258 }
259 else {
260 Ctrl->D.min = GMT_IMG_MINLAT_80;
261 Ctrl->D.max = GMT_IMG_MAXLAT_80;
262 }
263 break;
264 case 'E':
265 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
266 Ctrl->E.active = true;
267 break;
268 case 'G':
269 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
270 Ctrl->G.active = true;
271 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
272 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
273 break;
274 case 'm':
275 if (gmt_M_compat_check (GMT, 4)) /* Warn and fall through to 'I' on purpose */
276 GMT_Report (API, GMT_MSG_COMPAT, "-m<inc> is deprecated; use -I<inc> instead.\n");
277 else {
278 n_errors += gmt_default_error (GMT, opt->option);
279 break;
280 }
281 /* Intentionally fall through */
282 case 'I':
283 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
284 Ctrl->I.active = true;
285 L = strlen (opt->arg);
286 if (strchr ("ms", opt->arg[L])) { /* Valid minute or second unit */
287 if (opt->arg[L] == 's') sec = true;
288 if (opt->arg[L] == 'm') min = true;
289 opt->arg[L] = '\0';
290 }
291 if ((sscanf (opt->arg, "%lf", &Ctrl->I.value)) != 1) {
292 n_errors++;
293 GMT_Report (API, GMT_MSG_ERROR, "Option -I requires a positive value.\n");
294 }
295 if (sec) {
296 Ctrl->I.value /= 60.0;
297 opt->arg[L] = 's';
298 }
299 else if (min)
300 opt->arg[L] = 'm';
301 break;
302 case 'M':
303 n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active);
304 Ctrl->M.active = true;
305 break;
306 case 'N':
307 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
308 Ctrl->N.active = true;
309 Ctrl->N.value = atoi (opt->arg);
310 break;
311 case 'S':
312 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
313 Ctrl->S.active = true;
314 if (sscanf (opt->arg, "%lf", &Ctrl->S.value) != 1) Ctrl->S.mode = 1;
315 break;
316 case 'T':
317 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
318 Ctrl->T.active = true;
319 if ((sscanf (opt->arg, "%d", &Ctrl->T.value)) != 1) {
320 n_errors++;
321 GMT_Report (API, GMT_MSG_ERROR, "Option -T requires an output type 0-3.\n");
322 }
323 break;
324 case 'W':
325 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
326 Ctrl->W.active = true;
327 if ((sscanf (opt->arg, "%lf", &Ctrl->W.value)) != 1) {
328 n_errors++;
329 GMT_Report (API, GMT_MSG_ERROR, "Option -W requires a longitude >= 360.0.\n");
330 }
331 break;
332 default: /* Report bad options */
333 n_errors += gmt_default_error (GMT, opt->option);
334 break;
335 }
336 }
337
338 n_errors += gmt_M_check_condition (GMT, Ctrl->In.file == NULL, "Must specify input imgfile name.\n");
339 n_errors += gmt_M_check_condition (GMT, n_files > 1, "More than one world image file name given.\n");
340 n_errors += gmt_M_check_condition (GMT, !GMT->common.R.active[RSET], "Must specify -R option.\n");
341 n_errors += gmt_M_check_condition (GMT, GMT->common.R.active[RSET] && (GMT->common.R.wesn[XLO] >= GMT->common.R.wesn[XHI] || GMT->common.R.wesn[YLO] >= GMT->common.R.wesn[YHI]), "Must specify -R with west < east and south < north.\n");
342 n_errors += gmt_M_check_condition (GMT, !Ctrl->G.active || Ctrl->G.file == NULL, "Must specify output grid file name with -G.\n");
343 n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && (Ctrl->D.min <= -90 || Ctrl->D.max >= 90.0 || Ctrl->D.max <= Ctrl->D.min), "Min/max latitudes are invalid.\n");
344 n_errors += gmt_M_check_condition (GMT, Ctrl->T.value < 0 || Ctrl->T.value > 3, "Must specify output type in the range 0-3 with -T.\n");
345 n_errors += gmt_M_check_condition (GMT, Ctrl->W.active && Ctrl->W.value < 360.0, "Requires a maximum longitude >= 360.0 with -W.\n");
346 n_errors += gmt_M_check_condition (GMT, Ctrl->I.active && Ctrl->I.value <= 0.0, "Requires a positive value with -I.\n");
347 n_errors += gmt_M_check_condition (GMT, Ctrl->F.active && !Ctrl->M.active, "Option -F requires -M.\n");
348 n_errors += gmt_M_check_condition (GMT, Ctrl->E.active && Ctrl->M.active, "Option -E cannot be used with -M.\n");
349 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && Ctrl->N.value < 1, "Option -N requires an integer > 1.\n");
350
351 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
352 }
353
img2grd_gud_fwd(double y)354 GMT_LOCAL double img2grd_gud_fwd (double y) {
355 /* The Forward Gudermannian function. Given y, the distance
356 * from the Equator to a latitude on a spherical Mercator map
357 * developed from a sphere of unit radius, returns the latitude
358 * in radians. Should be called with -oo < y < +oo. Returned
359 * value will be in -M_PI_2 < value < +M_PI_2. */
360
361 return(2.0 * atan(exp(y)) - M_PI_2);
362 }
363
img2grd_gud_inv(double phi)364 GMT_LOCAL double img2grd_gud_inv (double phi) {
365 /* The Inverse Gudermannian function. Given phi, a latitude
366 * in radians, returns the distance from the Equator to this
367 * latitude on a Mercator map tangent to a sphere of unit
368 * radius. Should be called with -M_PI_2 < phi < +M_PI_2.
369 * Returned value will be in -oo < value < +oo. */
370
371 return(log(tan(M_PI_4 + 0.5 * phi)));
372 }
373
img2grd_lat_to_ypix(double lat,struct GMT_IMG_COORD * coord)374 GMT_LOCAL double img2grd_lat_to_ypix (double lat, struct GMT_IMG_COORD *coord) {
375 /* Given Latitude in degrees and pointer to coordinate struct,
376 * return (double) coordinate from top edge of input img file
377 * measured downward in coordinate pixels. */
378
379 return(coord->nytop - coord->radius * img2grd_gud_inv(lat*D2R));
380 }
381
img2grd_ypix_to_lat(double ypix,struct GMT_IMG_COORD * coord)382 GMT_LOCAL double img2grd_ypix_to_lat (double ypix, struct GMT_IMG_COORD *coord) {
383 /* Given Y coordinate, measured downward from top edge of
384 * input img file in pixels, and pointer to coordinate struct,
385 * return Latitude in degrees. */
386
387 return( R2D * img2grd_gud_fwd( (coord->nytop - ypix) / coord->radius) );
388 }
389
img2grd_setup_coord(struct GMT_CTRL * GMT,struct GMT_IMG_RANGE * r,struct GMT_IMG_COORD * c)390 GMT_LOCAL int img2grd_setup_coord (struct GMT_CTRL *GMT, struct GMT_IMG_RANGE *r, struct GMT_IMG_COORD *c) {
391 /* Given the RANGE info, set up the COORD values. Return (-1) on failure;
392 * 0 on success. */
393
394 if (r->maxlon < 360.0) {
395 GMT_Report (GMT->parent, GMT_MSG_ERROR, "ERROR from img2grd_setup_coord: Cannot handle maxlon < 360.\n");
396 return (-1);
397 }
398
399 c->nxcol = irint (r->maxlon * 60.0 / r->mpixel);
400 c->nx360 = irint (360.0 * 60.0 / r->mpixel);
401 c->radius = c->nx360 / (2.0 * M_PI);
402 c->nytop = irint (c->radius * img2grd_gud_inv(r->maxlat*D2R) );
403 c->nyrow = c->nytop - irint (c->radius * img2grd_gud_inv(r->minlat*D2R) );
404
405 return (0);
406 }
407
408 #define bailout(code) {gmt_M_free_options (mode); return (code);}
409 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
410
GMT_img2grd(void * V_API,int mode,void * args)411 EXTERN_MSC int GMT_img2grd (void *V_API, int mode, void *args) {
412 int error = 0;
413 unsigned int navgsq, navg; /* navg by navg pixels are averaged if navg > 1; else if navg == 1 do nothing */
414 unsigned int first, n_columns, n_rows, iout, jout, jinstart, jinstop, k, kk, ion, jj, iin, jin2, ii, kstart, *ix = NULL;
415 int jin, iinstart, iinstop;
416
417 uint64_t ij;
418 int16_t tempint;
419
420 size_t n_expected; /* Expected items per row */
421
422 double west, east, south, north, wesn[4], toplat, botlat, dx;
423 double south2, north2, rnavgsq, csum, dsum, left, bottom, inc[2];
424
425 int16_t *row = NULL;
426 uint16_t *u2 = NULL;
427
428 char infile[PATH_MAX] = {""}, cmd[GMT_BUFSIZ] = {""}, input[GMT_VF_LEN] = {""}, output[PATH_MAX] = {""};
429 char z_units[GMT_GRID_UNIT_LEN80] = {""}, exact_R[GMT_LEN256] = {""};
430
431 FILE *fp = NULL;
432
433 struct GMT_IMG_COORD imgcoord;
434 struct GMT_IMG_RANGE imgrange = { GMT_IMG_MAXLON, GMT_IMG_MINLAT, GMT_IMG_MAXLAT, GMT_IMG_MPIXEL };
435 struct GMT_GRID *Merc = NULL;
436 struct IMG2GRD_CTRL *Ctrl = NULL;
437 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
438 struct GMT_OPTION *options = NULL;
439 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
440
441 /*----------------------- Standard module initialization and parsing ----------------------*/
442
443 if (API == NULL) return (GMT_NOT_A_SESSION);
444 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
445 options = GMT_Create_Options (API, mode, args);
446 if (API->error)
447 return (API->error); /* Set or get option list */
448
449 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
450
451 /* Parse the command-line arguments */
452
453 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 */
454 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
455 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
456 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
457
458 /*---------------------------- This is the img2grd main code ----------------------------*/
459
460 gmt_set_pad (GMT, 2U); /* Ensure space for BCs in case an API passed pad == 0 */
461
462 /* Set up default settings if not specified */
463
464 if (Ctrl->W.active) imgrange.maxlon = Ctrl->W.value;
465 if (Ctrl->I.active) imgrange.mpixel = Ctrl->I.value;
466 if (Ctrl->D.active) {
467 imgrange.minlat = Ctrl->D.min;
468 imgrange.maxlat = Ctrl->D.max;
469 }
470
471 first = gmt_download_file_if_not_found (GMT, Ctrl->In.file, 0); /* Deal with downloadable GMT data sets first */
472
473 if (!gmt_getdatapath (GMT, &Ctrl->In.file[first], infile, R_OK)) {
474 GMT_Report (API, GMT_MSG_ERROR, "img file %s not found\n", &Ctrl->In.file[first]);
475 Return (GMT_GRDIO_FILE_NOT_FOUND);
476 }
477
478 if (! (Ctrl->I.active && Ctrl->D.active)) {
479 int min;
480 double lat = 0.0;
481 struct stat buf;
482 if (stat (infile, &buf)) Return (GMT_GRDIO_STAT_FAILED); /* Inquiry about file failed somehow */
483
484 switch (buf.st_size) { /* Known sizes are 1 or 2 min at lat_max = ~72 or 80 */
485 case GMT_IMG_NLON_1M*GMT_IMG_NLAT_1M_80*GMT_IMG_ITEMSIZE:
486 if (lat == 0.0) lat = GMT_IMG_MAXLAT_80;
487 min = 1;
488 break;
489 case GMT_IMG_NLON_1M*GMT_IMG_NLAT_1M_72*GMT_IMG_ITEMSIZE:
490 if (lat == 0.0) lat = GMT_IMG_MAXLAT_72;
491 min = 1;
492 break;
493 case GMT_IMG_NLON_2M*GMT_IMG_NLAT_2M_80*GMT_IMG_ITEMSIZE:
494 if (lat == 0.0) lat = GMT_IMG_MAXLAT_80;
495 min = 2;
496 break;
497 case GMT_IMG_NLON_2M*GMT_IMG_NLAT_2M_72*GMT_IMG_ITEMSIZE:
498 if (lat == 0.0) lat = GMT_IMG_MAXLAT_72;
499 min = 2;
500 break;
501 case GMT_IMG_NLON_4M*GMT_IMG_NLAT_4M_72*GMT_IMG_ITEMSIZE: /* Test grid only */
502 if (lat == 0.0) lat = GMT_IMG_MAXLAT_72;
503 min = 4;
504 break;
505 default:
506 if (lat == 0.0) return (GMT_GRDIO_BAD_IMG_LAT);
507 min = (buf.st_size > GMT_IMG_NLON_2M*GMT_IMG_NLAT_2M_80*GMT_IMG_ITEMSIZE) ? 1 : 2;
508 if (!Ctrl->I.active) GMT_Report (API, GMT_MSG_ERROR, "img file %s has unusual size - grid increment defaults to %d min\n", infile, min);
509 break;
510 }
511 if (!Ctrl->D.active) {imgrange.minlat = -lat; imgrange.maxlat = +lat;}
512 if (!Ctrl->I.active) imgrange.mpixel = Ctrl->I.value = (double) min;
513 GMT_Report (API, GMT_MSG_INFORMATION, "img file %s determined to have an increment of %d min and a latitude bound at +/- %g\n", infile, min, lat);
514 }
515
516 strcpy (z_units, "meters, mGal, Eotvos, micro-radians or Myr, depending on img file and -S.");
517 if (Ctrl->S.mode) { /* Guess the scaling */
518 if (strstr (infile, "age")) {
519 Ctrl->S.value = 0.01;
520 strcpy (z_units, "Myr");
521 GMT_Report (API, GMT_MSG_INFORMATION, "img file %s determined to be crustal ages with scale 0.01.\n", infile);
522 }
523 if (strstr (infile, "topo")) {
524 Ctrl->S.value = 1.0;
525 strcpy (z_units, "meter");
526 GMT_Report (API, GMT_MSG_INFORMATION, "img file %s determined to be bathymetry with scale 1.\n", infile);
527 }
528 else if (strstr (infile, "grav")) {
529 Ctrl->S.value = 0.1;
530 strcpy (z_units, "mGal");
531 GMT_Report (API, GMT_MSG_INFORMATION, "img file %s determined to be free-air anomalies with scale 0.1.\n", infile);
532 }
533 else if (strstr (infile, "curv") || strstr (infile, "vgg")) {
534 Ctrl->S.value = (Ctrl->I.value == 2.0) ? 0.05 : 0.02;
535 strcpy (z_units, "Eotvos");
536 GMT_Report (API, GMT_MSG_INFORMATION, "img file %s determined to be VGG anomalies with scale %g.\n", infile, Ctrl->S.value);
537 }
538 else {
539 GMT_Report (API, GMT_MSG_WARNING, "Unable to determine data type for img file %s; scale not used.\n", infile);
540 Ctrl->S.value = 1.0;
541 }
542 }
543 if (Ctrl->T.value == 3) strcpy (z_units, "T/F, one or more constraints fell in this pixel.");
544
545 if (img2grd_setup_coord (GMT, &imgrange, &imgcoord) ) {
546 GMT_Report (API, GMT_MSG_ERROR, "Failure in img coordinate specification [-I -W or -D].\n");
547 Return (GMT_RUNTIME_ERROR);
548 }
549 else if (Ctrl->N.active && (imgcoord.nx360%Ctrl->N.value != 0 || imgcoord.nyrow%Ctrl->N.value != 0) ) {
550 GMT_Report (API, GMT_MSG_ERROR, "Option -N: Bad choice of navg. Must divide %d and %d\n", imgcoord.nx360, imgcoord.nyrow);
551 Return (GMT_RUNTIME_ERROR);
552 }
553
554 if ((fp = gmt_fopen (GMT, infile, "rb")) == NULL) {
555 GMT_Report (API, GMT_MSG_ERROR, "Cannot open %s for binary read.\n", infile);
556 Return (GMT_RUNTIME_ERROR);
557 }
558
559 gmt_set_geographic (GMT, GMT_IN);
560 gmt_set_cartesian (GMT, GMT_OUT); /* Since output is no longer lon/lat */
561
562 /* Keep original -R for later */
563
564 west = wesn[XLO] = GMT->common.R.wesn[XLO];
565 east = wesn[XHI] = GMT->common.R.wesn[XHI];
566 south = wesn[YLO] = GMT->common.R.wesn[YLO];
567 north = wesn[YHI] = GMT->common.R.wesn[YHI];
568
569 navg = Ctrl->N.value;
570
571 /* Expected edges of input image based on coordinate initialization (might not exactly match user spec) */
572
573 toplat = img2grd_ypix_to_lat (0.0, &imgcoord);
574 botlat = img2grd_ypix_to_lat ((double)imgcoord.nyrow, &imgcoord);
575 dx = 1.0 / ((double)imgcoord.nx360 / 360.0);
576 inc[GMT_X] = inc[GMT_Y] = dx * navg;
577
578 GMT_Report (API, GMT_MSG_INFORMATION, "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);
579
580 if (toplat < wesn[YHI]) {
581 GMT_Report (API, GMT_MSG_WARNING, "Your top latitude (%.12g) lies outside top latitude of input (%.12g) - now truncated.\n", wesn[YHI], toplat);
582 wesn[YHI] = toplat - GMT_CONV8_LIMIT; /* To ensure proper round-off in calculating n_rows */
583 }
584 if (botlat > wesn[YLO]) {
585 GMT_Report (API, GMT_MSG_WARNING, "Your bottom latitude (%.12g) lies outside bottom latitude of input (%.12g) - now truncated.\n", wesn[YLO], botlat);
586 wesn[YLO] = botlat + GMT_CONV8_LIMIT; /* To ensure proper round-off in calculating n_rows */
587 }
588
589 /* Re-adjust user's -R so that it falls on pixel coordinate boundaries */
590
591 jinstart = navg * (urint (floor (img2grd_lat_to_ypix (wesn[YHI], &imgcoord) / navg)));
592 jinstop = navg * (urint (ceil (img2grd_lat_to_ypix (wesn[YLO], &imgcoord) / navg)));
593 /* jinstart <= jinputrow < jinstop */
594 n_rows = (int)(jinstop - jinstart) / navg;
595 north2 = wesn[YHI] = img2grd_ypix_to_lat ((double)jinstart, &imgcoord);
596 south2 = wesn[YLO] = img2grd_ypix_to_lat ((double)jinstop, &imgcoord);
597
598 iinstart = navg * (urint (floor (wesn[XLO]/inc[GMT_X])));
599 iinstop = navg * (urint (ceil (wesn[XHI]/inc[GMT_X])));
600 /* iinstart <= ipixelcol < iinstop, but modulo all with imgcoord.nx360 */
601 /* Reset left and right edges of user area */
602 wesn[XLO] = iinstart * dx;
603 wesn[XHI] = iinstop * dx;
604 n_columns = (int)(iinstop - iinstart) / navg;
605
606 sprintf (exact_R, "%.16g/%.16g/%.16g/%.16g", wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI]);
607 GMT_Report (API, GMT_MSG_INFORMATION, "To fit [averaged] input, your %s is adjusted to the exact region -R%s.\n", Ctrl->In.file, exact_R);
608 GMT_Report (API, GMT_MSG_INFORMATION, "The output grid size will be %d by %d pixels.\n", n_columns, n_rows);
609
610 /* Set iinstart so that it is non-negative, for use to index pixels. */
611 while (iinstart < 0) iinstart += imgcoord.nx360;
612
613 /* Set navgsq, rnavgsq, for the averaging */
614 navgsq = navg * navg;
615 rnavgsq = 1.0 / navgsq;
616
617 /* Set up header with Mercatorized dimensions assuming -Jm1i */
618 if (Ctrl->F.active || !Ctrl->M.active) { /* Backwards support for old behavior in -M, or being internally consistent with central meridian */
619 wesn[XLO] = 0.0;
620 wesn[XHI] = n_columns * inc[GMT_X];
621 wesn[YLO] = 0.0;
622 wesn[YHI] = n_rows * inc[GMT_Y];
623 left = wesn[XLO];
624 bottom = wesn[YLO];
625 }
626 else {
627 int equator = irint (img2grd_lat_to_ypix (0.0, &imgcoord));
628 wesn[XLO] = iinstart * inc[GMT_X];
629 wesn[XHI] = wesn[XLO] + n_columns * inc[GMT_X];
630 wesn[YHI] = (imgcoord.nyrow - (int)jinstart - equator) * inc[GMT_Y];
631 wesn[YLO] = wesn[YHI] - n_rows * inc[GMT_Y];
632 left = bottom = 0.0;
633 if (wesn[XHI] > 360.0) {
634 wesn[XHI] -= 360.0;
635 wesn[XLO] -= 360.0;
636 }
637 else if (west < 0.0 && east < 0.0 && west >= -180.0 && wesn[XLO] > 0.0) { /* Gave reasonable negative region, honor it */
638 wesn[XHI] -= 360.0;
639 wesn[XLO] -= 360.0;
640 }
641 }
642 GMT_Report (API, GMT_MSG_DEBUG, "Allocate Grid container for Mercator data\n");
643 if ((Merc = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, wesn, inc, \
644 GMT_GRID_PIXEL_REG, GMT_NOTSET, NULL)) == NULL) {
645 fclose (fp);
646 Return (API->error);
647 }
648
649 if (Ctrl->M.active) {
650 snprintf (Merc->header->x_units, GMT_GRID_UNIT_LEN80, "Spherical Mercator projected Longitude, -Jm1, length from %.12g", left);
651 snprintf (Merc->header->y_units, GMT_GRID_UNIT_LEN80, "Spherical Mercator projected Latitude, -Jm1, length from %.12g", bottom);
652 snprintf (Merc->header->remark, GMT_GRID_REMARK_LEN160, "%s%.12g/%.12g/%.12g/%.12g", GMT_IMG_REMARK, wesn[XLO], wesn[XHI], south2, north2);
653 }
654 else {
655 sprintf (Merc->header->x_units, "longitude [degrees_east]");
656 sprintf (Merc->header->y_units, "latitude [degrees_north]");
657 }
658 strncpy (Merc->header->z_units, z_units, GMT_GRID_UNIT_LEN80-1);
659 strcpy (Merc->header->title, "Data from Altimetry");
660 Merc->header->z_min = DBL_MAX; Merc->header->z_max = -DBL_MAX;
661 /* Now malloc some space for integer pixel index, and int16_t data buffer. */
662
663 row = gmt_M_memory (GMT, NULL, navg * imgcoord.nxcol, int16_t);
664 ix = gmt_M_memory (GMT, NULL, navgsq * Merc->header->n_columns, unsigned int);
665
666 /* Load ix with the index to the correct column, for each output desired. This helps for Greenwich,
667 also faster averaging of the file, etc. Note for averaging each n by n block is looped in turn. */
668
669 if (navg > 1) {
670 for (iout = k = 0; iout < Merc->header->n_columns; iout++) {
671 ion = iout * navg;
672 for (jin2 = 0; jin2 < navg; jin2++) {
673 jj = jin2 * imgcoord.nxcol;
674 for (iin = 0; iin < navg; iin++) {
675 ii = (iin + iinstart + ion) % imgcoord.nx360;
676 ix[k] = ii + jj;
677 k++;
678 }
679 }
680 }
681 }
682 else {
683 for (iout = 0; iout < Merc->header->n_columns; iout++) ix[iout] = (iout + iinstart) % imgcoord.nx360;
684 }
685
686
687 /* Now before beginning data loop, fseek if needed. */
688 if (jinstart > 0 && jinstart < (unsigned int)imgcoord.nyrow) fseek (fp, 2LL * imgcoord.nxcol * jinstart, SEEK_SET);
689
690 /* Now loop over output points, reading and handling data as needed */
691
692 n_expected = navg * imgcoord.nxcol;
693 for (jout = 0; jout < Merc->header->n_rows; jout++) {
694 jin = jinstart + navg * jout;
695 ij = gmt_M_ijp (Merc->header, jout, 0); /* Left-most index of this row */
696 if (jin < 0 || jin >= imgcoord.nyrow) { /* Outside latitude range; set row to NaNs */
697 for (iout = 0; iout < Merc->header->n_columns; iout++, ij++) Merc->data[ij] = GMT->session.f_NaN;
698 continue;
699 }
700 if ((fread (row, sizeof (int16_t), n_expected, fp) ) != n_expected) {
701 GMT_Report (API, GMT_MSG_ERROR, "Read failure at jin = %d.\n", jin);
702 gmt_M_free (GMT, ix); gmt_M_free (GMT, row);
703 fclose (fp);
704 Return (GMT_DATA_READ_ERROR);
705 }
706
707 #ifndef WORDS_BIGENDIAN
708 u2 = (uint16_t *)row;
709 for (iout = 0; iout < navg * imgcoord.nxcol; iout++)
710 u2[iout] = bswap16 (u2[iout]);
711 #endif
712
713 for (iout = 0, kstart = 0; iout < Merc->header->n_columns; iout++, ij++, kstart += navgsq) {
714 if (navg) {
715 csum = dsum = 0.0;
716 for (k = 0, kk = kstart; k < navgsq; k++, kk++) {
717 tempint = row[ix[kk]];
718 if (Ctrl->T.value && abs (tempint) % 2 != 0) {
719 csum += 1.0;
720 tempint--;
721 }
722 dsum += (double) tempint;
723 }
724 csum *= rnavgsq;
725 dsum *= rnavgsq;
726 }
727 else {
728 tempint = row[ix[iout]];
729 if (Ctrl->T.value && abs (tempint) %2 != 0) {
730 csum = 1.0;
731 tempint--;
732 }
733 else
734 csum = 0.0;
735 dsum = (double) tempint;
736 }
737
738 if (Ctrl->S.active) dsum *= Ctrl->S.value;
739
740 switch (Ctrl->T.value) {
741 case 0:
742 case 1:
743 Merc->data[ij] = (gmt_grdfloat) dsum;
744 break;
745 case 2:
746 Merc->data[ij] = (gmt_grdfloat)((csum >= 0.5) ? dsum : GMT->session.f_NaN);
747 break;
748 case 3:
749 Merc->data[ij] = (gmt_grdfloat)csum;
750 break;
751 }
752
753
754 if (Ctrl->T.value != 2 || csum >= 0.5) {
755 if (Merc->header->z_min > Merc->data[ij]) Merc->header->z_min = Merc->data[ij];
756 if (Merc->header->z_max < Merc->data[ij]) Merc->header->z_max = Merc->data[ij];
757 }
758 }
759 }
760 fclose (fp);
761
762 gmt_M_free (GMT, ix);
763 gmt_M_free (GMT, row);
764
765 if (!Ctrl->E.active) { /* Update R history with exact region found above */
766 /* Because the Mercator grid's equidistant y-nodes are not equidistant when converted back to geographic coordinates,
767 * the corresponding south and north coordinates for the outside of the pixels will generally not match the given
768 * south and north boundaries in the specified -R setting. Because this is only an issue in this program we replace
769 * the given -R with the exact -R in the gmt.history file so that any subsequent module call using -R shorthand will
770 * get the actual (exect) region and not the requested (initial, approximate) region. */
771 int id = gmt_get_option_id (0, "R"); /* The -R[P] item in the history array */
772 if (GMT->current.setting.run_mode == GMT_MODERN) id++; /* Instead pick the -RG item under modern mode since this is not a plotting tool */
773 gmt_M_str_free (GMT->init.history[id]); /* Free the previous history string set during parsing */
774 GMT->init.history[id] = strdup (exact_R); /* Replace it with the exact region */
775 }
776
777 /* We now have the Mercator grid in Grid. */
778
779 GMT_Report (API, GMT_MSG_INFORMATION, "Created %d by %d Mercatorized grid file. Min, Max values are %.8g %.8g\n", Merc->header->n_columns, Merc->header->n_rows, Merc->header->z_min, Merc->header->z_max);
780 if (Ctrl->M.active) { /* Write out the Mercator grid and return, no projection needed */
781 gmt_set_pad (GMT, API->pad); /* Reset to session default pad before output */
782 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Merc)) Return (API->error);
783 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Merc) != GMT_NOERROR) {
784 Return (API->error);
785 }
786 Return (GMT_NOERROR);
787 }
788
789 /* Here we need to reproject the data, and Merc becomes a temporary grid */
790
791 GMT_Report (API, GMT_MSG_INFORMATION, "Undo the implicit spherical Mercator -Jm1i projection.\n");
792 /* Preparing source and destination for GMT_grdproject */
793 /* a. Register the Mercator grid to be the source read by GMT_grdproject by passing a pointer */
794 GMT_Report (API, GMT_MSG_DEBUG, "Open Mercator Grid as grdproject virtual input\n");
795 if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_IN|GMT_IS_REFERENCE, Merc, input) != GMT_NOERROR) {
796 Return (API->error);
797 }
798 /* b. If -E: Register a grid struct Geo to be the destination allocated and written to by GMT_grdproject, else write to -G<file> */
799 if (Ctrl->E.active) { /* Since we will resample again, register a memory location for the result */
800 GMT_Report (API, GMT_MSG_DEBUG, "Register memory Grid container as grdproject output\n");
801 if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT|GMT_IS_REFERENCE, NULL, output) != GMT_NOERROR) {
802 Return (API->error);
803 }
804 }
805 else /* The output here is the final result */
806 strncpy (output, Ctrl->G.file, PATH_MAX-1);
807 sprintf (cmd, "-R%g/%g/%g/%g -Jm1i -I %s -G%s --PROJ_ELLIPSOID=Sphere --PROJ_LENGTH_UNIT=inch --GMT_HISTORY=readonly", west, east, south2, north2, input, output);
808 GMT_Report (API, GMT_MSG_DEBUG, "Calling grdproject %s.\n", cmd);
809 if (GMT_Call_Module (API, "grdproject", GMT_MODULE_CMD, cmd)!= GMT_NOERROR) { /* Inverse project the grid or fail */
810 Return (API->error);
811 }
812 /* Close the virtual file */
813 if (GMT_Close_VirtualFile (API, input) != GMT_NOERROR) {
814 Return (API->error);
815 }
816 /* Destroy the memory as well */
817 if (GMT_Destroy_Data (API, &Merc) != GMT_NOERROR) {
818 Return (API->error);
819 }
820 if (Ctrl->E.active) { /* Resample again using the given -R and the dx/dy in even minutes */
821 /* Preparing source and destination for GMT_grdsample */
822 /* a. Register the Geographic grid returned by GMT_grdproject to be the source read by GMT_grdsample by passing a pointer */
823 struct GMT_GRID *Geo = NULL;
824 GMT_Report (API, GMT_MSG_DEBUG, "Read Geo Grid container as grdproject virtual output\n");
825 if ((Geo = GMT_Read_VirtualFile (API, output)) == NULL) {
826 Return (API->error);
827 }
828 /* Close the virtual file */
829 if (GMT_Close_VirtualFile (API, output) != GMT_NOERROR) {
830 Return (API->error);
831 }
832 strcpy (Geo->header->title, "Data from Altimetry");
833 strncpy (Geo->header->z_units, z_units, GMT_GRID_UNIT_LEN80-1);
834 sprintf (Geo->header->x_units, "longitude [degrees_east]");
835 sprintf (Geo->header->y_units, "latitude [degrees_north]");
836 GMT_Report (API, GMT_MSG_DEBUG, "Open Geo Grid container as grdsample virtual input\n");
837 if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_IN|GMT_IS_REFERENCE, Geo, input) != GMT_NOERROR) {
838 Return (API->error);
839 }
840 sprintf (cmd, "-R%g/%g/%g/%g -I%gm %s -G%s -fg --GMT_HISTORY=readonly", west, east, south, north, Ctrl->I.value, input, Ctrl->G.file);
841 GMT_Report (API, GMT_MSG_INFORMATION, "Calling grdsample with args %s\n", cmd);
842 if (GMT_Call_Module (API, "grdsample", GMT_MODULE_CMD, cmd) != GMT_NOERROR) { /* Resample the grid or fail */
843 Return (API->error);
844 }
845 /* Close the virtual file */
846 if (GMT_Close_VirtualFile (API, input) != GMT_NOERROR) {
847 Return (API->error);
848 }
849 }
850
851 Return (GMT_NOERROR);
852 }
853