-
Notifications
You must be signed in to change notification settings - Fork 375
/
Copy pathgmt_grdio.c
4045 lines (3679 loc) · 186 KB
/
gmt_grdio.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*--------------------------------------------------------------------
*
* Copyright (c) 1991-2025 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
* See LICENSE.TXT file for copying and redistribution conditions.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation; version 3 or any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* Contact info: www.generic-mapping-tools.org
*--------------------------------------------------------------------*/
/*
*
* G M T _ G R D I O . C R O U T I N E S
*
* Generic routines that take care of all low-level gridfile input/output.
*
* Author: Paul Wessel
* Date: 1-JAN-2010
* Version: 5
* 64-bit Compliant: Yes
*
* NOTE: We do not support grids that have more rows or columns that can fit
* in a 32-bit integer, hence all linear dimensions (row, col, etc.) are addressed
* with 32-bit ints. However, 1-D array indices (e.g., ij = row*n_columns + col) are
* addressed with 64-bit integers.
*
* A) List of exported gmt_* functions available to modules and libraries via gmt_dev.h:
*
* gmt_adjust_loose_wesn
* gmt_change_grdreg
* gmt_change_grid_history
* gmt_copy_gridheader
* gmt_create_grid
* gmt_cube_info_syntax
* gmt_cube_pad_off
* gmt_cube_vminmax
* gmt_decode_cube_h_info
* gmt_decode_grd_h_info
* gmt_duplicate_grid
* gmt_free_grid
* gmt_free_header
* gmt_get_active_layers
* gmt_get_grd_command
* gmt_get_grd_remark
* gmt_get_grd_title
* gmt_get_grid
* gmt_grd_RI_verify
* gmt_grd_coord
* gmt_grd_detrend
* gmt_grd_domains_match
* gmt_grd_dump
* gmt_grd_flip_vertical
* gmt_grd_format_decoder
* gmt_grd_get_format
* gmt_grd_info_syntax
* gmt_grd_init
* gmt_grd_is_global
* gmt_grd_is_polar
* gmt_grd_layout
* gmt_grd_minmax
* gmt_grd_mux_demux
* gmt_grd_pad_off
* gmt_grd_pad_on
* gmt_grd_pad_status
* gmt_grd_pad_zero
* gmt_grd_prep_io
* gmt_grd_set_cartesian
* gmt_grd_set_datapadding
* gmt_grd_set_ij_inc
* gmt_grd_setregion
* gmt_grd_shift
* gmt_grd_zminmax
* gmt_grdcube_info_syntax
* gmt_grid_perimeter
* gmt_img_sanitycheck
* gmt_raster_type
* gmt_read_img
* gmt_scale_and_offset_f
* gmt_set_R_from_grd
* gmt_set_grddim
* gmt_set_grdinc
* gmt_set_outgrid
* gmt_update_grd_info
*
* B) List of exported gmtlib_* functions available to libraries via gmt_internals.h:
*
* gmtlib_create_cube
* gmtlib_duplicate_cube
* gmtlib_found_url_for_gdal
* gmtlib_free_cube
* gmtlib_free_cube_ptr
* gmtlib_free_grid_ptr
* gmtlib_get_grdtype
* gmtlib_get_matrixtype
* gmtlib_grd_data_size
* gmtlib_grd_get_units
* gmtlib_grd_real_interleave
* gmtlib_grd_set_units
* gmtlib_init_complex
* gmtlib_read_grd
* gmtlib_read_grd_info
* gmtlib_read_image
* gmtlib_read_image_info
* gmtlib_write_grd
* gmtlib_write_grd_info
* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
#include "gmt_dev.h"
#include "gmt_internals.h"
struct GRD_PAD { /* Local structure */
double wesn[4];
unsigned int pad[4];
};
/* Local functions */
/*! gmt_M_grd_get_size computes grid size including the padding, and doubles it if complex values */
GMT_LOCAL size_t gmtgrdio_grd_get_size (struct GMT_GRID_HEADER *h) {
size_t L = (((h->complex_mode & GMT_GRID_IS_COMPLEX_MASK) > 0) + 1ULL) * h->mx * h->my;
if (L % 2) L++; /* Make it an even number just in case an array is passed to an operator that expects real,imag pairs */
return (L);
}
int gmt_grd_layout (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, unsigned int complex_mode, unsigned int direction) {
/* Checks or sets the array arrangement for a complex array */
size_t needed_size; /* Space required to hold both components of a complex grid */
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
if ((complex_mode & GMT_GRID_IS_COMPLEX_MASK) == 0) return GMT_OK; /* Regular, non-complex grid, nothing special to do */
needed_size = 2ULL * ((size_t)header->mx) * ((size_t)header->my); /* For the complex array */
if (header->size < needed_size) {
GMT_Report (GMT->parent, GMT_MSG_ERROR, "Complex grid not large enough to hold both components!.\n");
return GMT_DIM_TOO_SMALL;
}
if (direction == GMT_IN) { /* About to read in a complex component; another one might have been read in earlier */
if (HH->arrangement == GMT_GRID_IS_INTERLEAVED) {
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Demultiplexing complex grid before reading can take place.\n");
gmt_grd_mux_demux (GMT, header, grid, GMT_GRID_IS_SERIAL);
}
if ((header->complex_mode & GMT_GRID_IS_COMPLEX_MASK) == GMT_GRID_IS_COMPLEX_MASK) { /* Already have both component; this will overwrite one of them */
unsigned int type = (complex_mode & GMT_GRID_IS_COMPLEX_REAL) ? 0 : 1;
char *kind[2] = {"read", "imaginary"};
GMT_Report (GMT->parent, GMT_MSG_WARNING, "Overwriting previously stored %s component in complex grid.\n", kind[type]);
}
header->complex_mode |= complex_mode; /* Update the grids complex mode */
}
else { /* About to write out a complex component */
if ((header->complex_mode & GMT_GRID_IS_COMPLEX_MASK) == 0) { /* Not a complex grid */
GMT_Report (GMT->parent, GMT_MSG_ERROR, "Asking to write out complex components from a non-complex grid.\n");
return GMT_WRONG_MATRIX_SHAPE;
}
if ((header->complex_mode & GMT_GRID_IS_COMPLEX_REAL) && (header->complex_mode & GMT_GRID_IS_COMPLEX_REAL) == 0) {
/* Programming error: Requesting to write real components when there are none */
GMT_Report (GMT->parent, GMT_MSG_ERROR, "Complex grid has no real components that can be written to file.\n");
return GMT_WRONG_MATRIX_SHAPE;
}
else if ((header->complex_mode & GMT_GRID_IS_COMPLEX_IMAG) && (header->complex_mode & GMT_GRID_IS_COMPLEX_IMAG) == 0) {
/* Programming error: Requesting to write imag components when there are none */
GMT_Report (GMT->parent, GMT_MSG_ERROR, "Complex grid has no imaginary components that can be written to file.\n");
return GMT_WRONG_MATRIX_SHAPE;
}
if (HH->arrangement == GMT_GRID_IS_INTERLEAVED) { /* Must first demultiplex the grid */
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Demultiplexing complex grid before writing can take place.\n");
gmt_grd_mux_demux (GMT, header, grid, GMT_GRID_IS_SERIAL);
}
}
/* header->arrangement might now have been changed accordingly */
return GMT_OK;
}
GMT_LOCAL void gmtgrdio_grd_parse_xy_units (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, char *file, unsigned int direction) {
/* Decode the optional +u|U<unit> and determine scales */
enum gmt_enum_units u_number;
unsigned int mode = 0;
char *c = NULL, *name = NULL, *f = NULL;
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (h);
if (gmt_M_is_geographic (GMT, direction)) return; /* Does not apply to geographic data */
name = (file) ? file : HH->name;
if ((f = gmt_strrstr (name, ".grd")) || (f = gmt_strrstr (file, ".nc")))
c = gmtlib_last_valid_file_modifier (GMT->parent, f, "Uu");
else
c = gmtlib_last_valid_file_modifier (GMT->parent, name, "Uu");
if (c == NULL) return; /* Did not find any modifier */
mode = (c[1] == 'u') ? 0 : 1;
u_number = gmtlib_get_unit_number (GMT, c[2]); /* Convert char unit to enumeration constant for this unit */
if (u_number == GMT_IS_NOUNIT) {
GMT_Report (GMT->parent, GMT_MSG_WARNING, "Grid file x/y unit specification %s was unrecognized (part of file name?) and is ignored.\n", c);
return;
}
/* Got a valid unit */
HH->xy_unit_to_meter[direction] = GMT->current.proj.m_per_unit[u_number]; /* Converts unit to meters */
if (mode) HH->xy_unit_to_meter[direction] = 1.0 / HH->xy_unit_to_meter[direction]; /* Wanted the inverse */
HH->xy_unit[direction] = u_number; /* Unit ID */
HH->xy_adjust[direction] |= 1; /* Says we have successfully parsed and readied the x/y scaling */
HH->xy_mode[direction] = mode;
c[0] = '\0'; /* Chop off the unit specification from the file name */
}
GMT_LOCAL void gmtgrdio_grd_xy_scale (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, unsigned int direction) {
unsigned int k;
/* Apply the scaling of wesn,inc as given by the header's xy_* settings.
* After reading a grid it will have wesn/inc in meters.
* Before writing a grid, it may have units changed back to original units
* or scaled to another set of units */
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (h);
if (direction == GMT_IN) {
if (HH->xy_adjust[direction] == 0) return; /* Nothing to do */
if (HH->xy_adjust[GMT_IN] & 2) return; /* Already scaled them */
for (k = 0; k < 4; k++) h->wesn[k] *= HH->xy_unit_to_meter[GMT_IN];
for (k = 0; k < 2; k++) h->inc[k] *= HH->xy_unit_to_meter[GMT_IN];
HH->xy_adjust[GMT_IN] = 2; /* Now the grid is ready for use and in meters */
if (HH->xy_mode[direction])
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Input grid file x/y unit was converted from meters to %s after reading.\n", GMT->current.proj.unit_name[HH->xy_unit[GMT_IN]]);
else
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Input grid file x/y unit was converted from %s to meters after reading.\n", GMT->current.proj.unit_name[HH->xy_unit[GMT_IN]]);
}
else if (direction == GMT_OUT) { /* grid x/y are assumed to be in meters */
if (HH->xy_adjust[GMT_OUT] & 1) { /* Was given a new unit for output */
for (k = 0; k < 4; k++) h->wesn[k] /= HH->xy_unit_to_meter[GMT_OUT];
for (k = 0; k < 2; k++) h->inc[k] /= HH->xy_unit_to_meter[GMT_OUT];
HH->xy_adjust[GMT_OUT] = 2; /* Now we are ready for writing */
if (HH->xy_mode[GMT_OUT])
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Output grid file x/y unit was converted from %s to meters before writing.\n", GMT->current.proj.unit_name[HH->xy_unit[GMT_OUT]]);
else
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Output grid file x/y unit was converted from meters to %s before writing.\n", GMT->current.proj.unit_name[HH->xy_unit[GMT_OUT]]);
}
else if (HH->xy_adjust[GMT_IN] & 2) { /* Just undo old scaling */
for (k = 0; k < 4; k++) h->wesn[k] /= HH->xy_unit_to_meter[GMT_IN];
for (k = 0; k < 2; k++) h->inc[k] /= HH->xy_unit_to_meter[GMT_IN];
HH->xy_adjust[GMT_IN] -= 2; /* Now it is back to where we started */
if (HH->xy_mode[GMT_OUT])
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Output grid file x/y unit was reverted back to %s from meters before writing.\n", GMT->current.proj.unit_name[HH->xy_unit[GMT_IN]]);
else
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Output grid file x/y unit was reverted back from meters to %s before writing.\n", GMT->current.proj.unit_name[HH->xy_unit[GMT_IN]]);
}
}
}
/* Routines to see if a particular grd file format is specified as part of filename. */
GMT_LOCAL void gmtgrdio_expand_filename (struct GMT_CTRL *GMT, const char *file, char *fname) {
bool found;
unsigned int i;
size_t f_length, length;
if (GMT->current.setting.io_gridfile_shorthand) { /* Look for matches */
f_length = strlen (file);
for (i = 0, found = false; !found && i < GMT->session.n_shorthands; ++i) {
length = strlen (GMT->session.shorthand[i].suffix);
found = (length > f_length) ? false : !strncmp (&file[f_length - length], GMT->session.shorthand[i].suffix, length);
}
if (found) { /* file ended in a recognized shorthand extension */
--i;
sprintf (fname, "%s=%s", file, GMT->session.shorthand[i].format);
}
else
strcpy (fname, file);
}
else /* Simply copy the full name */
strcpy (fname, file);
}
enum Grid_packing_mode {
k_grd_pack = 0, /* scale and offset before writing to disk */
k_grd_unpack /* remove scale and offset after reading packed data */
};
GMT_LOCAL void gmtgrdio_pack_grid (struct GMT_CTRL *Ctrl, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, unsigned pack_mode) {
size_t n_representations = 0; /* number of distinct values >= 0 that a signed integral type can represent */
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
if (pack_mode == k_grd_pack && (HH->z_scale_autoadjust || HH->z_offset_autoadjust)) {
switch (Ctrl->session.grdformat[header->type][1]) {
case 'b':
n_representations = 128; /* exp2 (8 * sizeof (int8_t)) / 2 */
break;
case 's':
n_representations = 32768; /* exp2 (8 * sizeof (int16_t)) / 2 */
break;
case 'i':
/* A single precision float's significand has a precision of 24 bits.
* In order to avoid round-off errors, we must not use all 2^32
* n_representations of an int32_t. */
n_representations = 0x1000000; /* exp2 (24) */
break;
/* default: do not auto-scale floating point */
}
}
if (n_representations != 0) {
/* Calculate auto-scale and offset */
gmt_grd_zminmax (Ctrl, header, grid); /* Calculate z_min/z_max */
if (HH->z_offset_autoadjust) {
/* shift to center values around 0 but shift only by integral value */
double z_range = header->z_max - header->z_min;
if (isfinite (z_range))
header->z_add_offset = rint(z_range / 2.0 + header->z_min);
}
if (HH->z_scale_autoadjust) {
/* scale z-range to use all n_representations */
double z_max = header->z_max - header->z_add_offset;
double z_min = fabs(header->z_min - header->z_add_offset);
double z_0_n_range = MAX (z_max, z_min); /* use [0,n] range because of signed int */
--n_representations; /* subtract 1 for NaN value */
if (isnormal (z_0_n_range))
header->z_scale_factor = z_0_n_range / n_representations;
}
}
/* Do actual packing/unpacking: */
switch (pack_mode) {
case k_grd_unpack:
gmt_scale_and_offset_f (Ctrl, grid, header->size, header->z_scale_factor, header->z_add_offset);
/* Adjust z-range in header: */
header->z_min = header->z_min * header->z_scale_factor + header->z_add_offset;
header->z_max = header->z_max * header->z_scale_factor + header->z_add_offset;
if (header->z_scale_factor < 0.0) gmt_M_double_swap (header->z_min, header->z_max);
break;
case k_grd_pack:
gmt_scale_and_offset_f (Ctrl, grid, header->size, 1.0/header->z_scale_factor, -header->z_add_offset/header->z_scale_factor);
break;
default:
assert (false); /* gmtgrdio_pack_grid() called with illegal pack_mode */
}
}
GMT_LOCAL int gmtgrdio_parse_grd_format_scale_old (struct GMT_CTRL *Ctrl, struct GMT_GRID_HEADER *header, char *format) {
/* parses format string after =-suffix: ff/scale/offset/invalid
* ff: can be one of [abcegnrs][bsifd]
* scale: can be any non-zero normalized number or 'a' for scale and
* offset auto-adjust, defaults to 1.0 if omitted
* offset: can be any finite number or 'a' for offset auto-adjust, defaults to 0 if omitted
* invalid: can be any finite number, defaults to NaN if omitted
* scale and offset may be left empty (e.g., ns//a will auto-adjust the offset only)
*/
char type_code[3];
char *p;
int err; /* gmt_M_err_trap */
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
/* decode grid type */
strncpy (type_code, format, 2);
type_code[2] = '\0';
if (type_code[0] == '/') /* user passed a scale with no id code =/scale[...]. Assume "nf" */
header->type = GMT_GRID_IS_NF;
else {
err = gmt_grd_format_decoder (Ctrl, type_code, &header->type); /* update header type id */
if (err != GMT_NOERROR)
return err;
}
/* parse scale/offset/invalid if any */
p = strchr (format, '/');
if (p != NULL && *p) {
++p;
/* parse scale */
if (*p == 'a')
HH->z_scale_autoadjust = HH->z_offset_autoadjust = true;
else
sscanf (p, "%lf", &header->z_scale_factor);
}
else
return GMT_NOERROR;
p = strchr (p, '/');
if (p != NULL && *p) {
++p;
/* parse offset */
if (*p != 'a') {
HH->z_offset_autoadjust = false;
sscanf (p, "%lf", &header->z_add_offset);
}
else
HH->z_offset_autoadjust = true;
}
else
return GMT_NOERROR;
p = strchr (p, '/');
if (p != NULL && *p) {
++p;
/* parse invalid value */
#ifdef DOUBLE_PRECISION_GRID
sscanf (p, "%lf", &header->nan_value);
#else
sscanf (p, "%f", &header->nan_value);
#endif
/* header->nan_value should be of same type as (float)*grid to avoid
* round-off errors. For example, =gd///-3.4028234e+38:gtiff, would fail
* otherwise because the GTiff'd NoData values are of type double but the
* grid is truncated to float.
* Don't allow infitiy: */
if (!isfinite (header->nan_value))
header->nan_value = (gmt_grdfloat)NAN;
}
return GMT_NOERROR;
}
GMT_LOCAL int gmtgrdio_parse_grd_format_scale_new (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, char *format) {
/* parses format string after = suffix: ff[+d<divisor>][+s<scale>][+o<offset>/][+n<invalid>]
* ff: can be one of [abcegnrs][bsifd]
* scale: can be any non-zero normalized number or 'a' for scale and
* offset auto-adjust, defaults to 1.0 if omitted
* offset: can be any finite number or 'a' for offset auto-adjust, defaults to 0 if omitted
* invalid: can be any finite number, defaults to NaN if omitted
* E.g., ns+oa will auto-adjust the offset only.
*/
char type_code[3];
char p[GMT_BUFSIZ] = {""};
unsigned int pos = 0, uerr = 0;
int err; /* gmt_M_err_trap */
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
/* decode grid type */
strncpy (type_code, format, 2);
type_code[2] = '\0';
if (type_code[0] == '+') /* User passed a scale, offset, or nan-modifier with no id code, e.g., =+s<scale>[...]. Assume "nf" format */
header->type = GMT_GRID_IS_NF;
else { /* Match given code with known codes */
err = gmt_grd_format_decoder (GMT, type_code, &header->type); /* update header type id */
if (err != GMT_NOERROR)
return err;
}
while (gmt_getmodopt (GMT, 0, format, "bdnos", &pos, p, &uerr) && uerr == 0) { /* Looking for +b, +d, +n, +o, +s */
switch (p[0]) {
case 'b': /* bands - not parsed here */
break;
case 'd': /* parse divisor */
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtgrdio_parse_grd_format_scale_new: Setting divisor as %s\n", &p[1]);
header->z_scale_factor = 1.0 / atof (&p[1]); /* Just convert to a scale */
HH->z_scale_given = true;
break;
case 'n': /* Nan value */
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtgrdio_parse_grd_format_scale_new: Using %s to represent missing value (NaN)\n", &p[1]);
header->nan_value = (gmt_grdfloat)atof (&p[1]);
/* header->nan_value should be of same type as (gmt_grdfloat)*grid to avoid
* round-off errors. For example, =gd+n-3.4028234e+38:gtiff, would fail
* otherwise because the GTiff'd NoData values are of type double but the
* grid is truncated to gmt_grdfloat. Don't allow infinity: */
if (!isfinite (header->nan_value))
header->nan_value = (gmt_grdfloat)NAN;
break;
case 'o': /* parse offset */
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtgrdio_parse_grd_format_scale_new: Setting offset as %s\n", &p[1]);
if (p[1] != 'a') {
HH->z_offset_autoadjust = false;
HH->z_offset_given = true;
header->z_add_offset = atof (&p[1]);
}
else
HH->z_offset_autoadjust = true;
break;
case 's': /* parse scale */
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtgrdio_parse_grd_format_scale_new: Setting scale as %s\n", &p[1]);
if (p[1] == 'a') /* This sets both scale and offset to auto */
HH->z_scale_autoadjust = HH->z_offset_autoadjust = true;
else {
header->z_scale_factor = atof (&p[1]);
HH->z_scale_given = true;
}
break;
default: /* These are caught in gmt_getmodopt so break is just for Coverity */
break;
}
}
return (uerr) ? GMT_NOT_A_VALID_MODIFIER : GMT_NOERROR;
}
GMT_LOCAL int gmtgrdio_parse_grd_format_scale (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, char *format) {
int code;
if (gmt_found_modifier (GMT, format, "bdnos"))
code = gmtgrdio_parse_grd_format_scale_new (GMT, header, format);
else
code = gmtgrdio_parse_grd_format_scale_old (GMT, header, format);
return (code);
}
GMT_LOCAL bool gmtgrdio_eq (double this, double that, double inc) {
/* this and that are the same value if less than 10e-4 * inc apart */
return ((fabs (this - that) / inc) < GMT_CONV4_LIMIT);
}
GMT_LOCAL int gmtgrdio_padspace (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, double *wesn, bool image, unsigned int *pad, struct GRD_PAD *P) {
/* When padding is requested it is usually used to set boundary conditions based on
* two extra rows/columns around the domain of interest. BCs like natural or periodic
* can then be used to fill in the pad. However, if the domain is taken from a grid
* whose full domain exceeds the region of interest we are better off using the extra
* data to fill those pad rows/columns. Thus, this function tries to determine if the
* input grid has the extra data we need to fill the BC pad with observations
* P returns the correct region and pad to use regardless if function returns true or false */
bool wrap;
unsigned int side, n_sides = 0;
double wesn2[4];
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
gmt_M_unused(GMT);
/* First copy over original settings to the Pad structure */
gmt_M_memset (P, 1, struct GRD_PAD); /* Initialize to zero */
gmt_M_memcpy (P->pad, pad, 4, int); /* Duplicate the pad */
gmt_M_memcpy (P->wesn, header->wesn, 4, double); /* Copy the header boundaries */
if (!wesn) return (false); /* No subset requested */
if (wesn[XLO] == wesn[XHI] && wesn[YLO] == wesn[YHI]) return (false); /* Subset not set */
if (gmtgrdio_eq (wesn[XLO], header->wesn[XLO], header->inc[GMT_X]) && gmtgrdio_eq (wesn[XHI], header->wesn[XHI], header->inc[GMT_X])
&& gmtgrdio_eq (wesn[YLO], header->wesn[YLO], header->inc[GMT_Y]) && gmtgrdio_eq (wesn[YHI], header->wesn[YHI], header->inc[GMT_Y]))
return (false); /* Subset equals whole area */
gmt_M_memcpy (P->wesn, wesn, 4, double); /* Copy the subset boundaries */
if (pad[XLO] == 0 && pad[XHI] == 0 && pad[YLO] == 0 && pad[YHI] == 0) return (false); /* No padding requested */
if (!GMT->current.io.grid_padding) return (false); /* Not requested */
/* Determine if data exist for a pad on all four sides. If not we give up */
wrap = !image && gmt_grd_is_global (GMT, header); /* If global grid wrap then we cannot be outside */
if ((wesn2[XLO] = wesn[XLO] - pad[XLO] * header->inc[GMT_X]) < header->wesn[XLO] && !wrap) /* Cannot extend west/xmin */
{ n_sides++; wesn2[XLO] = wesn[XLO]; }
else /* OK to load left pad with data */
P->pad[XLO] = 0;
if ((wesn2[XHI] = wesn[XHI] + pad[XHI] * header->inc[GMT_X]) > header->wesn[XHI] && !wrap) /* Cannot extend east/xmax */
{ n_sides++; wesn2[XHI] = wesn[XHI]; }
else /* OK to load right pad with data */
P->pad[XHI] = 0;
if ((wesn2[YLO] = wesn[YLO] - pad[YLO] * header->inc[GMT_Y]) < header->wesn[YLO]) /* Cannot extend south/ymin */
{ n_sides++; wesn2[YLO] = wesn[YLO]; }
else /* OK to load bottom pad with data */
P->pad[YLO] = 0;
if ((wesn2[YHI] = wesn[YHI] + pad[YHI] * header->inc[GMT_Y]) > header->wesn[YHI]) /* Cannot extend north/ymax */
{ n_sides++; wesn2[YHI] = wesn[YHI]; }
else /* OK to load top pad with data */
P->pad[YHI] = 0;
if (n_sides == 4) return (false); /* No can do */
/* Here we know that there is enough input data to fill some or all of the BC pad with actual data values */
/* We have temporarily set padding to zero (since the pad is now part of the region) for those sides we can extend */
/* Temporarily enlarge the region so it now includes the padding we need */
gmt_M_memcpy (P->wesn, wesn2, 4, double);
/* Set BC */
for (side = 0; side < 4; side++) {
if (P->pad[side] == 0)
HH->BC[side] = GMT_BC_IS_DATA;
}
return (true); /* Return true so the calling function can take appropriate action */
}
void gmt_grd_set_datapadding (struct GMT_CTRL *GMT, bool set) {
/* Changes the value of GMT->current.io.grid_padding.
* If set = true we turn padding on, else off. */
GMT->current.io.grid_padding = set;
}
GMT_LOCAL void gmtgrdio_handle_pole_averaging (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, gmt_grdfloat f_value, int pole) {
uint64_t node;
openmp_int col = 0;
char *name[3] = {"south", "", "north"};
gmt_M_unused(GMT);
if (pole == -1)
node = gmt_M_ijp (header, header->n_rows-1, 0); /* First node at S pole */
else
node = gmt_M_ijp (header, 0, 0); /* First node at N pole */
if (gmt_M_type (GMT, GMT_OUT, GMT_Z) == GMT_IS_AZIMUTH || gmt_M_type (GMT, GMT_OUT, GMT_Z) == GMT_IS_ANGLE) { /* Must average azimuths */
uint64_t orig = node;
double s, c, sum_s = 0.0, sum_c = 0.0;
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Average %d angles at the %s pole\n", header->n_columns, name[pole+1]);
for (col = 0; col < (openmp_int)header->n_columns; col++, node++) {
sincosd ((double)grid[node], &s, &c);
sum_s += s; sum_c += c;
}
f_value = (gmt_grdfloat)atan2d (sum_s, sum_c);
node = orig;
}
for (col = 0; col < (openmp_int)header->n_columns; col++, node++) grid[node] = f_value;
}
GMT_LOCAL void gmtgrdio_grd_check_consistency (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid) {
/* Enforce before writing a grid that periodic grids with repeating columns
* agree on the node values in those columns; if different replace with average.
* This only affects geographic grids of 360-degree extent with gridline registration.
* Also, if geographic grid with gridline registration, if the N or S pole row is present
* we ensure that they all have identical values, otherwise replace by mean value */
openmp_int row = 0, col = 0;
unsigned int we_conflicts = 0, p_conflicts = 0;
uint64_t left = 0, right = 0, node = 0;
if (header->registration == GMT_GRID_PIXEL_REG) return; /* Not gridline registered */
if (gmt_M_is_cartesian (GMT, GMT_OUT)) return; /* Not geographic */
if (header->wesn[YLO] == -90.0) { /* Check consistency of S pole duplicates */
double sum;
node = gmt_M_ijp (header, header->n_rows-1, 0); /* First node at S pole */
sum = grid[node++];
p_conflicts = 0;
for (col = 1; col < (openmp_int)header->n_columns; col++, node++) {
if (grid[node] != grid[node-1]) p_conflicts++;
sum += grid[node];
}
if (p_conflicts) {
gmt_grdfloat f_value = (gmt_grdfloat)(sum / header->n_columns);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected %u inconsistent values at south pole. Values fixed by setting all to average row value.\n", p_conflicts);
gmtgrdio_handle_pole_averaging (GMT, header, grid, f_value, -1);
}
}
if (header->wesn[YHI] == +90.0) { /* Check consistency of N pole duplicates */
double sum;
node = gmt_M_ijp (header, 0, 0); /* First node at N pole */
sum = grid[node++];
p_conflicts = 0;
for (col = 1; col < (openmp_int)header->n_columns; col++, node++) {
if (grid[node] != grid[node-1]) p_conflicts++;
sum += grid[node];
}
if (p_conflicts) {
gmt_grdfloat f_value = (gmt_grdfloat)(sum / header->n_columns);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected %u inconsistent values at north pole. Values fixed by setting all to average row value.\n", p_conflicts);
gmtgrdio_handle_pole_averaging (GMT, header, grid, f_value, +1);
}
}
if (!gmt_M_360_range (header->wesn[XLO], header->wesn[XHI])) return; /* Not 360-degree range */
for (row = 0; row < (openmp_int)header->n_rows; row++) {
left = gmt_M_ijp (header, row, 0); /* Left node */
right = left + header->n_columns - 1; /* Right node */
if (grid[right] != grid[left]) {
grid[right] = grid[left];
we_conflicts++;
}
}
if (we_conflicts)
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected %u inconsistent values along periodic east boundary of grid. Values fixed by duplicating west boundary.\n", we_conflicts);
}
GMT_LOCAL void gmtgrdio_grd_wipe_pad (struct GMT_CTRL *GMT, struct GMT_GRID *G) {
/* Reset padded areas to 0. */
unsigned int row;
size_t ij0;
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (G->header);
if (HH->arrangement == GMT_GRID_IS_INTERLEAVED) {
GMT_Report (GMT->parent, GMT_MSG_ERROR, "Calling gmtgrdio_grd_wipe_pad on interleaved complex grid! Programming error?\n");
return;
}
if (G->header->pad[YHI]) gmt_M_memset (G->data, G->header->mx * G->header->pad[YHI], gmt_grdfloat); /* Wipe top pad */
if (G->header->pad[YLO]) { /* Wipe bottom pad */
ij0 = gmt_M_ij (G->header, G->header->my - G->header->pad[YLO], 0); /* Index of start of bottom pad */
gmt_M_memset (&(G->data[ij0]), G->header->mx * G->header->pad[YLO], gmt_grdfloat);
}
if (G->header->pad[XLO] == 0 && G->header->pad[XHI] == 0) return; /* Nothing to do */
for (row = G->header->pad[YHI]; row < G->header->my - G->header->pad[YLO]; row++) { /* Wipe left and right pad which is trickier */
ij0 = gmt_M_ij (G->header, row, 0); /* Index of this row's left column (1st entry in west pad) */
if (G->header->pad[XLO]) gmt_M_memset (&(G->data[ij0]), G->header->pad[XLO], gmt_grdfloat);
ij0 += (G->header->mx - G->header->pad[XHI]); /* Start of this rows east pad's 1st column */
if (G->header->pad[XHI]) gmt_M_memset (&(G->data[ij0]), G->header->pad[XHI], gmt_grdfloat);
}
}
GMT_LOCAL void gmtgrdio_pad_grd_off_sub (struct GMT_GRID *G, gmt_grdfloat *data) {
/* Remove the current grid pad and shuffle all rows to the left */
uint64_t ijp, ij0;
openmp_int row;
for (row = 0; row < (openmp_int)G->header->n_rows; row++) {
ijp = gmt_M_ijp (G->header, row, 0); /* Index of start of this row's first column in padded grid */
ij0 = gmt_M_ij0 (G->header, row, 0); /* Index of start of this row's first column in unpadded grid */
gmt_M_memcpy (&(data[ij0]), &(data[ijp]), G->header->n_columns, gmt_grdfloat); /* Only copy the n_columns data values */
}
}
GMT_LOCAL void gmtgrdio_pad_grd_on_sub (struct GMT_CTRL *GMT, struct GMT_GRID *G, struct GMT_GRID_HEADER *h_old, gmt_grdfloat *data) {
/* Use G for dimensions but operate on data array which points to either the real or imaginary section */
uint64_t ij_new, ij_old;
openmp_int row, col, start_last_new_row, end_last_old_row;
/* See if the index of start of last new row exceeds index of last node in old grid */
start_last_new_row = gmt_M_get_nm (GMT, G->header->pad[YHI] + G->header->n_rows - 1, G->header->pad[XLO] + G->header->n_columns + G->header->pad[XHI]) + G->header->pad[XLO];
end_last_old_row = gmt_M_get_nm (GMT, h_old->pad[YHI] + h_old->n_rows - 1, h_old->pad[XLO] + h_old->n_columns + h_old->pad[XHI]) + h_old->pad[XLO] + h_old->n_columns - 1;
if (start_last_new_row > end_last_old_row && (start_last_new_row - end_last_old_row) > (openmp_int)G->header->n_columns) { /* May copy whole rows from bottom to top */
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtgrdio_pad_grd_on_sub can copy row-by-row\n");
for (row = (openmp_int)G->header->n_rows; row > 0; row--) {
ij_new = gmt_M_ijp (G->header, row-1, 0); /* Index of this row's first column in new padded grid */
ij_old = gmt_M_ijp (h_old, row-1, 0); /* Index of this row's first column in old padded grid */
gmt_M_memcpy (&(data[ij_new]), &(data[ij_old]), G->header->n_columns, gmt_grdfloat);
}
}
else { /* Must do it from bottom to top on a per node basis */
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtgrdio_pad_grd_on_sub must copy node-by-node\n");
ij_new = gmt_M_ijp (G->header, G->header->n_rows-1, G->header->n_columns-1); /* Index of this row's last column in new padded grid */
ij_old = gmt_M_ijp (h_old, h_old->n_rows-1, h_old->n_columns-1); /* Index of this row's last column in old padded grid */
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtgrdio_pad_grd_on_sub: last ij_new = %d, last ij_old = %d\n", (int)ij_new, (int)ij_old);
if (ij_new > ij_old) { /* Can go back to front */
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtgrdio_pad_grd_on_sub: Must loop from end to front\n");
for (row = (openmp_int)G->header->n_rows; row > 0; row--) {
ij_new = gmt_M_ijp (G->header, row-1, G->header->n_columns-1); /* Index of this row's last column in new padded grid */
ij_old = gmt_M_ijp (h_old, row-1, h_old->n_columns-1); /* Index of this row's last column in old padded grid */
for (col = 0; col < (openmp_int)G->header->n_columns; col++)
data[ij_new--] = data[ij_old--];
}
}
else {
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtgrdio_pad_grd_on_sub: Must loop from front to end\n");
for (row = 0; row < (openmp_int)G->header->n_rows; row++) {
ij_new = gmt_M_ijp (G->header, row, 0); /* Index of this row's last column in new padded grid */
ij_old = gmt_M_ijp (h_old, row, 0); /* Index of this row's last column in old padded grid */
for (col = 0; col < (openmp_int)G->header->n_columns; col++)
data[ij_new++] = data[ij_old++];
}
}
}
gmtgrdio_grd_wipe_pad (GMT, G); /* Set pad areas to 0 */
}
GMT_LOCAL void gmtgrdio_pad_grd_zero_sub (struct GMT_GRID *G, gmt_grdfloat *data) {
openmp_int row, col;
unsigned int nx1;
uint64_t ij_f, ij_l;
if (G->header->pad[YHI]) gmt_M_memset (data, G->header->pad[YHI] * G->header->mx, gmt_grdfloat); /* Zero the top pad */
nx1 = G->header->n_columns - 1; /* Last column */
gmt_M_row_loop (NULL, G, row) {
ij_f = gmt_M_ijp (G->header, row, 0); /* Index of first column this row */
ij_l = gmt_M_ijp (G->header, row, nx1); /* Index of last column this row */
for (col = 1; col <= (openmp_int)G->header->pad[XLO]; col++) data[ij_f-col] = 0.0f; /* Zero the left pad at this row */
for (col = 1; col <= (openmp_int)G->header->pad[XHI]; col++) data[ij_l+col] = 0.0f; /* Zero the left pad at this row */
}
if (G->header->pad[YLO]) {
int pad = G->header->pad[XLO];
ij_f = gmt_M_ijp (G->header, G->header->n_rows, -pad); /* Index of first column of bottom pad */
gmt_M_memset (&(data[ij_f]), G->header->pad[YLO] * G->header->mx, gmt_grdfloat); /* Zero the bottom pad */
}
}
GMT_LOCAL struct GMT_GRID_HEADER *gmtgrdio_duplicate_gridheader (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h) {
/* Duplicates a grid header. */
struct GMT_GRID_HEADER *hnew = gmt_get_header (GMT);
gmt_copy_gridheader (GMT, hnew, h);
return (hnew);
}
/*----------------------------------------------------------|
* Public functions that are part of the GMT Devel library |
*----------------------------------------------------------|
*/
#ifdef DEBUG
/* Uncomment this to have gmt_grd_dump be called and do something */
/* #define GMT_DUMPING */
#ifdef GMT_DUMPING
void gmt_grd_dump (struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, bool is_complex, char *txt) {
unsigned int row, col;
uint64_t k = 0U;
fprintf (stderr, "Dump [%s]:\n---------------------------------------------\n", txt);
for (row = 0; row < header->my; row++) {
if (is_complex)
for (col = 0; col < header->mx; col++, k+= 2) fprintf (stderr, "(%g,%g)\t", grid[k], grid[k+1]);
else
for (col = 0; col < header->mx; col++, k++) fprintf (stderr, "%g\t", grid[k]);
fprintf (stderr, "\n");
}
fprintf (stderr, "---------------------------------------------\n");
}
#else /* Just a dummy */
void gmt_grd_dump (struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, bool is_complex, char *txt) {
gmt_M_unused(header); gmt_M_unused(grid); gmt_M_unused(is_complex); gmt_M_unused(txt);
/* Nothing */
}
#endif
#endif
void gmt_copy_gridheader (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *to, struct GMT_GRID_HEADER *from) {
/* Destination must exist */
struct GMT_GRID_HEADER_HIDDEN *Hfrom = gmt_get_H_hidden (from), *Hto = gmt_get_H_hidden (to);
gmt_M_unused(GMT);
if (GMT->parent->internal && to->ProjRefWKT) gmt_M_str_free (to->ProjRefWKT); /* Since we will duplicate via from */
if (GMT->parent->internal && to->ProjRefPROJ4) gmt_M_str_free (to->ProjRefPROJ4); /* Since we will duplicate via from */
if (Hto->pocket) gmt_M_str_free (Hto->pocket); /* Since we will duplicate via from */
if (Hto->title) gmt_M_str_free (Hto->title); /* Since we will duplicate via from */
if (Hto->command) gmt_M_str_free (Hto->command); /* Since we will duplicate via from */
if (Hto->remark) gmt_M_str_free (Hto->remark); /* Since we will duplicate via from */
if (Hto->cpt) gmt_M_str_free (Hto->cpt); /* Since we will duplicate via from */
gmt_M_memcpy (to, from, 1, struct GMT_GRID_HEADER); /* Copies full contents but also duplicates the hidden address */
to->hidden = Hto; /* Restore the original hidden address in to */
gmt_M_memcpy (to->hidden, from->hidden, 1, struct GMT_GRID_HEADER_HIDDEN); /* Copies full contents of hidden area */
/* Must deal with six pointers individually */
if (from->ProjRefWKT) to->ProjRefWKT = strdup (from->ProjRefWKT);
if (from->ProjRefPROJ4) to->ProjRefPROJ4 = strdup (from->ProjRefPROJ4);
if (Hfrom->pocket) Hto->pocket = strdup (Hfrom->pocket);
if (Hfrom->title) Hto->title = strdup (Hfrom->title);
if (Hfrom->command) Hto->command = strdup (Hfrom->command);
if (Hfrom->remark) Hto->remark = strdup (Hfrom->remark);
if (Hfrom->cpt) Hto->cpt = strdup (Hfrom->cpt);
}
/*! gmt_grd_is_global returns true for a geographic grid with exactly 360-degree range (with or without repeating column) */
bool gmt_grd_is_global (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h) {
bool global;
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (h);
gmt_M_unused(GMT);
global = (HH->grdtype == GMT_GRID_GEOGRAPHIC_EXACT360_NOREPEAT || HH->grdtype == GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT);
return (global);
}
/*! gmt_grd_is_polar returns true for a geographic grid with exactly 180-degree range in latitude (with or without repeating column) */
bool gmt_grd_is_polar (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h) {
if (!gmt_M_y_is_lat (GMT, GMT_IN))
return false;
if (gmt_M_180_range (h->wesn[YHI], h->wesn[YLO]))
return true;
if (fabs (h->n_rows * h->inc[GMT_Y] - 180.0) < GMT_CONV4_LIMIT)
return true;
return false;
}
void gmt_set_R_from_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
/* When no -R was given we will inherit the region from the grid. However,
* many grids are hobbled by not clearly specifying they are truly global grids.
* What frequently happens is that gridnode-registered grids omit the repeating
* column in the east, leading to regions such as -R0/359/-90/90 for a 1-degree grid.
* Since these are clearly global we do now want to pass 0/359 to the projection
* machinery but 0/360. Hence we test if the grid is truly global and make this decision. */
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
gmt_M_memcpy (GMT->common.R.wesn, header->wesn, 4, double); /* Initially we set -R as is from grid header */
if (HH->grdtype != GMT_GRID_GEOGRAPHIC_EXACT360_NOREPEAT) return; /* Nothing to do */
if (!gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]) && fabs (header->n_columns * header->inc[GMT_X] - 360.0) < GMT_CONV4_LIMIT) {
/* The w/e need to state the complete 360 range: Let east = 360 + west */
GMT->common.R.wesn[XHI] = GMT->common.R.wesn[XLO] + 360.0;
}
if (!gmt_M_180_range (GMT->common.R.wesn[YLO], GMT->common.R.wesn[YHI]) && fabs (header->n_rows * header->inc[GMT_Y] - 180.0) < GMT_CONV4_LIMIT) {
/* The s/n need to state the complete 180 range */
GMT->common.R.wesn[YLO] = -90.0;
GMT->common.R.wesn[YHI] = +90.0;
}
}
void gmtlib_grd_get_units (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
/* Set input data types for columns 0, 1 and 2 based on unit strings for
grid coordinates x, y and z.
When "Time": transform the data scale and offset to match the current time system.
*/
unsigned int i;
char string[3][GMT_LEN256], *units = NULL;
double scale = 1.0, offset = 0.0;
struct GMT_TIME_SYSTEM time_system;
/* Copy unit strings */
strncpy (string[0], header->x_units, GMT_GRID_UNIT_LEN80);
strncpy (string[1], header->y_units, GMT_GRID_UNIT_LEN80);
strncpy (string[2], header->z_units, GMT_GRID_UNIT_LEN80);
/* Parse the unit strings one by one */
for (i = 0; i < 3; i++) {
/* Skip parsing when input data type is already set */
if (GMT->current.io.col_type[GMT_IN][i] & GMT_IS_GEO) continue;
if (GMT->current.io.col_type[GMT_IN][i] & GMT_IS_RATIME) {
GMT->current.proj.xyz_projection[i] = GMT_TIME;
continue;
}
/* Change name of variable and unit to lower case for comparison */
gmt_str_tolower (string[i]);
/* PW: This test was <= 360 but if you have grids from 340-375 and it says longitude then we should not make it Cartesian */
if ((!strncmp (string[i], "longitude", 9U) || strstr (string[i], "degrees_e")) && (header->wesn[XLO] > -360.0 && header->wesn[XHI] < 720.0)) {
/* Input data type is longitude ... if */
if ((header->wesn[XHI] - header->wesn[XLO]) > 360)
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Grid x units says it is longitude but it spans > 360 so is set to Cartesian\n");
else
gmt_set_column_type (GMT, GMT_IN, i, GMT_IS_LON);
}
else if ((!strncmp (string[i], "latitude", 8U) || strstr (string[i], "degrees_n")) && (header->wesn[YLO] >= -90.0 && header->wesn[YHI] <= 90.0)) {
/* Input data type is latitude */
gmt_set_column_type (GMT, GMT_IN, i, GMT_IS_LAT);
}
else if (!strcmp (string[i], "time") || !strncmp (string[i], "time [", 6U)) {
/* Input data type is time */
gmt_set_column_type (GMT, GMT_IN, i, GMT_IS_RELTIME);
GMT->current.proj.xyz_projection[i] = GMT_TIME;
/* Determine coordinates epoch and units (default is internal system) */
gmt_M_memcpy (&time_system, &GMT->current.setting.time_system, 1, struct GMT_TIME_SYSTEM);
units = strchr (string[i], '[');
if (!units || gmt_get_time_system (GMT, ++units, &time_system) || gmt_init_time_system_structure (GMT, &time_system))
GMT_Report (GMT->parent, GMT_MSG_WARNING, "Time units [%s] in grid not recognized, defaulting to %s settings.\n", units, GMT_SETTINGS_FILE);
/* Determine scale between grid and internal time system, as well as the offset (in internal units) */
scale = time_system.scale * GMT->current.setting.time_system.i_scale;
offset = (time_system.rata_die - GMT->current.setting.time_system.rata_die) + (time_system.epoch_t0 - GMT->current.setting.time_system.epoch_t0);
offset *= GMT_DAY2SEC_F * GMT->current.setting.time_system.i_scale;
/* Scale data scale and extremes based on scale and offset */
if (i == 0) {
header->wesn[XLO] = header->wesn[XLO] * scale + offset;
header->wesn[XHI] = header->wesn[XHI] * scale + offset;
header->inc[GMT_X] *= scale;
}
else if (i == 1) {
header->wesn[YLO] = header->wesn[YLO] * scale + offset;
header->wesn[YHI] = header->wesn[YHI] * scale + offset;
header->inc[GMT_Y] *= scale;
}
else {
header->z_add_offset = header->z_add_offset * scale + offset;
header->z_scale_factor *= scale;
}
}
}
}
double * gmt_grd_coord (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, int dir) {
/* Allocate, compute, and return the x- or y-coordinates for a grid */
unsigned int k;
double *coord = NULL;
assert (dir == GMT_X || dir == GMT_Y);
if (dir == GMT_X) {
if ((coord = gmt_M_memory (GMT, NULL, header->n_columns, double)) == NULL) return NULL;
for (k = 0; k < header->n_columns; k++) coord[k] = gmt_M_grd_col_to_x (GMT, k, header);
}
else if (dir == GMT_Y) {
if ((coord = gmt_M_memory (GMT, NULL, header->n_rows, double)) == NULL) return NULL;
for (k = 0; k < header->n_rows; k++) coord[k] = gmt_M_grd_row_to_y (GMT, k, header);
}
return (coord);
}
void gmtlib_grd_real_interleave (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *data) {
/* Sub-function since this step is also needed in the specific case when a module
* calls GMT_Read_Data with GMT_GRID_IS_COMPLEX_REAL but input comes via an
* external interface (MATLAB, Python) and thus has not been multiplexed yet.
* We assume the data array is already large enough to hold the double-sized grid.
*/
uint64_t row, col, col_1, col_2, left_node_1, left_node_2;
gmt_M_unused(GMT);
/* Here we have a grid with RRRRRR..._________ and want R_R_R_R_... */
for (row = header->my; row > 0; row--) { /* Going from last to first row */
left_node_1 = gmt_M_ij (header, row-1, 0); /* Start of row in RRRRR layout */
left_node_2 = 2 * left_node_1; /* Start of same row in R_R_R_ layout */
for (col = header->mx, col_1 = col - 1, col_2 = 2*col - 1; col > 0; col--, col_1--) { /* Go from right to left */
data[left_node_2+col_2] = 0.0f; col_2--; /* Set the Imag component to zero */
data[left_node_2+col_2] = data[left_node_1+col_1]; col_2--;
}
}
}
void gmt_grd_mux_demux (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *data, unsigned int desired_mode) {
/* Multiplex and demultiplex complex grids.
* Complex grids are read/written by dealing with just one component: real or imag.
* Thus, at read|write the developer must specify which component (GMT_CMPLX_REAL|IMAG).
* For fast disk i/o we read complex data in serial. I.e., if we ask for GMT_CMPLX_REAL
* then the array will contain RRRRR....________, where ______ is unused space for the
* imaginary components. Likewise, if we requested GMT_CMPLX_IMAG then the array will
* be returned as _______...IIIIIII....
* Operations like FFTs typically required the data to be interleaved, i.e., in the
* form RIRIRIRI.... Then, when the FFT work is done and we wish to write out the
* result we will need to demultiplex the array back to its serial RRRRR....IIIII
* format before writing takes place.
* gmt_grd_mux_demux performs either multiplex or demultiplex, depending on desired_mode.
* If grid is not complex then we just return doing nothing.
* Note: At this point the grid is mx * my and we visit all the nodes, including the pads.
* hence we use header->mx/my and gmt_M_ij below.
*/
uint64_t row, col, col_1, col_2, left_node_1, left_node_2, offset, ij, ij2;
gmt_grdfloat *array = NULL;
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
if (! (desired_mode == GMT_GRID_IS_INTERLEAVED || desired_mode == GMT_GRID_IS_SERIAL)) {
GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmt_grd_mux_demux called with inappropriate mode - skipped.\n");
return;
}
if ((header->complex_mode & GMT_GRID_IS_COMPLEX_MASK) == 0) return; /* Nuthin' to do */
if (HH->arrangement == desired_mode) return; /* Already has the right layout */
/* In most cases we will actually have RRRRR...______ or _____...IIIII..
* which means half the array is empty and it is easy to shuffle. However,
* in the case with actual RRRR...IIII there is no simple way to do this inplace;
* see http://stackoverflow.com/questions/1777901/array-interleaving-problem */
if (desired_mode == GMT_GRID_IS_INTERLEAVED) { /* Must multiplex the grid */
if ((header->complex_mode & GMT_GRID_IS_COMPLEX_MASK) == GMT_GRID_IS_COMPLEX_MASK) {
/* Transform from RRRRR...IIIII to RIRIRIRIRI... */
/* Implement properly later; for now waste memory by duplicating [memory is cheap and plentiful] */
/* One advantage is that the padding is all zero by virtue of the allocation */
if ((array = gmt_M_memory_aligned (GMT, NULL, header->size, gmt_grdfloat)) == NULL) return;
offset = header->size / 2; /* Position of 1st row in imag portion of RRRR...IIII... */
for (row = 0; row < header->my; row++) { /* Going from first to last row */
for (col = 0; col < header->mx; col++) {
ij = gmt_M_ij (header, row, col); /* Position of an 'R' in the RRRRR portion */
ij2 = 2 * ij;