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 * grd.h contains the definition for a GMT-SYSTEM Version >= 2 grd file 19 * 20 * grd is stored in rows going from west (xmin) to east (xmax) 21 * first row in file has yvalue = north (ymax). 22 * This is SCANLINE orientation. 23 * 24 * Author: Paul Wessel 25 * Date: 1-JAN-2010 26 * Version: 6 API 27 */ 28 29 /*! 30 * \file gmt_grd.h 31 * \brief Definition for a GMT-SYSTEM Version >= 2 grd file 32 */ 33 34 #ifndef GMT_GRID_H 35 #define GMT_GRID_H 36 37 /* netcdf convention */ 38 #define GMT_NC_CONVENTION "CF-1.7" 39 40 enum GMT_enum_grdtype { /* Used in gmt_grdio.c and gmt_nc.c */ 41 /* Special cases of geographic grids with periodicity */ 42 GMT_GRID_CARTESIAN=0, /* Cartesian data, no periodicity involved */ 43 GMT_GRID_GEOGRAPHIC_LESS360, /* x is longitude, but range is < 360 degrees */ 44 GMT_GRID_GEOGRAPHIC_EXACT360_NOREPEAT, /* x is longitude, range is 360 degrees, no repeat node */ 45 GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT, /* x is longitude, range is 360 degrees, gridline registered and repeat node at 360*/ 46 GMT_GRID_GEOGRAPHIC_MORE360 /* x is longitude, and range exceeds 360 degrees */ 47 }; 48 49 /*! Grid layout for complex grids */ 50 enum gmt_enum_grdlayout { /* Needed by some modules */ 51 GMT_GRID_IS_SERIAL = 0, /* Grid is RRRRRR...[IIIIII...] */ 52 GMT_GRID_IS_INTERLEAVED = 1}; /* Grid is RIRIRIRI... - required layout for FFT */ 53 54 enum gmt_enum_grid_nans { /* Flags to tell if a grid has NaNs */ 55 GMT_GRID_NOT_CHECKED = 0, 56 GMT_GRID_NO_NANS = 1, 57 GMT_GRID_HAS_NANS = 2}; 58 /* 59 * GMT's internal representation of grids is north-up, i.e., the index of the 60 * least dimension (aka y or lat) increases to the south. NetCDF files are 61 * usually written bottom-up, i.e., the index of the least dimension increases 62 * from south to north (k_nc_start_south): 63 * 64 * k_nc_start_north: k_nc_start_south: 65 * 66 * y ^ 0 1 2 -------> x 67 * | 3 4 5 | 0 1 2 68 * | 6 7 8 | 3 4 5 69 * -------> x y V 6 7 8 70 */ 71 72 /*! Order of rows in z variable */ 73 enum Netcdf_row_order { /* Used in gmt_grdio.c and gmt_nc.c */ 74 k_nc_start_north = -1, /* The least dimension (i.e., lat or y) decreases */ 75 k_nc_start_south = 1 /* The least dimension (i.e., lat or y) increases */ 76 }; 77 78 enum Netcdf_chunksize { /* Used in gmt_grdio.c and gmt_nc.c */ 79 k_netcdf_io_classic = 0, /* netCDF classic format */ 80 k_netcdf_io_chunked_auto /* netCDF 4 auto-determined optimal chunk size */ 81 }; 82 83 /* The array wesn in the header has a name that indicates the order (west, east, south, north). 84 * However, to avoid using confusing indices 0-5 we define very brief constants 85 * XLO, XHI, YLO, YHI, ZLO, ZHI that should be used instead: */ 86 enum gmt_enum_wesnids { 87 XLO = 0, /* Index for west or xmin value */ 88 XHI, /* Index for east or xmax value */ 89 YLO, /* Index for south or ymin value */ 90 YHI, /* Index for north or ymax value */ 91 ZLO, /* Index for zmin value */ 92 ZHI /* Index for zmax value */ 93 }; 94 95 /* These macros should be used to convert between (column,row) and (x,y). It will eliminate 96 * one source of typos and errors, and since macros are done at compilation time there is no 97 * overhead. Note: gmt_M_x_to_col does not need n_columns but we included it for symmetry reasons. 98 * gmt_M_y_to_row must first compute j', the number of rows in the increasing y-direction (to 99 * match the sense of truncation used for x) then we revert to row number increasing down 100 * by flipping: j = n_rows - 1 - j'. 101 * Note that input col, row _may_ be negative, hence we do the cast to (int) here. */ 102 103 #define gmt_M_x_to_col(x,x0,dx,off,n_columns) (irint((((x) - (x0)) / (dx)) - (off))) 104 #define gmt_M_y_to_row(y,y0,dy,off,n_rows) ((n_rows) - 1 - irint(((((y) - (y0)) / (dy)) - (off)))) 105 #define gmt_M_col_to_x(C,col,x0,x1,dx,off,n_columns) (((int)(col) == (int)((n_columns)-1)) ? (x1) - (off) * (dx) : (x0) + ((col) + (off)) * (dx)) 106 #define gmt_M_row_to_y(C,row,y0,y1,dy,off,n_rows) (((int)(row) == (int)((n_rows)-1)) ? (y0) + (off) * (dy) : (y1) - ((row) + (off)) * (dy)) 107 108 /*! The follow macros simplify using the 4 above macros when all info is in the struct header h. */ 109 110 #define gmt_M_grd_col_to_x(C,col,h) gmt_M_col_to_x(C,col,h->wesn[XLO],h->wesn[XHI],h->inc[GMT_X],h->xy_off,h->n_columns) 111 #define gmt_M_grd_row_to_y(C,row,h) gmt_M_row_to_y(C,row,h->wesn[YLO],h->wesn[YHI],h->inc[GMT_Y],h->xy_off,h->n_rows) 112 #define gmt_M_grd_x_to_col(C,x,h) gmt_M_x_to_col(x,h->wesn[XLO],h->inc[GMT_X],h->xy_off,h->n_columns) 113 #define gmt_M_grd_y_to_row(C,y,h) gmt_M_y_to_row(y,h->wesn[YLO],h->inc[GMT_Y],h->xy_off,h->n_rows) 114 115 /*! These macros calculate the number of nodes in x or y for the increment dx, dy */ 116 117 #define gmt_M_get_n(C,min,max,inc,off) (urint ((((max) - (min)) / (inc)) + 1 - (off)) ) 118 #define gmt_M_get_inc(C,min,max,n,off) (((n) + (off) - 1 == 0) ? ((max) - (min)) : ((max) - (min)) / ((n) + (off) - 1)) 119 120 /*! The follow macros simplify using the 2 above macros when all info is in the struct header */ 121 122 #define gmt_M_grd_get_nx(C,h) (gmt_M_get_n(C,h->wesn[XLO],h->wesn[XHI],h->inc[GMT_X],h->registration)) 123 #define gmt_M_grd_get_ny(C,h) (gmt_M_get_n(C,h->wesn[YLO],h->wesn[YHI],h->inc[GMT_Y],h->registration)) 124 125 /*! The follow macros gets the full length or rows and columns when padding is considered (i.e., mx and my) */ 126 127 #define gmt_M_grd_get_nxpad(h,pad) ((h->n_columns) + pad[XLO] + pad[XHI]) 128 #define gmt_M_grd_get_nypad(h,pad) ((h->n_rows) + pad[YLO] + pad[YHI]) 129 130 /*! 64-bit-safe macros to return the number of points in the grid given its dimensions */ 131 132 #define gmt_M_get_nm(C,n_columns,n_rows) (((uint64_t)(n_columns)) * ((uint64_t)(n_rows))) 133 #define gmt_M_grd_get_nm(h) (((uint64_t)(h->n_columns)) * ((uint64_t)(h->n_rows))) 134 135 /*! gmt_M_grd_setpad copies the given pad into the header */ 136 137 #define gmt_M_grd_setpad(C,h,newpad) memcpy ((h)->pad, newpad, 4*sizeof(unsigned int)) 138 139 /* Calculate 1-D index a[ij] corresponding to 2-D array a[row][col], with 64-bit precision. 140 * Use gmt_M_ijp when array is padded by BC rows/cols, else use gmt_M_ij0. In all cases 141 * we pass the interior row,col as padding is added by the macro. Note that row,col may 142 * be negative as we seek to address nodes within the padding itself. Hence calculations 143 * use int64_t for signed integers, but cast final index to uint64_t. Finally, there is 144 * gmt_M_ijpgi which is gmt_M_ijp for when there are more than 1 band (it uses h->n_bands). */ 145 146 /*! IJP macro using h and the pad info */ 147 #define gmt_M_ijp(h,row,col) ((uint64_t)(((int64_t)(row)+(int64_t)h->pad[YHI])*((int64_t)h->mx)+(int64_t)(col)+(int64_t)h->pad[XLO])) 148 /*! IJ0 macro using h but ignores the pad info */ 149 #define gmt_M_ij0(h,row,col) ((uint64_t)(((int64_t)(row))*((int64_t)h->n_columns)+(int64_t)(col))) 150 /*! IJ macro using h but treats the entire grid with pad as no-pad grid, i.e. using mx as width */ 151 #define gmt_M_ij(h,row,col) ((uint64_t)(((int64_t)(row))*((int64_t)h->mx)+(int64_t)(col))) 152 /*! IJPGI macro using h and the pad info that works for either grids (n_bands = 1) or images (n_bands = 1,3,4) */ 153 #define gmt_M_ijpgi(h,row,col) ((uint64_t)(((int64_t)(row)+(int64_t)h->pad[YHI])*((int64_t)h->mx*(int64_t)h->n_bands)+((int64_t)(col)+(int64_t)h->pad[XLO])*(int64_t)h->n_bands)) 154 155 /*! Obtain row and col from index */ 156 #define gmt_M_col(h,ij) (((ij) % h->mx) - h->pad[XLO]) 157 #define gmt_M_row(h,ij) (((ij) / h->mx) - h->pad[YHI]) 158 159 /* To set up a standard double for-loop over rows and columns to visit all nodes in a padded array by computing the node index, use gmt_M_grd_loop */ 160 /* Note: All arguments must be actual variables and not expressions. 161 * Note: that input col, row _may_ be signed, hence we do the cast to (int) here. */ 162 163 #define gmt_M_row_loop(C,G,row) for (row = 0; row < (int)G->header->n_rows; row++) 164 #define gmt_M_col_loop(C,G,row,col,ij) for (col = 0, ij = gmt_M_ijp (G->header, row, 0); col < (int)G->header->n_columns; col++, ij++) 165 #define gmt_M_grd_loop(C,G,row,col,ij) gmt_M_row_loop(C,G,row) gmt_M_col_loop(C,G,row,col,ij) 166 /*! Just a loop over columns */ 167 #define gmt_M_col_loop2(C,G,col) for (col = 0; col < (int)G->header->n_columns; col++) 168 169 /* The usage could be: 170 gmt_M_grd_loop (GMT, Grid, row, col, node) fprintf (stderr, "Value at row = %d and col = %d is %g\n", row, col, Grid->data[node]); 171 */ 172 /* The gmt_M_y_is_outside macro returns true if y is outside the given domain. 173 * For gmt_x_is_outside, see the function in gmt_support.c since we must also deal with longitude periodicity. 174 */ 175 176 /*! gmt_M_is_subset is true if wesn is set and wesn cuts through the grid region */ 177 #define gmt_M_is_subset(C,h,R) (R[XHI] > R[XLO] && R[YHI] > R[YLO] && (R[XLO] > h->wesn[XLO] || R[XHI] < h->wesn[XHI] || R[YLO] > h->wesn[YLO] || R[YHI] < h->wesn[YHI])) 178 /*! gmt_M_grd_same_region is true if two grids have the exact same regions */ 179 #define gmt_M_grd_same_region(C,G1,G2) (G1->header->wesn[XLO] == G2->header->wesn[XLO] && G1->header->wesn[XHI] == G2->header->wesn[XHI] && G1->header->wesn[YLO] == G2->header->wesn[YLO] && G1->header->wesn[YHI] == G2->header->wesn[YHI]) 180 /*! gmt_M_grd_same_inc is true if two grids have the exact same grid increments */ 181 #define gmt_M_grd_same_inc(C,G1,G2) (G1->header->inc[GMT_X] == G2->header->inc[GMT_X] && G1->header->inc[GMT_Y] == G2->header->inc[GMT_Y]) 182 /*! gmt_M_grd_equal_inc is true if a grid has approximately the same x and y grid increments (within 1e-6) */ 183 #define gmt_M_grd_equal_xy_inc(C,G) (fabs (1.0 - G->header->inc[GMT_X] / G->header->inc[GMT_Y]) < 1.0e-6) 184 /*! GMT_grd_same_dim is true if two grids have the exact same dimensions and registrations */ 185 #define gmt_M_grd_same_shape(C,G1,G2) (G1->header->n_columns == G2->header->n_columns && G1->header->n_rows == G2->header->n_rows && G1->header->registration == G2->header->registration) 186 /*! gmt_M_grd_same_pad is true if two grids have the exact same pad */ 187 #define gmt_M_grd_same_pad(C,G1,G2) (G1->header->pad[XLO] == G2->header->pad[XLO] && G1->header->pad[XHI] == G2->header->pad[XHI] && G1->header->pad[YLO] == G2->header->pad[YLO] && G1->header->pad[YHI] == G2->header->pad[YHI]) 188 /*! gmt_M_y_is_outside is true if y is outside the given range */ 189 #define gmt_M_y_is_outside(C,y,bottom,top) ((gmt_M_is_dnan(y) || (y) < bottom || (y) > top) ? true : false) 190 191 /*! gmt_M_grd_duplicate_column is true for geographical global grid where first and last data columns are identical */ 192 #define gmt_M_grd_duplicate_column(C,h,way) (C->current.io.col_type[way][GMT_X] == GMT_IS_LON && gmt_M_360_range (h->wesn[XHI], h->wesn[XLO]) && h->registration == GMT_GRID_NODE_REG) 193 194 #endif /* GMT_GRID_H */ 195