1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * Author: Paul Wessel
19 * Date: 8-NOV-2017
20 * Version: 6 API
21 *
22 * Brief synopsis: grd2kml reads a single grid and makes a Google Earth
23 * image quadtree. Optionally, supply an intensity grid (or auto-derive it)
24 * and a CPT (or use default table), and request a contour overlay.
25 * If contours are not requested we can write the PNG tiles directly without
26 * going via a PostScript plot.
27 */
28
29 #include "gmt_dev.h"
30
31 #define THIS_MODULE_CLASSIC_NAME "grd2kml"
32 #define THIS_MODULE_MODERN_NAME "grd2kml"
33 #define THIS_MODULE_LIB "core"
34 #define THIS_MODULE_PURPOSE "Create KML image quadtree from single grid"
35 #define THIS_MODULE_KEYS "<G{,CC(,IG(,WD("
36 #define THIS_MODULE_NEEDS ""
37 #define THIS_MODULE_OPTIONS "-Vfn"
38
39 /* Note: If -n is given here it is automatically set in any module called below, such as grdimage */
40
41 struct GRD2KML_CTRL {
42 struct GRD2KM_In {
43 bool active;
44 char *file;
45 } In;
46 struct GRD2KML_A { /* -Ag|s|a[<altitude>] [g] */
47 bool active;
48 int mode;
49 double altitude;
50 } A;
51 struct GRD2KML_C { /* -C<cpt> or -C<color1>,<color2>[,<color3>,...][+i<dz>] */
52 bool active;
53 double dz;
54 char *file;
55 } C;
56 struct GRD2KML_D { /* -D[+s][+d] [DEBUG ONLY, NOT DOCUMENTED] */
57 bool active;
58 bool single;
59 bool dump;
60 } D;
61 struct GRD2KML_E { /* -E<url> */
62 bool active;
63 char *url;
64 } E;
65 struct GRD2KML_F { /* -F<filter> */
66 bool active;
67 char filter;
68 } F;
69 struct GRD2KML_H { /* -H<scale> */
70 bool active;
71 int factor;
72 } H;
73 struct GRD2KML_N { /* -N<prefix> */
74 bool active;
75 char *prefix;
76 } N;
77 struct GRD2KML_I { /* -I[<intensfile>|<value>|<modifiers>] */
78 bool active;
79 bool constant;
80 bool derive;
81 double value;
82 char *azimuth; /* Default azimuth(s) for shading */
83 char *file;
84 char *method; /* Default scaling method */
85 } I;
86 struct GRD2KML_L { /* -L<size> */
87 bool active;
88 unsigned int size;
89 } L;
90 struct GRD2KML_S { /* -S[n] */
91 bool active;
92 unsigned int extra;
93 } S;
94 struct GRD2KML_T { /* -T<title> */
95 bool active;
96 char *title;
97 } T;
98 struct GRD2KML_W { /* -W<contour_file> */
99 bool active;
100 char *file;
101 double scale; /* Scaling of pen width modifier */
102 double cutoff; /* Ignore contours whose pen is < this width in points */
103 } W;
104 };
105
106 /* Structure used to keep track of which tile and its 4 possible underlings */
107 struct GMT_QUADTREE {
108 unsigned int level, q, type;
109 unsigned int row, col;
110 char tag[16];
111 char *region;
112 double wesn[4];
113 struct GMT_QUADTREE *next[4];
114 };
115
New_Ctrl(struct GMT_CTRL * GMT)116 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
117 struct GRD2KML_CTRL *C;
118
119 C = gmt_M_memory (GMT, NULL, 1, struct GRD2KML_CTRL);
120
121 /* Initialize values whose defaults are not 0/false/NULL */
122 C->F.filter = 'g';
123 C->I.method = strdup ("t1"); /* Default normalization for shading when -I is used */
124 C->L.size = 512; /* Default tile size unless global grids [360] */
125 C->W.scale = M_SQRT2;
126 C->W.cutoff = 0.1;
127 return (C);
128 }
129
Free_Ctrl(struct GMT_CTRL * GMT,struct GRD2KML_CTRL * C)130 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRD2KML_CTRL *C) { /* Deallocate control structure */
131 if (!C) return;
132 gmt_M_str_free (C->In.file);
133 gmt_M_str_free (C->C.file);
134 gmt_M_str_free (C->N.prefix);
135 gmt_M_str_free (C->I.file);
136 gmt_M_str_free (C->I.azimuth);
137 gmt_M_str_free (C->I.method);
138 gmt_M_str_free (C->T.title);
139 gmt_M_str_free (C->W.file);
140 gmt_M_free (GMT, C);
141 }
142
usage(struct GMTAPI_CTRL * API,int level)143 static int usage (struct GMTAPI_CTRL *API, int level) {
144 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
145 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
146 GMT_Usage (API, 0, "usage: %s %s -N<name> [-Aa|g|s[<altitude>]] [-C<cpt>] [-E<url>] [-F<filter>] "
147 "[-H<scale>] [-I[<intensgrid>|<value>|<modifiers>]] [-L<size>] [-S[<extra>]] [-T<title>] [%s] "
148 "[-W<contfile>|<pen>[+s<scl>/<limit>]] [%s] [%s] [%s]\n", name, GMT_INGRID, GMT_V_OPT, GMT_f_OPT, GMT_n_OPT, GMT_PAR_OPT);
149
150 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
151
152 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
153 GMT_Usage (API, 1, "\n<grid>");
154 gmt_ingrid_syntax (API, 0, "Name of grid to be plotted");
155 GMT_Usage (API, -2, "Note: Grid z-values are in user units and will be "
156 "converted to colors via the CPT [%s].", API->GMT->current.setting.cpt);
157 GMT_Usage (API, 1, "\n-N<name>");
158 GMT_Usage (API, -2, "Set file name prefix for image directory and KML file. If the directory "
159 "already exists we will overwrite the files.");
160 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
161 GMT_Usage (API, 1, "\n-Aa|g|s[<altitude>]");
162 GMT_Usage (API, -2, "Altitude mode of the image layer, choose among three modes:");
163 GMT_Usage (API, 3, "a: Absolute altitude.");
164 GMT_Usage (API, 3, "g: Altitude relative to sea surface or ground.");
165 GMT_Usage (API, 3, "s: Altitude relative to sea floor or ground.");
166 GMT_Usage (API, -2, "Optionally, append fixed <altitude> [g0: Clamped to sea surface or ground].");
167 GMT_Usage (API, 1, "\n-C<cpt>");
168 GMT_Usage (API, -2, "Color palette file to convert grid values to colors. Optionally, instead give name of a master cpt "
169 "to automatically assign continuous colors over the data range [%s]; if so, "
170 "optionally append +i<dz> to quantize the range [the exact grid range]. "
171 "Another option is to specify -C<color1>,<color2>[,<color3>,...] to build a "
172 "linear continuous cpt from those colors automatically.", API->GMT->current.setting.cpt);
173 GMT_Usage (API, 1, "\n-E<url>");
174 GMT_Usage (API, -2, "To store all files remotely, give leading URL [local files only].");
175 GMT_Usage (API, 1, "\n-F<filter>");
176 GMT_Usage (API, -2, "Specify filter type used for down-sampling (width is set automatically):");
177 GMT_Usage (API, 3, "b: Boxcar - simple averaging of all points inside filter domain.");
178 GMT_Usage (API, 3, "c: Cosine arch - weighted averaging with cosine arc weights.");
179 GMT_Usage (API, 3, "g: Gaussian - weighted averaging with Gaussian weights [Default].");
180 GMT_Usage (API, 3, "m: Median - median (50%% quantile) value of all points.");
181 GMT_Usage (API, 1, "\n-H<scale>");
182 GMT_Usage (API, -2, "Do sub-pixel smoothing using factor <scale> [no sub-pixel smoothing]. Ignored if -W not set.");
183 GMT_Usage (API, 1, "\n-I[<intensgrid>|<value>|<modifiers>]");
184 GMT_Usage (API, -2, "Apply directional illumination. Append name of intensity grid file. "
185 " For a constant intensity (i.e., change the ambient light), append a single value. "
186 "To derive intensities from <grid> instead, use -I+d to accept the default values (see grdgradient for details) or be specific:");
187 GMT_Usage (API, 3, "+a Append <azimuth> of illumination [-45]");
188 GMT_Usage (API, 3, "+n Append <method> pf intensity calculation [t1]");
189 GMT_Usage (API, 1, "\n-L<size>");
190 GMT_Usage (API, -2, "Set tile size as a power of 2 [512; for global grids, we instead select 360].n");
191 GMT_Usage (API, 1, "\n-S[<extra>]");
192 GMT_Usage (API, -2, "Add extra interpolated levels [no extra layers].");
193 GMT_Usage (API, 1, "\n-T<title>");
194 GMT_Usage (API, -2, "Set title (document description) for the top-level KML.");
195 GMT_Option (API, "V");
196 GMT_Usage (API, 1, "\n-W<contfile>|<pen>[+s<scl>/<limit>]");
197 GMT_Usage (API, -2, "Give file with select contours and pens to overlay contours [no contours]. "
198 "If no file is given we assume it is a <pen> and to use the contours implied by the CPT file. "
199 "Pen widths apply at final tile resolution and are reduced by <scl> [1.1412] for each lower level. "
200 "If a contour's scaled width is less than <limit> [0.1] it will not be drawn.");
201 GMT_Option (API, "f,n,.");
202
203 return (GMT_MODULE_USAGE);
204 }
205
parse(struct GMT_CTRL * GMT,struct GRD2KML_CTRL * Ctrl,struct GMT_OPTION * options)206 static int parse (struct GMT_CTRL *GMT, struct GRD2KML_CTRL *Ctrl, struct GMT_OPTION *options) {
207
208 /* This parses the options provided to grdcut and sets parameters in CTRL.
209 * Any GMT common options will override values set previously by other commands.
210 * It also replaces any file names specified as input or output with the data ID
211 * returned when registering these sources/destinations with the API.
212 */
213
214 unsigned int n_errors = 0, n_files = 0;
215 char *c = NULL;
216 struct GMT_OPTION *opt = NULL;
217 struct GMTAPI_CTRL *API = GMT->parent;
218
219 for (opt = options; opt; opt = opt->next) { /* Process all the options given */
220
221 switch (opt->option) {
222 /* Common parameters */
223
224 case '<': /* Input files */
225 if (n_files++ > 0) {n_errors++; continue; }
226 Ctrl->In.active = true;
227 if (opt->arg[0]) Ctrl->In.file = strdup (opt->arg);
228 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file))) n_errors++;
229 break;
230
231 /* Processes program-specific parameters */
232
233 case 'A': /* Altitude mode */
234 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
235 Ctrl->A.active = true;
236 switch (opt->arg[0]) {
237 case 'g': Ctrl->A.mode = 0; break;
238 case 's': Ctrl->A.mode = 1; break;
239 case 'a': Ctrl->A.mode = 2; break;
240 default:
241 GMT_Report (API, GMT_MSG_ERROR, "Option -A: Append either a(bsolute), g(round) or s(eafloor) with optional <altitude>.\n");
242 n_errors++;
243 break;
244 }
245 if (opt->arg[1]) Ctrl->A.altitude = atof (&opt->arg[1]);
246 break;
247 case 'C': /* CPT */
248 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
249 Ctrl->C.active = true;
250 gmt_M_str_free (Ctrl->C.file);
251 if (opt->arg[0]) Ctrl->C.file = strdup (opt->arg);
252 gmt_cpt_interval_modifier (GMT, &(Ctrl->C.file), &(Ctrl->C.dz));
253 break;
254 case 'D': /* Debug options - may fade away when happy with the performance */
255 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
256 Ctrl->D.active = true;
257 if (strstr (opt->arg, "+s")) Ctrl->D.single = true; /* Write all files in a single directory instead of one directory per level */
258 if (strstr (opt->arg, "+d")) Ctrl->D.dump = true; /* Dump quadtree information to stdout */
259 break;
260 case 'E': /* Remove URL for all contents but top driver kml */
261 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
262 Ctrl->E.active = true;
263 gmt_M_str_free (Ctrl->E.url);
264 Ctrl->E.url = strdup (opt->arg);
265 break;
266 case 'F': /* Select filter type */
267 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
268 Ctrl->F.active = true;
269 if (strchr ("bcgm", opt->arg[0]))
270 Ctrl->F.filter = opt->arg[0];
271 else {
272 GMT_Report (API, GMT_MSG_ERROR, "Option -F: Choose among b, c, g, m!\n");
273 n_errors++;
274 }
275 break;
276 case 'H': /* RIP at a higher dpi, then downsample in gs */
277 n_errors += gmt_M_repeated_module_option (API, Ctrl->H.active);
278 Ctrl->H.active = true;
279 Ctrl->H.factor = atoi (opt->arg);
280 break;
281 case 'I': /* Here, intensity must be a grid file since we need to filter it */
282 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
283 Ctrl->I.active = true;
284 if (!strcmp (opt->arg, "+d")) /* Gave +d only, so derive intensities from input grid using default settings */
285 Ctrl->I.derive = true;
286 else if ((c = gmt_first_modifier (GMT, opt->arg, "an"))) { /* Want to control how grdgradient is run */
287 unsigned int pos = 0;
288 char p[GMT_BUFSIZ] = {""};
289 Ctrl->I.derive = true;
290 while (gmt_getmodopt (GMT, 'I', c, "an", &pos, p, &n_errors) && n_errors == 0) {
291 switch (p[0]) {
292 case 'a': gmt_M_str_free (Ctrl->I.azimuth); Ctrl->I.azimuth = strdup (&p[1]); break;
293 case 'n': gmt_M_str_free (Ctrl->I.method); Ctrl->I.method = strdup (&p[1]); break;
294 default: break; /* These are caught in gmt_getmodopt so break is just for Coverity */
295 }
296 }
297 }
298 else if (!opt->arg[0] || strstr (opt->arg, "+")) /* No argument or just +, so derive intensities from input grid using default settings */
299 Ctrl->I.derive = true;
300 else if (!gmt_access (GMT, opt->arg, R_OK)) /* Got a file */
301 Ctrl->I.file = strdup (opt->arg);
302 else if (gmt_M_file_is_remote (opt->arg)) /* Got a remote file */
303 Ctrl->I.file = strdup (opt->arg);
304 else if (opt->arg[0] && !gmt_not_numeric (GMT, opt->arg)) { /* Looks like a constant value */
305 Ctrl->I.value = atof (opt->arg);
306 Ctrl->I.constant = true;
307 }
308 else {
309 GMT_Report (API, GMT_MSG_ERROR, "Option -I: Requires a valid grid file or a constant\n");
310 n_errors++;
311 }
312 break;
313 case 'L': /* Tiles sizes */
314 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
315 Ctrl->L.active = true;
316 Ctrl->L.size = atoi (opt->arg);
317 break;
318 case 'N': /* File name prefix */
319 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
320 Ctrl->N.active = true;
321 gmt_M_str_free (Ctrl->N.prefix);
322 Ctrl->N.prefix = strdup (opt->arg);
323 break;
324 case 'Q': /* Deprecated colormasking option */
325 GMT_Report (API, GMT_MSG_ERROR, "Option -Q is deprecated as transparency is automatically detected\n");
326 break;
327 case 'S': /* Extra levels */
328 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
329 Ctrl->S.active = true;
330 Ctrl->S.extra = (opt->arg[0]) ? atoi (opt->arg) : 1;
331 break;
332 case 'T': /* Title */
333 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
334 Ctrl->T.active = true;
335 if (opt->arg[0]) Ctrl->T.title = strdup (opt->arg);
336 break;
337 case 'W': /* Contours and pens */
338 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
339 Ctrl->W.active = true;
340 if ((c = strstr (opt->arg, "+s"))) { /* Gave +s<scl> modifier to scale pen widths given via -C and optional cutoff pen width */
341 sscanf (&c[2], "%lf/%lg", &Ctrl->W.scale, &Ctrl->W.cutoff);
342 c[0] = '\0'; /* Temporarily chop off the modifier */
343 }
344 if (opt->arg[0]) {
345 gmt_M_str_free (Ctrl->W.file);
346 Ctrl->W.file = strdup (opt->arg);
347 }
348 if (c) c[0] = '+'; /* Restore */
349
350 break;
351
352 default: /* Report bad options */
353 n_errors += gmt_default_error (GMT, opt->option);
354 break;
355 }
356 }
357
358 n_errors += gmt_M_check_condition (GMT, n_files != 1, "Must specify a single grid file\n");
359 n_errors += gmt_M_check_condition (GMT, Ctrl->In.file == NULL, "Must specify a single grid file\n");
360 n_errors += gmt_M_check_condition (GMT, Ctrl->N.prefix == NULL, "Option -N: Must specify a prefix for naming usage.\n");
361 n_errors += gmt_M_check_condition (GMT, Ctrl->H.active && Ctrl->H.factor <= 1, "Option -H: Must specify an integer factor > 1.\n");
362 n_errors += gmt_M_check_condition (GMT, Ctrl->E.active && Ctrl->E.url == NULL, "Option -E: Must specify an URL.\n");
363 n_errors += gmt_M_check_condition (GMT, Ctrl->I.active && !Ctrl->I.constant && !Ctrl->I.file && !Ctrl->I.derive,
364 "Option -I: Must specify intensity file, value, or modifiers\n");
365
366 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
367 }
368
grd2kml_find_quad_above(struct GMT_QUADTREE ** Q,unsigned int n,unsigned int row,unsigned int col,unsigned int level)369 GMT_LOCAL int grd2kml_find_quad_above (struct GMT_QUADTREE **Q, unsigned int n, unsigned int row, unsigned int col, unsigned int level) {
370 /* Finds the quad entry that matches the row, col, level args */
371 unsigned int k;
372 for (k = 0; k < n; k++) {
373 if (Q[k]->level != level) continue;
374 if (Q[k]->row != row) continue;
375 if (Q[k]->col != col) continue;
376 return (int)k;
377 }
378 return -1; /* Very bad */
379 }
380
grd2kml_set_dirpath(bool single,char * url,char * prefix,unsigned int level,int dir,char * string)381 GMT_LOCAL void grd2kml_set_dirpath (bool single, char *url, char *prefix, unsigned int level, int dir, char *string) {
382 if (single) { /* Write everything into the prefix dir */
383 if (url && level == 0) /* Set the leading URL for zero-level first kml */
384 sprintf (string, "%s/%s/L%2.2d", url, prefix, level);
385 else /* Everything below is in same folder */
386 sprintf (string, "L%2.2d", level);
387 }
388 else { /* Write to separate level directories */
389 if (url && level == 0) /* Set the leading URL for zero-level first kml */
390 sprintf (string, "%s/%s/%2.2d/", url, prefix, level);
391 else if (dir == -1) /* Need to refer to another directory at same level as this one */
392 sprintf (string, "../%2.2d/", level);
393 else if (dir == 0) /* At current dir */
394 string[0] = '\0';
395 else /* Down in a dir */
396 sprintf (string, "%2.2d/", level);
397 }
398 }
399
grd2kml_halve_dimensions(double * inc,double * step)400 GMT_LOCAL void grd2kml_halve_dimensions (double *inc, double *step) {
401 /* We MUST use division by two since it is a quadtree */
402 *step /= 2;
403 *inc /= 2;
404 }
405
grd2kml_max_level(struct GMT_CTRL * GMT,bool global,struct GMT_GRID_HEADER * H,unsigned int size,unsigned int extra)406 GMT_LOCAL unsigned int grd2kml_max_level (struct GMT_CTRL *GMT, bool global, struct GMT_GRID_HEADER *H, unsigned int size, unsigned int extra) {
407 unsigned int level = 0;
408 if (global) {
409 unsigned int n = 1, f = 1, go = 1;
410 double inc = 1.0, step = 360, range;
411 GMT_Report (GMT->parent, GMT_MSG_NOTICE, "Level = %2.2d tile size = %gd grid inc = %gm n_tiles = %d\n", level, step, 60*inc, n);
412 f = 2;
413 do {
414 step /= f; inc /= f;
415 n *= f;
416 range = n * H->inc[GMT_X] * 360;
417 if (range >= (360.0-GMT_CONV6_LIMIT))
418 go = 0;
419 level++;
420 GMT_Report (GMT->parent, GMT_MSG_NOTICE, "Level = %2.2d tile size = %gd grid inc = %gm n_tiles = %d\n", level, step, 60*inc, n);
421 } while (go);
422 while (extra) { /* Add the extra levels */
423 step /= f; inc /= f;
424 n *= f;
425 level++;
426 GMT_Report (GMT->parent, GMT_MSG_NOTICE, "Level = %2.2d tile size = %gd grid inc = %gm n_tiles = %d\n", level, step, 60*inc, n);
427 extra--;
428 }
429 }
430 else {
431 unsigned int mx, my;
432 mx = urint (ceil ((double)H->n_columns / (double)size)) * size; /* Nearest image size in multiples of tile size */
433 my = urint (ceil ((double)H->n_rows / (double)size)) * size;
434 level = urint (ceil (log2 (MAX (mx, my) / (double)size))) + extra; /* Number of levels in the quadtree plus optional extra levels */
435 }
436 return level;
437 }
438
grd2kml_assert_tile_size(struct GMT_CTRL * GMT,bool global,bool active,unsigned int * size)439 GMT_LOCAL void grd2kml_assert_tile_size (struct GMT_CTRL *GMT, bool global, bool active, unsigned int *size) {
440 /* For global grids we select tile size of 360 given the 360-degree range of the file, else we use radix-2.
441 * We only change size if -L was not given, else we warn if the size is not ideal. */
442
443 if (global) { /* 360 degree range in longitude and <= 180 degree range in latitude */
444 if (active && *size != 360)
445 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Option -L: For global grids the tile size should ideally be 360; your size %d is not.\n", *size);
446 else if (!active) {
447 *size = 360;
448 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Option -L: For global grids we select a tile size of %d\n", *size);
449 }
450 }
451 else { /* A smaller region, so strictly radix 2 size */
452 if (active && fabs (log2 ((double)(*size)) - irint (log2 ((double)(*size)))) > GMT_CONV8_LIMIT)
453 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Option -L: Should ideally be radix 2; your size %d is not.\n", *size);
454 }
455 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Final tile size is %d.\n", *size);
456 }
457
grd2kml_coarsen_grid(struct GMT_CTRL * GMT,unsigned int level,char filter,unsigned int registration,double orig_inc,double inc,char * DataGrid,char * Zgrid,char * filt_report)458 int grd2kml_coarsen_grid (struct GMT_CTRL *GMT, unsigned int level, char filter, unsigned int registration, double orig_inc, double inc, char *DataGrid, char *Zgrid, char *filt_report) {
459 char s_int[GMT_LEN32] = {""}, fwidth[GMT_LEN32] = {""};
460 char cmd[GMT_LEN256] = {""};
461 int error;
462 double f;
463 /* We will filter the grid to make a coarser version that exactly matches the Ctrl->L.size expectations.
464 * Note: If the desired increment is smaller than the actual, we must sample the grid instead.
465 * Note: If gridline-registered we scale filter width by sqrt(2) but rounded up a bit to make sure
466 * it extends to the 4 corner nodes. */
467 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Level %d: Low-pass-filtering the grid(s)\n", level);
468 f = (registration == GMT_GRID_NODE_REG) ? 1e-4*ceil (1e4*M_SQRT2) : 1.0; /* Round up a bit */
469 sprintf (fwidth, "%.8g", f * inc);
470 sprintf (s_int, "%.16g", inc);
471 if (inc < orig_inc) { /* Resample original instead of filtering it down */
472 sprintf (filt_report, " [Resampled with -I%s]", s_int);
473 sprintf (cmd, "%s -I%s -rp -G%s", DataGrid, s_int, Zgrid);
474 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Running grdsample : %s\n", cmd);
475 if ((error = GMT_Call_Module (GMT->parent, "grdsample", GMT_MODULE_CMD, cmd)) != GMT_NOERROR)
476 return (GMT_RUNTIME_ERROR);
477 }
478 else { /* Filter the grid */
479 unsigned int k = 0;
480 char *kind[4] = {"Boxcar", "Cosine-taper", "Gaussian", "Median"};
481 switch (filter) {
482 case 'b': k = 0; break;
483 case 'c': k = 1; break;
484 case 'g': k = 2; break;
485 case 'm': k = 3; break;
486 }
487 sprintf (filt_report, " [%s filtered with -F%c%s -I%s]", kind[k], filter, fwidth, s_int);
488 sprintf (cmd, "%s -D0 -F%c%s -I%s -rp -G%s", DataGrid, filter, fwidth, s_int, Zgrid);
489 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Running grdfilter : %s\n", cmd);
490 if ((error = GMT_Call_Module (GMT->parent, "grdfilter", GMT_MODULE_CMD, cmd)) != GMT_NOERROR)
491 return (GMT_RUNTIME_ERROR);
492 }
493 return (GMT_NOERROR);
494 }
495
496 #define bailout(code) {gmt_M_free_options (mode); return (code);}
497 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
498
GMT_grd2kml(void * V_API,int mode,void * args)499 EXTERN_MSC int GMT_grd2kml (void *V_API, int mode, void *args) {
500 bool use_tile = false, z_extend = false, i_extend = false, tmp_cpt = false, global_lon, adjust = false;
501
502 int error = 0, kk, uniq, i_dir, dir, minLodPixels, maxLodPixels;
503
504 unsigned int level, max_level, dpi = 100, n = 0, k, nx, ny, row, col, n_skip, quad, im_type;
505 unsigned int n_NaN, n_NM, n_alloc = GMT_CHUNK, n_bummer = 0, n_tiles = 0, n_img[2] = {0,0};
506 unsigned int registration;
507
508 uint64_t node;
509
510 double factor = 2.0, dim, step, wesn[4], ext_wesn[4], inc;
511
512 char cmd[GMT_BUFSIZ] = {""}, level_dir[PATH_MAX] = {""}, Zgrid[PATH_MAX] = {""}, Igrid[PATH_MAX] = {""};
513 char W[GMT_LEN16] = {""}, E[GMT_LEN16] = {""}, S[GMT_LEN16] = {""}, N[GMT_LEN16] = {""}, file[PATH_MAX] = {""};
514 char DataGrid[PATH_MAX] = {""}, IntensGrid[PATH_MAX] = {""}, path[PATH_MAX] = {""}, filt_report[GMT_LEN128] = {""};
515 char region[GMT_LEN128] = {""}, ps_cmd[GMT_LEN128] = {""}, contour_file[GMT_VF_LEN] = {""}, K[4] = {""};
516 char box[GMT_LEN32] = {""}, grdimage[GMT_LEN256] = {""}, grdcontour[GMT_LEN256] = {""}, scalepen_arg[GMT_LEN32] = {""};
517
518 static char *kml_xmlns = "<kml xmlns=\"http://www.opengis.net/kml/2.2\" xmlns:gx=\"http://www.google.com/kml/ext/2.2\" xmlns:kml=\"http://www.opengis.net/kml/2.2\" xmlns:atom=\"http://www.w3.org/2005/Atom\">";
519 static char *alt_mode[3] = {"relativeToGround", "relativeToSeaFloor", "absolute"};
520 static char *ext[2] = {"jpg", "png"}, img_code[2] = {'j', 'G'}, *transp = " -Q";
521
522 FILE *fp = NULL;
523
524 struct GMT_DATASET *C = NULL;
525 struct GMT_QUADTREE **Q = NULL;
526 struct GRD2KML_CTRL *Ctrl = NULL;
527 struct GMT_GRID *G = NULL, *T = NULL;
528 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
529 struct GMT_OPTION *options = NULL;
530 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
531
532 /*----------------------- Standard module initialization and parsing ----------------------*/
533
534 if (API == NULL) return (GMT_NOT_A_SESSION);
535 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
536 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
537
538 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
539
540 /* Parse the command-line arguments */
541
542 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 */
543 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
544 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
545 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
546
547 /*---------------------------- This is the grd2kml main code ----------------------------*/
548
549 gmt_grd_set_datapadding (GMT, true); /* Turn on gridpadding when reading a subset */
550
551 uniq = (int)getpid(); /* Unique number for temporary files */
552
553 /* Read grid header only to determine dimensions and required levels for the Pyramid */
554 if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) {
555 Return (API->error);
556 }
557 if (!gmt_M_is_geographic (GMT, GMT_IN)) {
558 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Grid must be geographic (lon, lat)\n");
559 Return (GMT_RUNTIME_ERROR);
560 }
561 if (!gmt_M_grd_equal_xy_inc (GMT, G)) {
562 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Grid spacing must be the same in longitude and latitude!\n");
563 Return (GMT_RUNTIME_ERROR);
564 }
565
566 global_lon = gmt_M_360_range (G->header->wesn[XLO], G->header->wesn[XHI]);
567
568 if (!global_lon && G->header->registration == GMT_GRID_NODE_REG) {
569 /* It is better to have a pixel-registered grid for tiling so that we can use the highest resolution exactly at the highest level.
570 * We will convert input to a fake pixel grid by extending its domain by half the increments. No resampling takes place here. */
571 GMT_Destroy_Data (API, &G); /* Delete previous grid struct */
572 sprintf (file, "%s/grd2kml_pixeldata_tmp_%6.6d.grd", API->tmp_dir, uniq);
573 sprintf (cmd, "%s -T -G%s", Ctrl->In.file, file); /* Toggle registration */
574 GMT_Report (API, GMT_MSG_INFORMATION, "Convert data grid to pixel orientation\n");
575 if ((error = GMT_Call_Module (API, "grdedit", GMT_MODULE_CMD, cmd)) != GMT_NOERROR) {
576 Return (GMT_RUNTIME_ERROR);
577 }
578 adjust = true; /* Since may need to do the same for an intensity grid */
579 gmt_M_str_free (Ctrl->In.file);
580 Ctrl->In.file = strdup (file);
581 if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) {
582 Return (API->error);
583 }
584 }
585
586 grd2kml_assert_tile_size (GMT, global_lon, Ctrl->L.active, &Ctrl->L.size); /* Set or comment on tile size */
587
588 max_level = grd2kml_max_level (GMT, global_lon, G->header, Ctrl->L.size, Ctrl->S.extra);
589
590 nx = G->header->n_columns; ny = (global_lon) ? nx : G->header->n_rows; /* Dimensions of original grid, possibly made square for global grids */
591
592 dim = dpi * 0.0001 * Ctrl->L.size; /* Constant tile map size in inches for a fixed dpi of 100 yields PNGS of the requested dimension in -L */
593
594 /* Create the container quadtree directory first */
595 if (gmt_mkdir (Ctrl->N.prefix))
596 GMT_Report (API, GMT_MSG_INFORMATION, "Directory %s already exist - will overwrite files\n", Ctrl->N.prefix);
597
598 if (Ctrl->I.derive) { /* Auto-create single intensity grid from (possibly pixel-reregistered) data grid to ensure constant scaling */
599 sprintf (file, "%s/grd2kml_intensity_tmp_%6.6d.grd", API->tmp_dir, uniq);
600 Ctrl->I.file = strdup (file);
601 GMT_Report (API, GMT_MSG_INFORMATION, "Derive an intensity grid from data grid\n");
602 /* Prepare the grdgradient arguments using selected -A -N and the data region in effect */
603 sprintf (cmd, "%s -G%s -A%s -N%s --GMT_HISTORY=readonly", Ctrl->In.file, Ctrl->I.file, Ctrl->I.azimuth, Ctrl->I.method);
604 /* Call the grdgradient module */
605 GMT_Report (API, GMT_MSG_INFORMATION, "Calling grdgradient with args %s\n", cmd);
606 if (GMT_Call_Module (API, "grdgradient", GMT_MODULE_CMD, cmd))
607 Return (API->error);
608 }
609 else if (Ctrl->I.file) { /* Can read the intensity file and compare region etc. to data grid */
610 struct GMT_GRID *I = NULL;
611 double x_noise = GMT_CONV4_LIMIT * G->header->inc[GMT_X], y_noise = GMT_CONV4_LIMIT * G->header->inc[GMT_Y];
612 if (adjust) { /* Also make intensity grid pixel-registered */
613 sprintf (file, "%s/grd2kml_pixelintens_tmp_%6.6d.grd", API->tmp_dir, uniq);
614 sprintf (cmd, "%s -T -G%s", Ctrl->I.file, file);
615 GMT_Report (API, GMT_MSG_INFORMATION, "Convert intensity grid to pixel orientation\n");
616 if ((error = GMT_Call_Module (API, "grdedit", GMT_MODULE_CMD, cmd)) != GMT_NOERROR) {
617 Return (GMT_RUNTIME_ERROR);
618 }
619 gmt_M_str_free (Ctrl->I.file);
620 Ctrl->I.file = strdup (file);
621 }
622 if ((I = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->I.file, NULL)) == NULL) {
623 Return (API->error);
624 }
625 if (!gmt_M_grd_same_shape (GMT, G, I)) {
626 GMT_Report (API, GMT_MSG_ERROR, "Data grid and intensity grid are not of the same dimensions!\n");
627 Return (GMT_RUNTIME_ERROR);
628 }
629 if (fabs (G->header->wesn[XLO] - I->header->wesn[XLO]) > x_noise || fabs (G->header->wesn[XHI] - I->header->wesn[XHI]) > x_noise ||
630 fabs (G->header->wesn[YLO] - I->header->wesn[YLO]) > y_noise || fabs (G->header->wesn[YHI] - I->header->wesn[YHI]) > y_noise) {
631 GMT_Report (API, GMT_MSG_ERROR, "Data grid and intensity grid do not cover the same area!\n");
632 Return (GMT_RUNTIME_ERROR);
633 }
634 GMT_Destroy_Data (API, &I);
635 }
636
637 Q = gmt_M_memory (GMT, NULL, n_alloc, struct GMT_QUADTREE *);
638
639 /* Determine extended region required if using the largest multiple of original grid spacing */
640
641 registration = G->header->registration;
642 if (global_lon) { /* Make it a square 360x360 grid so the quadtree splitting works */
643 ext_wesn[XLO] = G->header->wesn[XLO];
644 ext_wesn[XHI] = G->header->wesn[XHI];
645 ext_wesn[YLO] = -180.0;
646 ext_wesn[YHI] = +180.0;
647 step = 360.0; /* Initial single tile is 360x360 degrees with increment 1 degree */
648 inc = 1.0;
649 }
650 else { /* Make a square Ctrl->size region centered on the input grid center */
651 double tile_size = Ctrl->L.size * G->header->inc[GMT_X]; /* Dimension of highest-resolution tile in degrees */
652 unsigned int width;
653
654 width = pow (2.0, max_level); /* Radix-2 number of the smallest tiles in both x and y */
655 step = width * tile_size; /* Square dimension of extended grid in degrees */
656 col = G->header->n_columns / 2; row = G->header->n_rows / 2; /* Half-way point in grid */
657 /* Extend the grid half-step away from the approximate center of the input grid */
658 ext_wesn[XLO] = gmt_M_grd_col_to_x (GMT, col, G->header);
659 if (registration == GMT_GRID_PIXEL_REG) ext_wesn[XLO] -= 0.5 * G->header->inc[GMT_X]; /* Must adjust since we don't want node location but cell boundary */
660 ext_wesn[XLO] -= step / 2;
661 ext_wesn[XHI] = ext_wesn[XLO] + step;
662 ext_wesn[YLO] = gmt_M_grd_row_to_y (GMT, row, G->header);
663 if (registration == GMT_GRID_PIXEL_REG) ext_wesn[YLO] -= 0.5 * G->header->inc[GMT_Y]; /* Must adjust since we don't want node location but cell boundary */
664 ext_wesn[YLO] -= step / 2;
665 ext_wesn[YHI] = ext_wesn[YLO] + step;
666 inc = step / Ctrl->L.size;
667 nx = ny = width * Ctrl->L.size;
668 }
669 if (ext_wesn[XLO] < G->header->wesn[XLO] || ext_wesn[XHI] > G->header->wesn[XHI] || ext_wesn[YLO] < G->header->wesn[YLO] || ext_wesn[YHI] > G->header->wesn[YHI]) {
670 /* Extend the original grid with NaNs so it is an exact multiple of largest grid stride at max level */
671 sprintf (DataGrid, "%s/grd2kml_extended_data_%6.6d.grd", API->tmp_dir, uniq);
672 sprintf (cmd, "%s -R%.16g/%.16g/%.16g/%.16g -N -G%s", Ctrl->In.file, ext_wesn[XLO], ext_wesn[XHI], ext_wesn[YLO], ext_wesn[YHI], DataGrid);
673 GMT_Report (API, GMT_MSG_INFORMATION, "Extend original data grid to multiple of largest grid spacing\n");
674 if ((error = GMT_Call_Module (API, "grdcut", GMT_MODULE_CMD, cmd)) != GMT_NOERROR) {
675 gmt_M_free (GMT, Q);
676 Return (GMT_RUNTIME_ERROR);
677 }
678 z_extend = true; /* We made a temp file we need to zap later */
679 if (Ctrl->I.active) { /* Also extend the intensity grid */
680 sprintf (IntensGrid, "%s/grd2kml_extended_intens_%6.6d.grd", API->tmp_dir, uniq);
681 sprintf (cmd, "%s -R%.16g/%.16g/%.16g/%.16g -N -G%s", Ctrl->I.file, ext_wesn[XLO], ext_wesn[XHI], ext_wesn[YLO], ext_wesn[YHI], IntensGrid);
682 GMT_Report (API, GMT_MSG_INFORMATION, "Extend intensity grid to multiple of largest grid spacing\n");
683 if ((error = GMT_Call_Module (API, "grdcut", GMT_MODULE_CMD, cmd)) != GMT_NOERROR) {
684 gmt_M_free (GMT, Q);
685 Return (GMT_RUNTIME_ERROR);
686 }
687 i_extend = true; /* We made a temp file we need to zap later */
688 }
689 if (ext_wesn[YLO] < -90.0 || ext_wesn[YHI] > 90.0) gmt_grd_set_cartesian (GMT, G->header, 2); /* We must treat the extended grid as Cartesian since exceeding latitude bounds */
690 }
691 else { /* No need to extend, use the input files as is */
692 strcpy (DataGrid, Ctrl->In.file);
693 if (Ctrl->I.active)
694 strcpy (IntensGrid, Ctrl->I.file);
695 }
696
697 if (!Ctrl->C.active || gmt_is_cpt_master (GMT, Ctrl->C.file)) { /* If no cpt given or just a master then we must compute a scaled one from the full-size grid and use it throughout */
698 char *cpt = gmt_cpt_default (API, Ctrl->C.file, Ctrl->In.file);
699 char cptfile[PATH_MAX] = {""};
700 struct GMT_PALETTE *P = NULL;
701 if ((P = gmt_get_palette (GMT, cpt, GMT_CPT_OPTIONAL, G->header->z_min, G->header->z_max, Ctrl->C.dz)) == NULL) {
702 GMT_Report (API, GMT_MSG_ERROR, "Failed to create a CPT\n");
703 Return (API->error); /* Well, that did not go well... */
704 }
705 if (cpt) gmt_M_str_free (cpt);
706 sprintf (cptfile, "%s/grd2kml_%d.cpt", API->tmp_dir, uniq);
707 if (GMT_Write_Data (API, GMT_IS_PALETTE, GMT_IS_FILE, GMT_IS_NONE, GMT_WRITE_NORMAL, NULL, cptfile, P) != GMT_NOERROR) {
708 Return (API->error);
709 }
710 Ctrl->C.active = tmp_cpt = true;
711 Ctrl->C.file = strdup (cptfile);
712 }
713
714 if (Ctrl->W.active) { /* Want to overlay contours given via file */
715 /* Contour file has records of <cvalue> [<angle>] C|A [<pen>]. With angle = 0 we want <cvalue> C <pen> records */
716 uint64_t c;
717 char line[GMT_LEN256] = {""};
718 struct GMT_DATASEGMENT *S = NULL;
719
720 if (!gmt_access (GMT, Ctrl->W.file, F_OK)) { /* Was given an actual file */
721 gmt_set_cartesian (GMT, GMT_IN);
722 if ((error = GMT_Set_Columns (API, GMT_IN, 1, GMT_COL_FIX)) != GMT_NOERROR) {
723 GMT_Report (API, GMT_MSG_ERROR, "Unable to specify number of columns (1) to read\n");
724 gmt_M_free (GMT, Q);
725 Return (GMT_RUNTIME_ERROR);
726 }
727 if ((C = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_READ_NORMAL, NULL, Ctrl->W.file, NULL)) == NULL) {
728 gmt_M_free (GMT, Q);
729 Return (GMT_RUNTIME_ERROR);
730 }
731 if (C->n_segments > 1 || C->n_records == 0 || C->table[0]->segment[0]->text == NULL) {
732 GMT_Report (API, GMT_MSG_ERROR, "Contour file has more than one segment, no records at all, or no text\n");
733 gmt_M_free (GMT, Q);
734 Return (GMT_RUNTIME_ERROR);
735 }
736 S = C->table[0]->segment[0];
737 for (c = 0; c < C->n_records; c++) { /* Must reformat the records to fit grdcontour requirements */
738 if (S->text[c] == NULL) {
739 GMT_Report (API, GMT_MSG_ERROR, "No text record found\n");
740 gmt_M_free (GMT, Q);
741 Return (GMT_RUNTIME_ERROR);
742 }
743 sprintf (line, "C %s", S->text[c]); /* Build the required record format for grdcontour */
744 gmt_M_str_free (S->text[c]); /* Free previous string */
745 S->text[c] = strdup (line); /* Update string */
746 }
747 }
748 else { /* Use contours from CPT file, with -W<pen> */
749 struct GMT_PALETTE *P = NULL;
750 uint64_t dim_c[4] = {1, 1, 0, 1};
751 if ((P = GMT_Read_Data (API, GMT_IS_PALETTE, GMT_IS_FILE, GMT_IS_NONE, GMT_READ_NORMAL, NULL, Ctrl->C.file, NULL)) == NULL) {
752 gmt_M_free (GMT, Q);
753 Return (API->error);
754 }
755 dim_c[GMT_ROW] = P->n_colors + 1; /* Number of contours implied by CPT */
756 if ((C = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_NONE, GMT_WITH_STRINGS, dim_c, NULL, NULL, 0, 0, NULL)) == NULL) {
757 gmt_M_free (GMT, Q);
758 Return (GMT_RUNTIME_ERROR);
759 }
760 S = C->table[0]->segment[0];
761 sprintf (line, "C %s", Ctrl->W.file); /* Build the required record format for grdcontour */
762 for (c = 0; c < P->n_colors; c++) { /* Do all the low boundaries */
763 S->data[GMT_X][c] = P->data[c].z_low;
764 gmt_M_str_free (S->text[c]); /* Free previous string */
765 S->text[c] = strdup (line); /* Update string */
766 }
767 S->data[GMT_X][c] = P->data[P->n_colors-1].z_high;
768 gmt_M_str_free (S->text[c]); /* Free previous string */
769 S->text[c] = strdup (line); /* Update string */
770 }
771 if (GMT_Open_VirtualFile (API, GMT_IS_DATASET, GMT_IS_NONE, GMT_IN|GMT_IS_REFERENCE, C, contour_file) != GMT_NOERROR) {
772 GMT_Report (API, GMT_MSG_ERROR, "Unable to create virtual file for contours\n");
773 gmt_M_free (GMT, Q);
774 GMT_Destroy_Data (API, &C);
775 Return (GMT_RUNTIME_ERROR);
776 }
777 strcpy (K, " -K"); /* Since now we must do a contour overlay */
778 }
779
780 /* Set up the constant parts of the grdimage command */
781 sprintf (grdimage, "-JX%3.2lfi -X0 -Y0%s -W -Ve --PS_MEDIA=%3.2lfix%3.2lfi", dim, K, dim, dim);
782 if (Ctrl->C.active) { strcat (grdimage, " -C"); strcat (grdimage, Ctrl->C.file); }
783 /* Set up the constant parts of the grdcontour command */
784 if (Ctrl->W.active) /* Overlay contours */
785 sprintf (grdcontour, "-JX%3.2lfi -O -C%s -Ve", dim, contour_file);
786
787 if (Ctrl->H.active) /* Do sub-pixel smoothing */
788 sprintf (ps_cmd, "-E100 -P -Ve -Z -H%d", Ctrl->H.factor);
789 else
790 sprintf (ps_cmd, "-E100 -P -Ve -Z");
791
792 GMT_Report (GMT->parent, GMT_MSG_NOTICE, "Extended grid size is by %d by %d, tile size is %d x %d\n", ny, nx, Ctrl->L.size, Ctrl->L.size);
793 n_NM = Ctrl->L.size * Ctrl->L.size; /* Number of nodes or pixels in a tile */
794
795 /* Loop over all the levels, starting at the top level (0) */
796 for (level = 0; level <= max_level; level++) {
797 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Level %d: Factor = %g Dim = %d x %d -> %d x %d\n",
798 level, factor, irint (factor * Ctrl->L.size), irint (factor * Ctrl->L.size), Ctrl->L.size, Ctrl->L.size);
799 /* Create the level directory */
800 if (!Ctrl->D.single) {
801 sprintf (level_dir, "%s/%2.2d", Ctrl->N.prefix, level);
802 if (gmt_mkdir (level_dir))
803 GMT_Report (API, GMT_MSG_INFORMATION, "Level directory %s already exist - overwriting files\n", level_dir);
804 }
805 if (level < max_level || registration == GMT_GRID_NODE_REG || !doubleAlmostEqual (inc, G->header->inc[GMT_X])) { /* Filter the data to match level resolution */
806 sprintf (Zgrid, "%s/grd2kml_Z_L%d_tmp_%6.6d.grd", API->tmp_dir, level, uniq);
807 if (grd2kml_coarsen_grid (GMT, level, Ctrl->F.filter, registration, G->header->inc[GMT_X], inc, DataGrid, Zgrid, filt_report)) {
808 gmt_M_free (GMT, Q);
809 if (Ctrl->W.active) GMT_Destroy_Data (API, &C);
810 Return (GMT_RUNTIME_ERROR);
811 }
812 if (Ctrl->I.active) { /* Also filter the intensity grid */
813 sprintf (Igrid, "%s/grd2kml_I_L%d_tmp_%6.6d.grd", API->tmp_dir, level, uniq);
814 if (grd2kml_coarsen_grid (GMT, level, Ctrl->F.filter, registration, G->header->inc[GMT_X], inc, IntensGrid, Igrid, filt_report)) {
815 gmt_M_free (GMT, Q);
816 if (Ctrl->W.active) GMT_Destroy_Data (API, &C);
817 Return (GMT_RUNTIME_ERROR);
818 }
819 }
820 }
821 else { /* Use as is for the highest resolution */
822 sprintf (filt_report, " [Original grid used]");
823 strcpy (Zgrid, DataGrid);
824 step = Ctrl->L.size * G->header->inc[GMT_X];
825 if (Ctrl->I.active) strcpy (Igrid, IntensGrid);
826 }
827 /* Loop over all rows at this level */
828 row = col = n_skip = 0;
829 wesn[YLO] = ext_wesn[YLO];
830 gmt_ascii_format_one (GMT, S, wesn[YLO], GMT_IS_FLOAT);
831
832 if (Ctrl->W.active) { /* Set pen width scale and overall cutoff pen size [0.1] in points */
833 double p = pow (Ctrl->W.scale, -(double)(max_level - level));
834 sprintf (scalepen_arg, " -W+s%g/%g", p, Ctrl->W.cutoff);
835 }
836 while (wesn[YLO] < (G->header->wesn[YHI]-G->header->inc[GMT_Y])) { /* Small correction to avoid issues due to round-off */
837 wesn[YHI] = wesn[YLO] + step; /* Top row may extend beyond grid and be transparent */
838 gmt_ascii_format_one (GMT, N, wesn[YHI], GMT_IS_FLOAT); /* GMT_IS_FLOAT and not GMT_IS_LAT since we may exceed 90 */
839 if (wesn[YHI] <= G->header->wesn[YLO] || wesn[YLO] >= G->header->wesn[YHI]) { /* Tile row outside data range */
840 row++; /* Onwards to next row */
841 wesn[YLO] = wesn[YHI];
842 strcpy (S, N);
843 n_skip++;
844 continue;
845 }
846 /* Loop over all columns at this level */
847 col = 0;
848 wesn[XLO] = ext_wesn[XLO];
849 gmt_ascii_format_one (GMT, W, wesn[XLO], GMT_IS_FLOAT);
850 while (wesn[XLO] < (G->header->wesn[XHI]-G->header->inc[GMT_X])) { /* Small correction to avoid issues due to round-off */
851 unsigned int trow, tcol;
852 wesn[XHI] = wesn[XLO] + step; /* So right column may extend beyond grid and be transparent */
853 gmt_ascii_format_one (GMT, E, wesn[XHI], GMT_IS_FLOAT);
854 if (wesn[XHI] <= G->header->wesn[XLO] || wesn[XLO] >= G->header->wesn[XHI]) { /* Tile outside x data range */
855 col++; /* Onwards to next column */
856 wesn[XLO] = wesn[XHI];
857 strcpy (W, E);
858 n_skip++;
859 continue;
860 }
861 /* Now we have the current tile region */
862 if ((T = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, wesn, Zgrid, NULL)) == NULL) {
863 GMT_Report (API, GMT_MSG_ERROR, "Unable to read in grid tile!\n");
864 gmt_M_free (GMT, Q);
865 if (Ctrl->W.active) GMT_Destroy_Data (API, &C);
866 Return (API->error);
867 }
868 /* Determine if we have any non-NaN data points inside this grid */
869 for (trow = n_NaN = 0; trow < T->header->n_rows; trow++) {
870 for (tcol = 0; tcol < T->header->n_columns; tcol++) {
871 node = gmt_M_ijp (T->header, trow, tcol);
872 if (gmt_M_is_fnan (T->data[node])) n_NaN++;
873 }
874 }
875 use_tile = (n_NaN < n_NM);
876 im_type = (n_NaN) ? 1 : 0;
877 if (use_tile) { /* Found data inside this tile, make plot and rasterize (if PostScript) */
878 char z_data[GMT_VF_LEN] = {""}, psfile[PATH_MAX] = {""};
879 /* Open the grid subset as a virtual file we can pass to grdimage */
880 if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_IN|GMT_IS_REFERENCE, T, z_data) == GMT_NOTSET) {
881 GMT_Report (API, GMT_MSG_ERROR, "Unable to open grid tile as virtual file!\n");
882 gmt_M_free (GMT, Q);
883 if (Ctrl->W.active) GMT_Destroy_Data (API, &C);
884 Return (API->error);
885 }
886 /* Will pass -W so grdimage will notify us if there was no valid image data imaged */
887 if (!Ctrl->W.active) {
888 /* Build the grdimage command to write an image directly via GDAL - no need to go via PostSCript */
889 char imagefile[PATH_MAX] = {""};
890 if (Ctrl->D.single)
891 sprintf (imagefile, "%s/L%2.2dR%3.3dC%3.3d.%s", Ctrl->N.prefix, level, row, col, ext[im_type]);
892 else
893 sprintf (imagefile, "%s/R%3.3dC%3.3d.%s", level_dir, row, col, ext[im_type]);
894 if (Ctrl->I.active) /* Must pass two grids */
895 sprintf (cmd, "%s %s -I%s -R%s/%s/%s/%s -A%s", grdimage, z_data, Igrid, W, E, S, N, imagefile);
896 else
897 sprintf (cmd, "%s %s -R%s/%s/%s/%s -A%s", grdimage, z_data, W, E, S, N, imagefile);
898 if (im_type) strcat (cmd, transp);
899 error = GMT_Call_Module (API, "grdimage", GMT_MODULE_CMD, cmd);
900 if (!(error == GMT_NOERROR || error == GMT_IMAGE_NO_DATA)) {
901 GMT_Report (API, GMT_MSG_ERROR, "Unable to create a direct PNG from grid\n");
902 gmt_M_free (GMT, Q);
903 Return (API->error);
904 }
905 }
906 else {
907 /* Build the grdimage command to make the PostScript plot */
908 sprintf (psfile, "%s/grd2kml_tile_tmp_%6.6d.ps", API->tmp_dir, uniq);
909 if (Ctrl->I.active) /* Must pass two grids */
910 sprintf (cmd, "%s %s -I%s -R%s/%s/%s/%s ->%s", grdimage, z_data, Igrid, W, E, S, N, psfile);
911 else
912 sprintf (cmd, "%s %s -R%s/%s/%s/%s ->%s", grdimage, z_data, W, E, S, N, psfile);
913 if (im_type) strcat (cmd, transp);
914 error = GMT_Call_Module (API, "grdimage", GMT_MODULE_CMD, cmd);
915 if (error == GMT_NOERROR && Ctrl->W.active) { /* Overlay contours */
916 sprintf (cmd, "%s %s -R%s/%s/%s/%s %s ->>%s", grdcontour, z_data, W, E, S, N, scalepen_arg, psfile);
917 GMT_Init_VirtualFile (API, 0, z_data); /* Read the same grid again */
918 GMT_Init_VirtualFile (API, 0, contour_file); /* Read the same contours again */
919 if ((error = GMT_Call_Module (API, "grdcontour", GMT_MODULE_CMD, cmd))) {
920 GMT_Report (API, GMT_MSG_ERROR, "Unable to overlay contours!\n");
921 gmt_M_free (GMT, Q);
922 GMT_Destroy_Data (API, &C);
923 Return (API->error);
924 }
925 }
926 }
927 GMT_Close_VirtualFile (API, z_data);
928 if (GMT_Destroy_Data (API, &T) != GMT_NOERROR) {
929 GMT_Report (API, GMT_MSG_ERROR, "Unable to free memory of grid tile!\n");
930 gmt_M_free (GMT, Q);
931 Return (API->error);
932 }
933 if (error == GMT_IMAGE_NO_DATA) { /* Must have found non-NaNs in the pad since the image is all NaN? */
934 GMT_Report (API, GMT_MSG_WARNING, "No image content for current tile (%d, %d, %d) - skipped\n", level, row, col);
935 gmt_remove_file (GMT, psfile);
936 n_bummer++;
937 }
938 else { /* Made a meaningful plot, time to rip. */
939 /* Create the psconvert command to convert the PS to transparent PNG */
940 sprintf (region, "%s/%s/%s/%s", W, E, S, N);
941 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Level %2.2d: Mapped tile %s\n", level, region);
942 if (Ctrl->W.active) { /* Make PS, now rasterize to PNG */
943 if (Ctrl->D.single)
944 sprintf (cmd, "%s -T%c -D%s -FL%2.2dR%3.3dC%3.3d %s", ps_cmd, img_code[im_type], Ctrl->N.prefix, level, row, col, psfile);
945 else
946 sprintf (cmd, "%s -T%c -D%s -FR%3.3dC%3.3d %s", ps_cmd, img_code[im_type], level_dir, row, col, psfile);
947 if (GMT_Call_Module (API, "psconvert", GMT_MODULE_CMD, cmd)) {
948 GMT_Report (API, GMT_MSG_ERROR, "Unable to rasterize current PNG tile!\n");
949 gmt_M_free (GMT, Q);
950 Return (API->error);
951 }
952 }
953 /* Update our list of tiles created */
954 Q[n] = gmt_M_memory (GMT, NULL, 1, struct GMT_QUADTREE);
955 Q[n]->row = row; Q[n]->col = col; Q[n]->level = level; Q[n]->type = im_type;
956 Q[n]->wesn[XLO] = wesn[XLO]; Q[n]->wesn[XHI] = wesn[XHI];
957 Q[n]->wesn[YLO] = wesn[YLO]; Q[n]->wesn[YHI] = wesn[YHI];
958 sprintf (Q[n]->tag, "L%2.2dR%3.3dC%3.3d", level, row, col);
959 Q[n]->region = strdup (region);
960 if (++n == n_alloc) { /* Extend the array */
961 n_alloc <<= 1;
962 Q = gmt_M_memory (GMT, Q, n_alloc, struct GMT_QUADTREE *);
963 }
964 n_img[im_type]++;
965 }
966 }
967 else { /* Just NaNs inside this tile */
968 GMT_Report (API, GMT_MSG_INFORMATION, "Level %2.2d: Tile %s/%s/%s/%s had no data - skipped\n", level, W, E, S, N);
969 n_skip++;
970 }
971 col++; /* Onwards to next column */
972 wesn[XLO] = wesn[XHI];
973 strcpy (W, E);
974 }
975 row++; /* Onwards to next row */
976 wesn[YLO] = wesn[YHI];
977 strcpy (S, N);
978 }
979 if (step > 1)
980 sprintf (box, "%g x %g d", step, step);
981 else
982 sprintf (box, "%g x %g m", 60*step, 60*step);
983 GMT_Report (GMT->parent, GMT_MSG_NOTICE, "Level %2.2d: Tile size: %18s Tiles: %3d by %3d = %5d %5d mapped %3d empty%s\n", level, box, row, col, row*col, row*col - n_skip, n_skip, filt_report);
984 if (level < max_level) { /* Delete the temporary filtered grid(s) */
985 gmt_remove_file (GMT, Zgrid);
986 if (Ctrl->I.active) gmt_remove_file (GMT, Igrid);
987 }
988 grd2kml_halve_dimensions (&inc, &step);
989 }
990 n_tiles = n;
991
992 if (n_bummer)
993 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Found %d tiles that passed the no-NaN test but gave a blank image (?)\n", n_bummer);
994
995 /* Clean up any temporary files */
996
997 if (z_extend && !access (DataGrid, F_OK))
998 gmt_remove_file (GMT, DataGrid);
999 if (i_extend && !access (IntensGrid, F_OK))
1000 gmt_remove_file (GMT, IntensGrid);
1001 if (Ctrl->I.derive)
1002 gmt_remove_file (GMT, Ctrl->I.file);
1003 if (tmp_cpt)
1004 gmt_remove_file (GMT, Ctrl->C.file);
1005
1006 /* Process quadtree links */
1007
1008 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Processes quadtree links for %d tiles.\n", n);
1009 Q = gmt_M_memory (GMT, Q, n, struct GMT_QUADTREE *); /* Final size */
1010 for (level = max_level; level > 0; level--) {
1011 for (k = 0; k < n; k++) {
1012 if (Q[k]->level != level) continue; /* Only deal with this level here */
1013 /* Determine the parent tile and the quad (0-3) we belong to */
1014 /* This is the parent row and col since we increase by a factor of 2 each time */
1015 row = Q[k]->row / 2; col = Q[k]->col / 2;
1016 /* The quad is given by comparing the high and low values of row, col */
1017 quad = 2 * (Q[k]->row - 2 * row) + (Q[k]->col - 2 * col);
1018 kk = grd2kml_find_quad_above (Q, n, row, col, level-1); /* kk is the parent of k */
1019 if (kk < 0) { /* THis can happen when the lower-level tile grazes one above it but there really are no data involved */
1020 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Tile %s: Unable to link tile for row = %d, col = %d at level %d to a parent (!?). Probably empty - skipped.\n", Q[k]->tag, row, col, level);
1021 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Tile %s: Region was %g/%g/%g/%g.\n", Q[k]->tag, Q[k]->wesn[XLO], Q[k]->wesn[XHI], Q[k]->wesn[YLO], Q[k]->wesn[YHI]);
1022 continue;
1023 }
1024 assert (quad < 4); /* Sanity check */
1025 Q[kk]->next[quad] = Q[k]; /* Do the linking */
1026 Q[kk]->q++; /* Count the links for this parent */
1027 }
1028 }
1029
1030 /* Create all the KML files in the quadtree with their links down the tree */
1031
1032 for (k = 0; k < n; k++) {
1033 if (Ctrl->D.dump) {
1034 printf ("%s [%s]:\t", Q[k]->tag, Q[k]->region);
1035 for (quad = 0; quad < 4; quad++) {
1036 if (Q[k]->next[quad] == NULL) continue;
1037 minLodPixels = (Q[k]->next[quad]->level == 0) ? 1 : Ctrl->L.size / 2;
1038 maxLodPixels = (Q[k]->next[quad]->level == max_level) ? -1 : 4*Ctrl->L.size;
1039 printf (" %c=%s [%d/%d]", 'A'+quad, Q[k]->next[quad]->tag, minLodPixels, maxLodPixels);
1040 }
1041 printf ("\n");
1042 }
1043 if (Q[k]->level == 0)
1044 sprintf (file, "%s/%s.kml", Ctrl->N.prefix, Ctrl->N.prefix);
1045 else if (Ctrl->D.single)
1046 sprintf (file, "%s/L%2.2dR%3.3dC%3.3d.kml", Ctrl->N.prefix, Q[k]->level, Q[k]->row, Q[k]->col);
1047 else
1048 sprintf (file, "%s/%2.2d/R%3.3dC%3.3d.kml", Ctrl->N.prefix, Q[k]->level, Q[k]->row, Q[k]->col);
1049 if ((fp = fopen (file, "w")) == NULL) {
1050 GMT_Report (API, GMT_MSG_ERROR, "Unable to create file : %s\n", file);
1051 Return (GMT_RUNTIME_ERROR);
1052 }
1053 /* First this tile's kml info and image href */
1054 minLodPixels = (Q[k]->level == 0) ? 1 : Ctrl->L.size / 2;
1055 maxLodPixels = (Q[k]->level == max_level) ? -1 : 4*Ctrl->L.size;
1056 grd2kml_set_dirpath (Ctrl->D.single, NULL, Ctrl->N.prefix, Q[k]->level, 1, path);
1057 fprintf (fp, "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n%s\n", kml_xmlns);
1058 if (Q[k]->level == 0) {
1059 char *cmd_args = GMT_Create_Cmd (API, options);
1060 fprintf (fp, " <Document>\n <name>%s</name>\n", Ctrl->N.prefix);
1061 fprintf (fp, " <!-- Produced by the Generic Mapping Tools [www.generic-mapping-tools.org] -->\n");
1062 fprintf (fp, " <!-- cmd: gmt %s %s -->\n", GMT->init.module_name, cmd_args);
1063 gmt_M_free (GMT, cmd_args);
1064 dir = 0;
1065 if (Ctrl->T.active)
1066 fprintf (fp, " <description>%s</description>\n\n", Ctrl->T.title);
1067 else
1068 fprintf (fp, " <description>GMT image quadtree representation of %s</description>\n\n", Ctrl->In.file);
1069 }
1070 else {
1071 fprintf (fp, " <Document>\n <name>%s.kml</name>\n", Q[k]->tag);
1072 fprintf (fp, " <description></description>\n\n");
1073 dir = 1;
1074 }
1075 fprintf (fp, " <Style>\n");
1076 fprintf (fp, " <ListStyle id=\"hideChildren\">\n");
1077 fprintf (fp, " <listItemType>checkHideChildren</listItemType>\n");
1078 fprintf (fp, " </ListStyle>\n");
1079 fprintf (fp, " </Style>\n");
1080 fprintf (fp, " <Region>\n <LatLonAltBox>\n");
1081 fprintf (fp, " <north>%.14g</north>\n", Q[k]->wesn[YHI]);
1082 fprintf (fp, " <south>%.14g</south>\n", Q[k]->wesn[YLO]);
1083 fprintf (fp, " <east>%.14g</east>\n", Q[k]->wesn[XHI]);
1084 fprintf (fp, " <west>%.14g</west>\n", Q[k]->wesn[XLO]);
1085 fprintf (fp, " </LatLonAltBox>\n");
1086 fprintf (fp, " <Lod>\n");
1087 fprintf (fp, " <minLodPixels>%d</minLodPixels>\n", minLodPixels);
1088 fprintf (fp, " <maxLodPixels>%d</maxLodPixels>\n", maxLodPixels);
1089 fprintf (fp, " </Lod>\n");
1090 fprintf (fp, " </Region>\n");
1091 fprintf (fp, " <GroundOverlay>\n");
1092 fprintf (fp, " <drawOrder>%d</drawOrder>\n", 10+2*Q[k]->level);
1093 grd2kml_set_dirpath (Ctrl->D.single, NULL, Ctrl->N.prefix, Q[k]->level, 1-dir, path);
1094 fprintf (fp, " <Icon>\n <href>%sR%3.3dC%3.3d.%s</href>\n </Icon>\n", path, Q[k]->row, Q[k]->col, ext[Q[k]->type]);
1095 fprintf (fp, " <LatLonBox>\n");
1096 fprintf (fp, " <north>%.14g</north>\n", Q[k]->wesn[YHI]);
1097 fprintf (fp, " <south>%.14g</south>\n", Q[k]->wesn[YLO]);
1098 fprintf (fp, " <east>%.14g</east>\n", Q[k]->wesn[XHI]);
1099 fprintf (fp, " <west>%.14g</west>\n", Q[k]->wesn[XLO]);
1100 fprintf (fp, " </LatLonBox>\n");
1101 fprintf (fp, " <altitudeMode>%s</altitudeMode><altitude>%g</altitude>\n", alt_mode[Ctrl->A.mode], Ctrl->A.altitude);
1102 fprintf (fp, " </GroundOverlay>\n");
1103 i_dir = (Q[k]->level == 0) ? 1 : -1;
1104 /* Now add up to 4 quad links */
1105 for (quad = 0; quad < 4; quad++) {
1106 if (Q[k]->next[quad] == NULL) continue;
1107 minLodPixels = (Q[k]->level == 0) ? 1 : Ctrl->L.size / 2;
1108 maxLodPixels = (Q[k]->level == max_level) ? -1 : 4*Ctrl->L.size;
1109
1110 grd2kml_set_dirpath (Ctrl->D.single, NULL, Ctrl->N.prefix, Q[k]->next[quad]->level, dir, path);
1111 fprintf (fp, "\n <NetworkLink>\n <name>%s</name>\n", Q[k]->next[quad]->tag);
1112 fprintf (fp, " <Region>\n <LatLonAltBox>\n");
1113 fprintf (fp, " <north>%.14g</north>\n", Q[k]->next[quad]->wesn[YHI]);
1114 fprintf (fp, " <south>%.14g</south>\n", Q[k]->next[quad]->wesn[YLO]);
1115 fprintf (fp, " <east>%.14g</east>\n", Q[k]->next[quad]->wesn[XHI]);
1116 fprintf (fp, " <west>%.14g</west>\n", Q[k]->next[quad]->wesn[XLO]);
1117 fprintf (fp, " </LatLonAltBox>\n");
1118 fprintf (fp, " <Lod>\n");
1119 fprintf (fp, " <minLodPixels>%d</minLodPixels>\n", minLodPixels);
1120 fprintf (fp, " <maxLodPixels>%d</maxLodPixels>\n", maxLodPixels);
1121 fprintf (fp, " </Lod>\n");
1122 fprintf (fp, " </Region>\n");
1123 grd2kml_set_dirpath (Ctrl->D.single, NULL, Ctrl->N.prefix, Q[k]->next[quad]->level, i_dir, path);
1124 fprintf (fp, " <Link>\n <href>%sR%3.3dC%3.3d.kml</href>\n", path, Q[k]->next[quad]->row, Q[k]->next[quad]->col);
1125 fprintf (fp, " <viewRefreshMode>onRegion</viewRefreshMode>\n");
1126 fprintf (fp, " </Link>\n");
1127 fprintf (fp, " </NetworkLink>\n");
1128 }
1129 fprintf (fp, " </Document>\n</kml>\n");
1130 fclose (fp);
1131
1132 gmt_M_str_free (Q[k]->region); /* Free this tile region */
1133 gmt_M_free (GMT, Q[k]); /* Free this tile information */
1134 }
1135 gmt_M_free (GMT, Q);
1136 GMT_Report (API, GMT_MSG_NOTICE, "Done: %d tiles (%d JPG and %d PNG) written to directory %s\n", n_tiles, n_img[0], n_img[1], Ctrl->N.prefix);
1137
1138 Return (GMT_NOERROR);
1139 }
1140