/*-------------------------------------------------------------------- * * Copyright (c) 1991-2020 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 _ N C . C R O U T I N E S * * Takes care of all grd input/output built on NCAR's NetCDF routines. * This version is intended to provide more general support for reading * NetCDF files that were not generated by GMT. At the same time, the grids * written by these routines are intended to be more conform COARDS conventions. * These routines are to eventually replace the older gmt_cdf_ routines. * * Most functions will return with error message if an internal error is returned. * There functions are only called indirectly via the GMT_* grdio functions. * * Author: Remko Scharroo * Date: 04-AUG-2005 * Version: 1 * * Added support for chunked I/O, Florian Wobbe, June 2012. * * Functions include: * * gmt_nc_read_grd_info: Read header from file * gmt_nc_read_grd: Read data set from file * gmt_nc_update_grd_info: Update header in existing file * gmt_nc_write_grd_info: Write header to new file * gmt_nc_write_grd: Write header and data set to new file * gmtlib_is_nc_grid: Determine if we have a nc grid * * Private functions: * gmtnc_setup_chunk_cache: Change the default HDF5 chunk cache settings * gmtnc_pad_grid: Add padding to a grid * gmtnc_unpad_grid: Remove padding from a grid * gmtnc_padding_copy: Fill padding by replicating the border cells * gmtnc_padding_zero: Fill padding with zeros * gmtnc_n_chunked_rows_in_cache Determines how many chunks to read at once * gmtnc_io_nc_grid Does the actual netcdf I/O * gmtnc_netcdf_libvers returns the netCDF library version * gmtnc_right_shift_grid * gmtnc_set_optimal_chunksize Determines the optimal chunksize * gmtnc_get_units * gmtnc_put_units * gmtnc_check_step * gmtnc_grd_info * gmtnc_grid_fix_repeat_col * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ #include "gmt_dev.h" #include "gmt_internals.h" #include /* Declaration modifier for netcdf DLL support * annoying: why can't netcdf.h do this on its own? */ #if defined WIN32 && ! defined NETCDF_STATIC #define DLL_NETCDF #endif /* HDF5 chunk cache: reasonable defaults assuming min. chunk size of 128x128 and type byte */ #define NC_CACHE_SIZE 33554432 /* 32MiB */ #define NC_CACHE_NELEMS 2053 /* prime > NC_CACHE_SIZE / (128*128*1byte) */ #define NC_CACHE_PREEMPTION 0.75f int gmt_cdf_grd_info (struct GMT_CTRL *GMT, int ncid, struct GMT_GRID_HEADER *header, char job); int gmt_cdf_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode); static int nc_libvers[] = {-1, -1, -1, -1}; /* holds the version of the netCDF library */ static char *regtype[2] = {"gridline", "pixel"}; GMT_LOCAL const int * gmtnc_netcdf_libvers (void) { static bool inquired = false; if (!inquired) { const char *vers_string = nc_inq_libvers(); sscanf (vers_string, "%d.%d.%d.%d", nc_libvers, nc_libvers+1, nc_libvers+2, nc_libvers+3); inquired = true; } return nc_libvers; /* return pointer to version array */ } /* netcdf I/O mode */ enum Netcdf_io_mode { k_put_netcdf = 0, k_get_netcdf }; /* Wrapper around gmt_nc_put_vara_grdfloat and gmt_nc_get_vara_grdfloat */ static inline int gmtnc_vara_grdfloat (int ncid, int varid, const size_t *startp, const size_t *countp, gmt_grdfloat *fp, unsigned io_mode) { if (io_mode == k_put_netcdf) /* write netcdf */ return gmt_nc_put_vara_grdfloat (ncid, varid, startp, countp, fp); /* read netcdf */ return gmt_nc_get_vara_grdfloat (ncid, varid, startp, countp, fp); } /* Wrapper around gmt_nc_put_varm_grdfloat and gmt_nc_get_varm_grdfloat */ static inline int gmtnc_varm_grdfloat (int ncid, int varid, const size_t *startp, const size_t *countp, const ptrdiff_t *stridep, const ptrdiff_t *imapp, gmt_grdfloat *fp, unsigned io_mode) { if (io_mode == k_put_netcdf) /* write netcdf */ return gmt_nc_put_varm_grdfloat (ncid, varid, startp, countp, stridep, imapp, fp); /* read netcdf */ return gmt_nc_get_varm_grdfloat (ncid, varid, startp, countp, stridep, imapp, fp); } /* Get number of chunked rows that fit into cache (32MiB) */ GMT_LOCAL int gmtnc_n_chunked_rows_in_cache (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, unsigned width, unsigned height, size_t *n_contiguous_chunk_rows, size_t *chunksize) { struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header); nc_type z_type; /* type of z variable */ size_t z_size; /* size of z variable */ size_t z_bytes; /* Number of bytes */ size_t width_t = (size_t)width; size_t height_t = (size_t)height; unsigned yx_dim[2] = {HH->xy_dim[1], HH->xy_dim[0]}; /* because xy_dim not row major */ int err, storage_in; gmt_M_err_trap (nc_inq_vartype(HH->ncid, HH->z_id, &z_type)); /* type of z */ gmt_M_err_trap (nc_inq_type(HH->ncid, z_type, NULL, &z_size)); /* size of z elements in bytes */ gmt_M_err_trap (nc_inq_var_chunking (HH->ncid, HH->z_id, &storage_in, chunksize)); if (storage_in != NC_CHUNKED) { /* default if NC_CONTIGUOUS */ chunksize[yx_dim[0]] = 128; /* 128 rows */ chunksize[yx_dim[1]] = width_t; /* all columns */ } z_bytes = height_t * height_t * z_size; if (z_bytes > NC_CACHE_SIZE) { /* memory needed for subset exceeds the cache size */ unsigned int level; size_t chunks_per_row = (size_t) ceil ((double)width / chunksize[yx_dim[1]]); *n_contiguous_chunk_rows = NC_CACHE_SIZE / (width_t * z_size) / chunksize[yx_dim[0]]; #ifdef NC4_DEBUG level = GMT_MSG_WARNING; #else level = GMT_MSG_DEBUG; #endif GMT_Report (GMT->parent, level, "processing at most %" PRIuS " (%" PRIuS "x%" PRIuS ") chunks at a time (%.1lf MiB)...\n", *n_contiguous_chunk_rows * chunks_per_row, *n_contiguous_chunk_rows, chunks_per_row, *n_contiguous_chunk_rows * z_size * width * chunksize[yx_dim[0]] / 1048576.0); } else *n_contiguous_chunk_rows = 0; /* all chunks fit into cache */ return GMT_NOERROR; } /* Read and write classic or chunked netcdf files */ GMT_LOCAL int gmtnc_io_nc_grid (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, unsigned dim[], unsigned origin[], size_t stride, unsigned io_mode, gmt_grdfloat* grid) { /* io_mode = k_get_netcdf: read a netcdf file to grid * io_mode = k_put_netcdf: write a grid to netcdf */ int status = NC_NOERR; unsigned width = dim[1], height = dim[0]; unsigned yx_dim[2]; /* because xy_dim is not row major! */ size_t width_t = (size_t)width; size_t height_t = (size_t)height; size_t chunksize[5]; /* chunksize of z */ size_t start[5] = {0,0,0,0,0}, count[5] = {1,1,1,1,1}; size_t n_contiguous_chunk_rows = 0; /* that are processed at once, 0 = all */ ptrdiff_t imap[5] = {1,1,1,1,1}; /* mapping between dims of netCDF and in-memory grid */ const ptrdiff_t onestride[5] = {1,1,1,1,1}; /* Passing this instead of NULL bypasses netCDF bug in 4.6.2 */ struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header); /* catch illegal io_mode in debug */ assert (io_mode == k_put_netcdf || io_mode == k_get_netcdf); #ifdef NC4_DEBUG GMT_Report (GMT->parent, GMT_MSG_WARNING, "%s n_columns:%u n_rows:%u x0:%u y0:%u y-order:%s\n", io_mode == k_put_netcdf ? "writing," : "reading,", dim[1], dim[0], origin[1], origin[0], HH->row_order == k_nc_start_south ? "S->N" : "N->S"); #endif /* set index of input origin */ yx_dim[0] = HH->xy_dim[1], yx_dim[1] = HH->xy_dim[0]; /* xy_dim not row major */ memcpy (start, HH->t_index, 3 * sizeof(size_t)); /* set lower dimensions first (e.g. layer) */ start[yx_dim[0]] = origin[0]; /* first row */ start[yx_dim[1]] = origin[1]; /* first col */ /* set mapping of complex grids or if reading a part of a grid */ imap[yx_dim[0]] = (stride == 0 ? width : stride); /* distance between each row */ imap[yx_dim[1]] = 1; /* distance between values in a row */ /* determine how many chunks to process at once */ gmtnc_n_chunked_rows_in_cache (GMT, header, width, height, &n_contiguous_chunk_rows, chunksize); if (n_contiguous_chunk_rows) { /* read/write grid in chunks to keep memory footprint low */ size_t remainder; #ifdef NC4_DEBUG unsigned row_num = 0; GMT_Report (GMT->parent, GMT_MSG_WARNING, "stride: %u width: %u\n", stride, width); #endif /* adjust row count, so that it ends on the bottom of a chunk */ count[yx_dim[0]] = chunksize[yx_dim[0]] * n_contiguous_chunk_rows; remainder = start[yx_dim[0]] % chunksize[yx_dim[0]]; count[yx_dim[0]] -= remainder; count[yx_dim[1]] = width_t; while ( (start[yx_dim[0]] + count[yx_dim[0]]) <= height_t && status == NC_NOERR) { #ifdef NC4_DEBUG GMT_Report (GMT->parent, GMT_MSG_WARNING, "chunked row #%u start-y:%" PRIuS " height:%" PRIuS "\n", ++row_num, start[yx_dim[0]], count[yx_dim[0]]); #endif /* get/put chunked rows */ if (stride) status = gmtnc_varm_grdfloat (HH->ncid, HH->z_id, start, count, onestride, imap, grid, io_mode); else status = gmtnc_vara_grdfloat (HH->ncid, HH->z_id, start, count, grid, io_mode); /* advance grid location and set new origin */ grid += count[yx_dim[0]] * ((stride == 0 ? width_t : stride)); start[yx_dim[0]] += count[yx_dim[0]]; if (remainder) { /* reset count to full chunk height */ count[yx_dim[0]] += remainder; remainder = 0; } } if ( start[yx_dim[0]] != height_t && status == NC_NOERR ) { /* get/put last chunked row */ count[yx_dim[0]] = height_t - start[yx_dim[0]] + origin[0]; #ifdef NC4_DEBUG GMT_Report (GMT->parent, GMT_MSG_WARNING, "chunked row #%u start-y:%" PRIuS " height:%" PRIuS "\n", ++row_num, start[yx_dim[0]], count[yx_dim[0]]); #endif if (stride) status = gmtnc_varm_grdfloat (HH->ncid, HH->z_id, start, count, onestride, imap, grid, io_mode); else status = gmtnc_vara_grdfloat (HH->ncid, HH->z_id, start, count, grid, io_mode); } } else { /* get/put whole grid contiguous */ count[yx_dim[0]] = height_t; count[yx_dim[1]] = width_t; if (stride) status = gmtnc_varm_grdfloat (HH->ncid, HH->z_id, start, count, onestride, imap, grid, io_mode); else status = gmtnc_vara_grdfloat (HH->ncid, HH->z_id, start, count, grid, io_mode); } return status; } GMT_LOCAL void gmtnc_get_units (struct GMT_CTRL *GMT, int ncid, int varid, char *name_units) { /* Get attributes long_name and units for given variable ID * and assign variable name if attributes are not available. * ncid, varid : as in nc_get_att_text * nameunit : long_name and units in form "long_name [units]" */ char name[GMT_GRID_UNIT_LEN80], units[GMT_GRID_UNIT_LEN80]; if (gmtlib_nc_get_att_text (GMT, ncid, varid, "long_name", name, GMT_GRID_UNIT_LEN80)) nc_inq_varname (ncid, varid, name); if (!gmtlib_nc_get_att_text (GMT, ncid, varid, "units", units, GMT_GRID_UNIT_LEN80) && units[0]) snprintf (name_units, GMT_GRID_UNIT_LEN80, "%s [%s]", name, units); else strncpy (name_units, name, GMT_GRID_UNIT_LEN80-1); } GMT_LOCAL void gmtnc_put_units (int ncid, int varid, char *name_units) { /* Put attributes long_name and units for given variable ID based on * string name_unit in the form "long_name [units]". * ncid, varid : as is nc_put_att_text * name_units : string in form "long_name [units]" */ bool copy = false; int i = 0, j = 0; char name[GMT_GRID_UNIT_LEN80], units[GMT_GRID_UNIT_LEN80]; strncpy (name, name_units, GMT_GRID_UNIT_LEN80-1); units[0] = '\0'; for (i = 0; i < GMT_GRID_UNIT_LEN80 && name[i]; i++) { if (name[i] == ']') copy = false, units[j] = '\0'; if (copy) units[j++] = name[i]; if (name[i] == '[') { name[i] = '\0'; if (i > 0 && name[i-1] == ' ') name[i-1] = '\0'; copy = true; } } if (name[0]) nc_put_att_text (ncid, varid, "long_name", strlen(name), name); if (units[0]) nc_put_att_text (ncid, varid, "units", strlen(units), units); if (strstr(units, "degrees_east")) { nc_put_att_text (ncid, varid, "standard_name", 9, "longitude"); nc_put_att_text (ncid, varid, "axis", 1, "X"); } else if (strstr(units, "degrees_north")) { nc_put_att_text (ncid, varid, "standard_name", 8, "latitude"); nc_put_att_text (ncid, varid, "axis", 1, "Y"); } } GMT_LOCAL int gmtnc_check_step (struct GMT_CTRL *GMT, uint32_t n, double *x, char *varname, char *file, bool save_xy_array) { /* Check if all steps in range are the same (within 0.1%). Returns 0 if OK and 1 if variable spacing */ double step, step_min, step_max; uint32_t i; unsigned int wlevel = (save_xy_array) ? GMT_MSG_INFORMATION : GMT_MSG_WARNING; if (n < 2) return 0; step_min = step_max = x[1]-x[0]; for (i = 2; i < n; i++) { step = x[i]-x[i-1]; if (step < step_min) step_min = step; if (step > step_max) step_max = step; } if (fabs (step_min-step_max)/(fabs (step_min)*0.5 + fabs (step_max)*0.5) > 0.001) { GMT_Report (GMT->parent, wlevel, "The step size of coordinate (%s) in grid %s is not constant.\n", varname, file); GMT_Report (GMT->parent, wlevel, "GMT will use a constant step size of %g; the original ranges from %g to %g.\n", (x[n-1]-x[0])/(n-1), step_min, step_max); return 1; } return 0; } GMT_LOCAL void gmtnc_set_optimal_chunksize (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) { /* For optimal performance, set the number of elements in a given chunk * dimension (n) to be the ceiling of the number of elements in that * dimension of the array variable (d) divided by a natural number N>1. * That is, set n = ceil (d / N). Using a chunk size slightly larger than * this value is also acceptable. For example: 129 = ceil (257 / 2). * Do NOT set n = floor (d / N), for example 128 = floor (257 / 2). */ double chunksize[2] = {128, 128}; /* default min chunksize */ const size_t min_chunk_pixels = (size_t)(chunksize[0] * chunksize[1]); /* min pixel count per chunk */ if (GMT->current.setting.io_nc4_chunksize[0] == k_netcdf_io_classic) /* no chunking with classic model */ return; if (GMT->current.setting.io_nc4_chunksize[0] != k_netcdf_io_chunked_auto && ( header->n_rows >= GMT->current.setting.io_nc4_chunksize[0] || header->n_columns >= GMT->current.setting.io_nc4_chunksize[1])) { /* if chunk size is smaller than grid size */ return; } /* here, chunk size is either k_netcdf_io_chunked_auto or the chunk size is * larger than grid size */ if ( header->nm < min_chunk_pixels ) { /* the grid dimension is too small for chunking to make sense. switch to * classic model */ GMT->current.setting.io_nc4_chunksize[0] = k_netcdf_io_classic; return; } /* adjust default chunk sizes for grids that have more than min_chunk_pixels * cells but less than chunksize (default 128) cells in one dimension */ if (header->n_rows < chunksize[0]) { chunksize[0] = header->n_rows; chunksize[1] = floor (min_chunk_pixels / chunksize[0]); } else if (header->n_columns < chunksize[1]) { chunksize[1] = header->n_columns; chunksize[0] = floor (min_chunk_pixels / chunksize[1]); } /* determine optimal chunk size in the range [chunksize,2*chunksize) */ GMT->current.setting.io_nc4_chunksize[0] = (size_t) ceil (header->n_rows / floor (header->n_rows / chunksize[0])); GMT->current.setting.io_nc4_chunksize[1] = (size_t) ceil (header->n_columns / floor (header->n_columns / chunksize[1])); } GMT_LOCAL bool gmtnc_not_obviously_global (double *we) { /* If range is not 360 and boundaries are not 0, 180, 360 then we pass */ if (!gmt_M_360_range (we[0], we[1])) return true; if (!(doubleAlmostEqualZero (we[0], 0.0) || doubleAlmostEqual (we[0], -180.0))) return true; return false; } GMT_LOCAL bool gmtnc_not_obviously_polar (double *se) { /* If range is not 180 and boundaries are not -90/90 then we pass */ if (!gmt_M_180_range (se[0], se[1])) return true; if (!(doubleAlmostEqual (se[0], -90.0) && doubleAlmostEqual (se[1], 90.0))) return true; return false; } GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, char job) { int j, err, has_vector, has_range, registration, var_spacing = 0; int old_fill_mode, status; double dummy[2], *xy = NULL; char dimname[GMT_GRID_UNIT_LEN80], coord[GMT_LEN8]; nc_type z_type; bool save_xy_array = !strncmp (GMT->init.module_name, "grd2xyz", 7U); struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header); /* Dimension ids, variable ids, etc.. */ int i, ncid, z_id = -1, ids[5] = {-1,-1,-1,-1,-1}, gm_id = -1, dims[5], nvars, ndims = 0; size_t lens[5], item[2]; /* If not yet determined, attempt to get the layer IDs from the variable name */ int n_t_index_by_value = 0; double t_value[3]; if (HH->varname[0]) { char *c = strpbrk (HH->varname, "(["); /* find first occurrence of ( or [ */ if (c != NULL && *c == '(') { n_t_index_by_value = sscanf (c+1, "%lf,%lf,%lf", &t_value[0], &t_value[1], &t_value[2]); *c = '\0'; /* t_index will be determined later from t_value when the nc-file is opened */ } else if (c != NULL && *c == '[') { sscanf (c+1, "%" SCNuS ",%" SCNuS ",%" SCNuS "", &HH->t_index[0], &HH->t_index[1], &HH->t_index[2]); *c = '\0'; } } /* Open NetCDF file */ if (!strcmp (HH->name,"=")) return (GMT_GRDIO_NC_NO_PIPE); switch (job) { case 'r': gmt_M_err_trap (nc_open (HH->name, NC_NOWRITE, &ncid)); break; case 'u': gmt_M_err_trap (nc_open (HH->name, NC_WRITE, &ncid)); gmt_M_err_trap (nc_set_fill (ncid, NC_NOFILL, &old_fill_mode)); break; default: /* create new nc-file */ gmtnc_set_optimal_chunksize (GMT, header); if (GMT->current.setting.io_nc4_chunksize[0] != k_netcdf_io_classic) { /* create chunked nc4 file */ gmt_M_err_trap (nc_create (HH->name, NC_NETCDF4, &ncid)); HH->is_netcdf4 = true; } else { /* create nc using classic data model */ gmt_M_err_trap (nc_create (HH->name, NC_CLOBBER, &ncid)); HH->is_netcdf4 = false; } gmt_M_err_trap (nc_set_fill (ncid, NC_NOFILL, &old_fill_mode)); break; } /* Retrieve or define dimensions and variables */ if (job == 'r' || job == 'u') { int kind; /* determine netCDF data model */ gmt_M_err_trap (nc_inq_format(ncid, &kind)); HH->is_netcdf4 = (kind == NC_FORMAT_NETCDF4 || kind == NC_FORMAT_NETCDF4_CLASSIC); /* First see if this is an old NetCDF formatted file */ if (!nc_inq_dimid (ncid, "xysize", &i)) return (gmt_cdf_grd_info (GMT, ncid, header, job)); /* Find first 2-dimensional (z) variable or specified variable */ if (!HH->varname[0]) { /* No specific named variable was requested; search for suitable matrix */ int ID = -1, dim = 0; gmt_M_err_trap (nc_inq_nvars (ncid, &nvars)); i = 0; while (i < nvars && z_id < 0) { /* Look for first 2D grid, with fallback to first higher-dimension (3-4D) grid if 2D not found */ gmt_M_err_trap (nc_inq_varndims (ncid, i, &ndims)); if (ndims == 2) /* Found the first 2-D grid */ z_id = i; else if (ID == -1 && ndims > 2 && ndims < 5) { /* Also look for higher-dim grid in case no 2-D */ ID = i; dim = ndims; } i++; } if (z_id < 0) { /* No 2-D grid found, check if we found a higher dimension cube */ if (ID == -1) return (GMT_GRDIO_NO_2DVAR); /* No we didn't */ z_id = ID; /* Pick the higher dimensioned cube instead, get its name, and warn */ nc_inq_varname (ncid, z_id, HH->varname); GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No 2-D array in file %s. Selecting first 2-D slice in the %d-D array %s\n", HH->name, dim, HH->varname); } } else if (nc_inq_varid (ncid, HH->varname, &z_id) == NC_NOERR) { /* Gave a named variable that exist, is it 2-4 D? */ gmt_M_err_trap (nc_inq_varndims (ncid, z_id, &ndims)); if (ndims < 2 || ndims > 5) return (GMT_GRDIO_BAD_DIM); } else /* No such named variable in the grid */ return (GMT_GRDIO_NO_VAR); if (z_id < 0) return (GMT_GRDIO_NO_2DVAR); /* Get the z data type and determine its dimensions */ gmt_M_err_trap (nc_inq_vartype (ncid, z_id, &z_type)); gmt_M_err_trap (nc_inq_vardimid (ncid, z_id, dims)); switch (z_type) { case NC_UBYTE: header->type = GMT_GRID_IS_NB; HH->orig_datatype = GMT_UCHAR; break; case NC_BYTE: header->type = GMT_GRID_IS_NB; HH->orig_datatype = GMT_CHAR; break; case NC_SHORT: header->type = GMT_GRID_IS_NS; HH->orig_datatype = GMT_SHORT; break; case NC_INT: header->type = GMT_GRID_IS_NI; HH->orig_datatype = GMT_INT; break; case NC_FLOAT: header->type = GMT_GRID_IS_NF; HH->orig_datatype = GMT_FLOAT; break; case NC_DOUBLE: header->type = GMT_GRID_IS_ND; HH->orig_datatype = GMT_DOUBLE; break; default: header->type = k_grd_unknown_fmt; break; } /* Get the ids of the x and y (and depth and time) coordinate variables */ for (i = 0; i < ndims; i++) { gmt_M_err_trap (nc_inq_dim (ncid, dims[i], dimname, &lens[i])); //if (nc_inq_varid (ncid, dimname, &ids[i])) return (GMT_GRDIO_NC_NOT_COARDS); if ((status = nc_inq_varid (ncid, dimname, &ids[i])) != NC_NOERR) GMT_Report (GMT->parent, GMT_MSG_WARNING, "\"%s\", %s\n\tIf something bad happens later, try importing via GDAL.\n", dimname, nc_strerror(status)); } HH->xy_dim[0] = ndims-1; HH->xy_dim[1] = ndims-2; /* Check if LatLon variable exists, then we may need to flip x and y */ if (nc_inq_varid (ncid, "LatLon", &i) == NC_NOERR) nc_get_var_int (ncid, i, HH->xy_dim); header->n_columns = (uint32_t) lens[HH->xy_dim[0]]; header->n_rows = (uint32_t) lens[HH->xy_dim[1]]; /* Check if the grid_mapping variable exists */ if (nc_inq_varid (ncid, "grid_mapping", &i) == NC_NOERR) gm_id = i; } /* if (job == 'r' || job == 'u') */ else { /* Define dimensions of z variable */ ndims = 2; HH->xy_dim[0] = 1; HH->xy_dim[1] = 0; strcpy (coord, (gmt_M_x_is_lon (GMT, GMT_OUT)) ? "lon" : (gmt_M_type (GMT, GMT_OUT, GMT_X) & GMT_IS_RATIME) ? "time" : "x"); gmt_M_err_trap (nc_def_dim (ncid, coord, (size_t) header->n_columns, &dims[1])); gmt_M_err_trap (nc_def_var (ncid, coord, NC_DOUBLE, 1, &dims[1], &ids[1])); strcpy (coord, (gmt_M_y_is_lat (GMT, GMT_OUT)) ? "lat" : (gmt_M_type (GMT, GMT_OUT, GMT_Y) & GMT_IS_RATIME) ? "time" : "y"); gmt_M_err_trap (nc_def_dim (ncid, coord, (size_t) header->n_rows, &dims[0])); gmt_M_err_trap (nc_def_var (ncid, coord, NC_DOUBLE, 1, &dims[0], &ids[0])); switch (header->type) { case GMT_GRID_IS_NB: z_type = NC_BYTE; break; case GMT_GRID_IS_NS: z_type = NC_SHORT; break; case GMT_GRID_IS_NI: z_type = NC_INT; break; case GMT_GRID_IS_NF: z_type = NC_FLOAT; break; case GMT_GRID_IS_ND: z_type = NC_DOUBLE; break; default: z_type = NC_NAT; } /* Variable name is given, or defaults to "z" */ if (!HH->varname[0]) strcpy (HH->varname, "z"); gmt_M_err_trap (nc_def_var (ncid, HH->varname, z_type, 2, dims, &z_id)); /* set deflation and chunking */ if (GMT->current.setting.io_nc4_chunksize[0] != k_netcdf_io_classic) { /* set chunk size */ gmt_M_err_trap (nc_def_var_chunking (ncid, z_id, NC_CHUNKED, GMT->current.setting.io_nc4_chunksize)); /* set deflation level and shuffle for x, y, and z variable */ if (GMT->current.setting.io_nc4_deflation_level) { gmt_M_err_trap (nc_def_var_deflate (ncid, ids[1], true, true, GMT->current.setting.io_nc4_deflation_level)); gmt_M_err_trap (nc_def_var_deflate (ncid, ids[0], true, true, GMT->current.setting.io_nc4_deflation_level)); gmt_M_err_trap (nc_def_var_deflate (ncid, z_id, true, true, GMT->current.setting.io_nc4_deflation_level)); } } /* GMT->current.setting.io_nc4_chunksize[0] != k_netcdf_io_classic */ } /* if (job == 'r' || job == 'u') */ HH->z_id = z_id; HH->ncid = ncid; /* Query or assign attributes */ if (job == 'u') gmt_M_err_trap (nc_redef (ncid)); if (job == 'r') { bool set_reg = true; double dx = 0, dy = 0; /* Get global information */ if (gmtlib_nc_get_att_text (GMT, ncid, NC_GLOBAL, "title", header->title, GMT_GRID_TITLE_LEN80)) gmtlib_nc_get_att_text (GMT, ncid, z_id, "long_name", header->title, GMT_GRID_TITLE_LEN80); if (gmtlib_nc_get_att_text (GMT, ncid, NC_GLOBAL, "history", header->command, GMT_GRID_COMMAND_LEN320)) gmtlib_nc_get_att_text (GMT, ncid, NC_GLOBAL, "source", header->command, GMT_GRID_COMMAND_LEN320); gmtlib_nc_get_att_text (GMT, ncid, NC_GLOBAL, "description", header->remark, GMT_GRID_REMARK_LEN160); header->registration = GMT_GRID_NODE_REG; if (!nc_get_att_int (ncid, NC_GLOBAL, "node_offset", &i)) { /* GMT wrote the registration in the grid */ header->registration = i; set_reg = false; /* Do not update it below since we know the registration */ } if (gm_id > 0) { size_t len; char *pch = NULL; gmt_M_err_trap (nc_inq_attlen (ncid, gm_id, "spatial_ref", &len)); /* Get attrib length */ gmt_M_str_free (header->ProjRefWKT); /* Make sure we didn't have a previously allocated one */ pch = gmt_M_memory (GMT, NULL, len+1, char); /* and allocate the needed space */ gmt_M_err_trap (nc_get_att_text (ncid, gm_id, "spatial_ref", pch)); header->ProjRefWKT = strdup (pch); /* Turn it into a strdup allocation to be compatible with other instances elsewhere */ gmt_M_free (GMT, pch); } /* Explanation for the logic below: Not all netCDF grids are proper COARDS grids and hence we sometime must guess * regarding the settings. The x and y coordinates may be written as arrays, which reflect the positions of the * nodes. There may also be the attributes actual_range which specifies the x range. These will differ if the * grid is pixel registered, hence when both are present we use this difference to detect a pixel grid. However, * some external tools such as xarray may slice a grid but not update the attributes. In this case the actual_range * may have an initial range that is no longer the case. We have added a check if these differ by more than a * half grid increment. If not then we can trust it. If actual_range is missing then we have to guess the registration * which we do by checking if the start coordinate is an integer multiple of the increment. If not, we guess pixel registration * but cannot know if this is the case unless the adjusted coordinates in x has a range of 360 and in y a range of 180. * Finally, if there is no array just the actual_range, then we cannot tell the registration from the range but try * and leave it as gridline registration. */ /* Create enough memory to store the x- and y-coordinate values */ xy = gmt_M_memory (GMT, NULL, MAX (header->n_columns, header->n_rows), double); /* Get information about x variable */ gmtnc_get_units (GMT, ncid, ids[HH->xy_dim[0]], header->x_units); /* Set default range to number of nodes in case nothing is found further down */ dummy[0] = 0.0, dummy[1] = (double) header->n_columns-1; /* Default */ registration = GMT_GRID_NODE_REG; /* Look for the x-coordinate vector */ if ((has_vector = !nc_get_var_double (ncid, ids[HH->xy_dim[0]], xy))) { var_spacing = gmtnc_check_step (GMT, header->n_columns, xy, header->x_units, HH->name, save_xy_array); if (save_xy_array && var_spacing) { if (GMT->current.io.nc_xarray) gmt_M_free (GMT, GMT->current.io.nc_xarray); GMT->current.io.nc_xarray = gmt_M_memory (GMT, NULL, header->n_columns, double); gmt_M_memcpy (GMT->current.io.nc_xarray, xy, header->n_columns, double); HH->var_spacing[GMT_X] = var_spacing; } dx = fabs (xy[1] - xy[0]); /* Grid spacing in x */ } /* Look for the x-coordinate range attributes */ has_range = (!nc_get_att_double (ncid, ids[HH->xy_dim[0]], "actual_range", dummy) || !nc_get_att_double (ncid, ids[HH->xy_dim[0]], "valid_range", dummy) || !(nc_get_att_double (ncid, ids[HH->xy_dim[0]], "valid_min", &dummy[0]) + nc_get_att_double (ncid, ids[HH->xy_dim[0]], "valid_max", &dummy[1]))); if (has_vector && has_range) { /* Has both so we can do a basic sanity check */ if (fabs (dummy[0] - xy[0]) > ((0.5+GMT_CONV5_LIMIT) * dx) || fabs (dummy[1] - xy[header->n_columns-1]) > ((0.5+GMT_CONV5_LIMIT) * dx)) { if (gmtnc_not_obviously_global (dummy)) { GMT_Report (GMT->parent, GMT_MSG_WARNING, "The x-coordinates and range attribute are in conflict; must rely on coordinates only\n"); dummy[0] = xy[0], dummy[1] = xy[header->n_columns-1]; has_range = false; /* Since useless information */ /* For registration, we have to assume that the actual range is an integer multiple of increment. * If so, then if the coordinates are off by 0.5*dx then we assume we have pixel registration */ if (set_reg && fabs (fmod (dummy[0], dx)) > ((0.5-GMT_CONV5_LIMIT) * dx)) { /* Pixel registration */ registration = header->registration = GMT_GRID_PIXEL_REG; dummy[0] -= 0.5 * dx; dummy[1] += 0.5 * dx; if (gmt_M_360_range (dummy[0], dummy[1])) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Grid x-coordinates after pixel registration adjustment have exactly 360 range\n"); } } else { /* Got clean global longitude range, stick with that information and ignore xy */ GMT_Report (GMT->parent, GMT_MSG_WARNING, "The x-coordinates and range attribute are in conflict but range is exactly 360; we rely on this range\n"); if (set_reg && (header->n_columns%2) == 0) { /* Pixel registration */ registration = header->registration = GMT_GRID_PIXEL_REG; GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Global longitudes, guessing registration to be %s since nx is odd\n", regtype[header->registration]); } } } else { /* Data seems OK; determine registration */ dummy[0] = xy[0], dummy[1] = xy[header->n_columns-1]; if (set_reg && (fabs(dummy[1] - dummy[0]) / fabs(xy[header->n_columns-1] - xy[0]) - 1.0 > 0.5 / (header->n_columns - 1))) registration = header->registration = GMT_GRID_PIXEL_REG; } } else if (has_vector) { /* No attribute for range, use coordinates */ dummy[0] = xy[0], dummy[1] = xy[header->n_columns-1]; if (set_reg && fabs (fmod (dummy[0], dx)) > ((0.5-GMT_CONV5_LIMIT) * dx)) { /* Most likely pixel registration since off by dx/2 */ registration = header->registration = GMT_GRID_PIXEL_REG; dummy[0] -= 0.5 * dx; dummy[1] += 0.5 * dx; if (gmt_M_360_range (dummy[0], dummy[1])) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Grid x-coordinates after pixel registration adjustment have exactly 360 range\n"); } GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No range attribute, guessing registration to be %s\n", regtype[header->registration]); } else { /* Only has the valid_range settings. If no registration set, and no dx available, guess based on nx */ if (set_reg && (header->n_columns%2) == 0) { /* Pixel registration */ registration = header->registration = GMT_GRID_PIXEL_REG; GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No x-coordinates, guessing registration to be %s since nx is odd\n", regtype[header->registration]); } } /* Determine grid step */ header->inc[GMT_X] = gmt_M_get_inc (GMT, dummy[0], dummy[1], header->n_columns, registration); if (gmt_M_is_dnan(header->inc[GMT_X])) header->inc[GMT_X] = 1.0; #ifdef NC4_DEBUG GMT_Report (GMT->parent, GMT_MSG_WARNING, "x registration: %u\n", header->registration); GMT_Report (GMT->parent, GMT_MSG_WARNING, "x dummy: %g %g\n", dummy[0], dummy[1]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "x[0] x[nx-1]: %g %g\n", xy[0], xy[header->n_columns-1]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "xinc: %g %g\n", header->inc[GMT_X]); #endif /* Extend x boundaries by half if we found pixel registration */ if (registration == GMT_GRID_NODE_REG && header->registration == GMT_GRID_PIXEL_REG) header->wesn[XLO] = dummy[0] - header->inc[GMT_X] / 2.0, header->wesn[XHI] = dummy[1] + header->inc[GMT_X] / 2.0; else /* Use as is */ header->wesn[XLO] = dummy[0], header->wesn[XHI] = dummy[1]; /* Get information about y variable */ gmtnc_get_units (GMT, ncid, ids[HH->xy_dim[1]], header->y_units); /* Set default range to number of nodes in case nothing is found further down */ dummy[0] = 0.0, dummy[1] = (double) header->n_rows-1; /* Default */ registration = GMT_GRID_NODE_REG; /* Read the y-coordinate vector (if available), otherwise just look for range attributes */ if ((has_vector = !nc_get_var_double (ncid, ids[HH->xy_dim[1]], xy))) { var_spacing = gmtnc_check_step (GMT, header->n_rows, xy, header->y_units, HH->name, save_xy_array); if (save_xy_array && var_spacing) { if (GMT->current.io.nc_yarray) gmt_M_free (GMT, GMT->current.io.nc_yarray); GMT->current.io.nc_yarray = gmt_M_memory (GMT, NULL, header->n_rows, double); gmt_M_memcpy (GMT->current.io.nc_yarray, xy, header->n_rows, double); HH->var_spacing[GMT_Y] = var_spacing; /* Flip-ud y-array since row = 0 is the last value */ gmt_grd_flip_vertical (GMT->current.io.nc_yarray, 1, header->n_rows, 0, sizeof(double)); } dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1]; dy = fabs (xy[1] - xy[0]); /* Grid spacing in y */ } has_range = (!nc_get_att_double (ncid, ids[HH->xy_dim[1]], "actual_range", dummy) || !nc_get_att_double (ncid, ids[HH->xy_dim[1]], "valid_range", dummy) || !(nc_get_att_double (ncid, ids[HH->xy_dim[1]], "valid_min", &dummy[0]) + nc_get_att_double (ncid, ids[HH->xy_dim[1]], "valid_max", &dummy[1]))); if (has_vector && has_range) { /* Has both so we can do a basic sanity check */ if (fabs (dummy[0] - xy[0]) > ((0.5+GMT_CONV5_LIMIT) * dy) || fabs (dummy[1] - xy[header->n_rows-1]) > ((0.5+GMT_CONV5_LIMIT) * dy)) { if (gmtnc_not_obviously_polar (dummy)) { GMT_Report (GMT->parent, GMT_MSG_WARNING, "The y-coordinates and range attribute are in conflict; must rely on coordinates only\n"); dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1]; has_range = false; /* Since useless information */ /* Registration was set using x values, so here we just check that we get the same result. * If the coordinates are off by 0.5*dy then we assume we have pixel registration */ if (fabs (fmod (dummy[0], dy)) > ((0.5-GMT_CONV5_LIMIT) * dy)) { /* Pixel registration? */ if (header->registration == GMT_GRID_NODE_REG) /* No, somehow messed up now */ GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]); else { /* Pixel registration confirmed */ dummy[0] -= 0.5 * dy; dummy[1] += 0.5 * dy; registration = GMT_GRID_PIXEL_REG; if (gmt_M_180_range (dummy[0], dummy[1])) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Grid y-coordinates after pixel registration adjustment have exactly 180 range\n"); } } } else { GMT_Report (GMT->parent, GMT_MSG_WARNING, "The y-coordinates and range attribute are in conflict but range is exactly 180; we rely on this range\n"); if ((header->n_rows%2) == 1 && header->registration == GMT_GRID_NODE_REG) /* Pixel registration? */ GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]); else registration = header->registration; } } else { /* Data seems OK; determine registration and set dummy from data coordinates */ dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1]; if ((fabs(dummy[1] - dummy[0]) / fabs(xy[header->n_rows-1] - xy[0]) - 1.0 > 0.5 / (header->n_rows - 1)) && header->registration == GMT_GRID_NODE_REG) GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]); } } else if (has_vector) { /* No attribute for range, use coordinates */ dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1]; if ((fabs (fmod (dummy[0], dy)) > ((0.5-GMT_CONV5_LIMIT) * dy))) { /* Most likely pixel registration since off by dy/2 */ if (header->registration == GMT_GRID_NODE_REG) /* No, somehow messed up now */ GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]); else { /* Pixel registration confirmed */ if (dummy[0] > dummy[1]) { /* Check for reverse order of y-coordinate */ gmt_M_double_swap (dummy[0], dummy[1]); HH->row_order = k_nc_start_north; } dummy[0] -= 0.5 * dy; dummy[1] += 0.5 * dy; registration = GMT_GRID_PIXEL_REG; if (gmt_M_180_range (dummy[0], dummy[1])) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Grid y-coordinates after pixel registration adjustment have exactly 180 range\n"); } } } else { /* Only has the valid_range settings. If no registration set, and no dy available, guess based on ny */ if ((header->n_rows%2) == 1 && header->registration == GMT_GRID_NODE_REG) /* Pixel registration? */ GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]); } /* Check for reverse order of y-coordinate */ if (dummy[0] > dummy[1]) { HH->row_order = k_nc_start_north; gmt_M_double_swap (dummy[0], dummy[1]); } else if (has_vector && (xy[0] > xy[header->n_rows-1]) && (dummy[0] > dummy[1])) /* Here the lat vector is top-down but range is bottum-up */ HH->row_order = k_nc_start_north; else if (HH->row_order != k_nc_start_north) /* If not already set to north in an above test */ HH->row_order = k_nc_start_south; /* Determine grid step */ header->inc[GMT_Y] = gmt_M_get_inc (GMT, dummy[0], dummy[1], header->n_rows, registration); if (gmt_M_is_dnan(header->inc[GMT_Y])) header->inc[GMT_Y] = 1.0; #ifdef NC4_DEBUG GMT_Report (GMT->parent, GMT_MSG_WARNING, "y registration: %u\n", header->registration); GMT_Report (GMT->parent, GMT_MSG_WARNING, "y dummy: %g %g\n", dummy[0], dummy[1]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "y[0] y[ny-1]: %g %g\n", xy[0], xy[header->n_rows-1]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "xinc: %g %g\n", header->inc[GMT_Y]); #endif /* Extend y boundaries by half if we found pixel registration */ if (registration == GMT_GRID_NODE_REG && header->registration == GMT_GRID_PIXEL_REG) header->wesn[YLO] = dummy[0] - header->inc[GMT_Y] / 2.0, header->wesn[YHI] = dummy[1] + header->inc[GMT_Y] / 2.0; else /* Use as is */ header->wesn[YLO] = dummy[0], header->wesn[YHI] = dummy[1]; gmt_M_free (GMT, xy); /* Done with the array */ /* Get information about z variable */ gmtnc_get_units (GMT, ncid, z_id, header->z_units); if (nc_get_att_double (ncid, z_id, "scale_factor", &header->z_scale_factor)) header->z_scale_factor = 1.0; if (nc_get_att_double (ncid, z_id, "add_offset", &header->z_add_offset)) header->z_add_offset = 0.0; #ifdef DOUBLE_PRECISION_GRID if (nc_get_att_double (ncid, z_id, "_FillValue", &header->nan_value)) nc_get_att_double (ncid, z_id, "missing_value", &header->nan_value); #else if (nc_get_att_float (ncid, z_id, "_FillValue", &header->nan_value)) nc_get_att_float (ncid, z_id, "missing_value", &header->nan_value); #endif if (!nc_get_att_double (ncid, z_id, "actual_range", dummy)) { /* z-limits need to be converted from actual to internal grid units. */ header->z_min = (dummy[0] - header->z_add_offset) / header->z_scale_factor; header->z_max = (dummy[1] - header->z_add_offset) / header->z_scale_factor; } else if (!nc_get_att_double (ncid, z_id, "valid_range", dummy) || !(nc_get_att_double (ncid, ids[HH->xy_dim[1]], "valid_min", &dummy[0]) + nc_get_att_double (ncid, ids[HH->xy_dim[1]], "valid_max", &dummy[1]))) { /* Valid range is already in packed units, so do not convert */ header->z_min = dummy[0], header->z_max = dummy[1]; } if (gmt_M_is_dnan (header->z_min) && gmt_M_is_dnan (header->z_max)) { GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "netCDF grid %s information has zmin = zmax = NaN. Reset to 0/0.\n", HH->name); header->z_min = header->z_max = 0.0; } { /* Get deflation and chunking info */ int storage_mode, shuffle = 0, deflate = 0, deflate_level = 0; size_t chunksize[5]; /* chunksize of z */ gmt_M_err_trap (nc_inq_var_chunking (ncid, z_id, &storage_mode, chunksize)); if (storage_mode == NC_CHUNKED) { HH->z_chunksize[0] = chunksize[dims[0]]; /* chunk size in vertical dimension */ HH->z_chunksize[1] = chunksize[dims[1]]; /* chunk size in horizontal dimension */ } else { /* NC_CONTIGUOUS */ HH->z_chunksize[0] = HH->z_chunksize[1] = 0; } if (HH->is_netcdf4) gmt_M_err_trap (nc_inq_var_deflate (ncid, z_id, &shuffle, &deflate, &deflate_level)); HH->z_shuffle = shuffle ? true : false; /* if shuffle filter is turned on */ HH->z_deflate_level = deflate ? deflate_level : 0; /* if deflate filter is in use */ } /* Determine t_index from t_value */ item[0] = 0; for (i = 0; i < n_t_index_by_value; i++) { item[1] = lens[i]-1; if (nc_get_att_double (ncid, ids[i], "actual_range", dummy)) { gmt_M_err_trap (nc_get_var1_double (ncid, ids[i], &item[0], &dummy[0])); gmt_M_err_trap (nc_get_var1_double (ncid, ids[i], &item[1], &dummy[1])); } if (item[1] != 0 && dummy[0] != dummy[1]) { /* avoid dvision by 0 */ double index = (t_value[i] - dummy[0]) / (dummy[1] - dummy[0]) * item[1]; if (index > 0) HH->t_index[i] = (size_t)index; } } #ifdef NC4_DEBUG GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->wesn: %g %g %g %g\n", header->wesn[XLO], header->wesn[XHI], header->wesn[YLO], header->wesn[YHI]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->registration:%u\n", header->registration); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->row_order: %s\n", HH->row_order == k_nc_start_south ? "S->N" : "N->S"); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->n_columns: %3d head->n_rows:%3d\n", header->n_columns, header->n_rows); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->inc: %g %g\n", header->inc[GMT_X], header->inc[GMT_Y]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->mx: %3d head->my:%3d\n", header->mx, header->my); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->nm: %3d head->size:%3d\n", header->nm, header->size); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->t-index %d,%d,%d\n", HH->t_index[0], HH->t_index[1], HH->t_index[2]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->pad xlo:%u xhi:%u ylo:%u yhi:%u\n", header->pad[XLO], header->pad[XHI], header->pad[YLO], header->pad[YHI]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->BC xlo:%u xhi:%u ylo:%u yhi:%u\n", HH->BC[XLO], HH->BC[XHI], HH->BC[YLO], HH->BC[YHI]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->grdtype:%u %u\n", HH->grdtype, GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT); #endif } else { /* Store global attributes */ unsigned int row, col; const int *nc_vers = gmtnc_netcdf_libvers(); GMT_Report (GMT->parent, GMT_MSG_DEBUG, "netCDF Library version: %d\n", *nc_vers); gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "Conventions", strlen(GMT_NC_CONVENTION), GMT_NC_CONVENTION)); if (header->title[0]) gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "title", strlen(header->title), header->title)); gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "history", strlen(header->command), header->command)); if (header->remark[0]) gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "description", strlen(header->remark), header->remark)); gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "GMT_version", strlen(GMT_VERSION), (const char *) GMT_VERSION)); if (header->registration == GMT_GRID_PIXEL_REG) { int reg = header->registration; gmt_M_err_trap (nc_put_att_int (ncid, NC_GLOBAL, "node_offset", NC_LONG, 1U, ®)); } else nc_del_att (ncid, NC_GLOBAL, "node_offset"); #ifdef HAVE_GDAL /* If we have projection information create a container variable named "grid_mapping" with an attribute "spatial_ref" that will hold the projection info in WKT format. GDAL and Mirone know use this info */ if ((header->ProjRefWKT != NULL) || (header->ProjRefPROJ4 != NULL)) { int id[1], dim[1]; if (header->ProjRefWKT == NULL || !header->ProjRefWKT[0]) { /* Must convert from proj4 string to WKT */ OGRSpatialReferenceH hSRS = OSRNewSpatialReference(NULL); if (header->ProjRefPROJ4 && (!strncmp(header->ProjRefPROJ4, "+unavailable", 4) || strlen(header->ProjRefPROJ4) <= 5)) { /* Silently jump out of here */ OSRDestroySpatialReference(hSRS); goto L100; } GMT_Report(GMT->parent, GMT_MSG_INFORMATION, "Proj4 string to be converted to WKT:\n\t%s\n", header->ProjRefPROJ4); if (OSRImportFromProj4(hSRS, header->ProjRefPROJ4) == CE_None) { char *pszPrettyWkt = NULL; OSRExportToPrettyWkt(hSRS, &pszPrettyWkt, false); header->ProjRefWKT = strdup(pszPrettyWkt); CPLFree(pszPrettyWkt); GMT_Report(GMT->parent, GMT_MSG_INFORMATION, "WKT converted from proj4 string:\n%s\n", header->ProjRefWKT); } else { header->ProjRefWKT = NULL; GMT_Report(GMT->parent, GMT_MSG_WARNING, "gmtnc_grd_info failed to convert the proj4 string\n%s\n to WKT\n", header->ProjRefPROJ4); } OSRDestroySpatialReference(hSRS); } if (header->ProjRefWKT != NULL) { /* It may be NULL if the above conversion failed */ if (nc_inq_varid(ncid, "grid_mapping", &id[0]) != NC_NOERR) { gmt_M_err_trap(nc_def_dim(ncid, "grid_mapping", 12U, &dim[0])); gmt_M_err_trap(nc_def_var(ncid, "grid_mapping", NC_CHAR, 1, dim, &id[0])); } gmt_M_err_trap(nc_put_att_text(ncid, id[0], "spatial_ref", strlen(header->ProjRefWKT), header->ProjRefWKT)); gmt_M_err_trap(nc_put_att_text(ncid, z_id, "grid_mapping", 12U, "grid_mapping")); /* Create attrib in z variable */ } } L100: #endif /* Avoid NaN increments */ if (gmt_M_is_dnan(header->inc[GMT_X])) header->inc[GMT_X] = 1.0; if (gmt_M_is_dnan(header->inc[GMT_Y])) header->inc[GMT_Y] = 1.0; /* Define x variable */ gmtnc_put_units (ncid, ids[HH->xy_dim[0]], header->x_units); dummy[0] = header->wesn[XLO], dummy[1] = header->wesn[XHI]; gmt_M_err_trap (nc_put_att_double (ncid, ids[HH->xy_dim[0]], "actual_range", NC_DOUBLE, 2U, dummy)); /* Define y variable */ gmtnc_put_units (ncid, ids[HH->xy_dim[1]], header->y_units); dummy[0] = header->wesn[YLO], dummy[1] = header->wesn[YHI]; gmt_M_err_trap (nc_put_att_double (ncid, ids[HH->xy_dim[1]], "actual_range", NC_DOUBLE, 2U, dummy)); /* When varname is given, and z_units is default, overrule z_units with varname */ if (HH->varname[0] && !strcmp (header->z_units, "z")) strncpy (header->z_units, HH->varname, GMT_GRID_UNIT_LEN80-1); /* Define z variable. Attempt to remove "scale_factor" or "add_offset" when no longer needed */ gmtnc_put_units (ncid, z_id, header->z_units); if (header->z_scale_factor != 1.0) { gmt_M_err_trap (nc_put_att_double (ncid, z_id, "scale_factor", NC_DOUBLE, 1U, &header->z_scale_factor)); } else if (job == 'u') nc_del_att (ncid, z_id, "scale_factor"); if (header->z_add_offset != 0.0) { gmt_M_err_trap (nc_put_att_double (ncid, z_id, "add_offset", NC_DOUBLE, 1U, &header->z_add_offset)); } else if (job == 'u') nc_del_att (ncid, z_id, "add_offset"); if (z_type != NC_FLOAT && z_type != NC_DOUBLE) header->nan_value = rintf (header->nan_value); /* round to integer */ if (job == 'u' && HH->is_netcdf4) { /* netCDF-4 has a bug and crash when * rewriting the _FillValue attribute in netCDF-4 files https://bugtracking.unidata.ucar.edu/browse/NCF-187 To work-around it we implement the renaming trick advised on NCF-133 Edit: This work-around should eventually be removed because the bug was fixed as of 2015-04-02 (any version >4.3.3.1 should be fine). */ gmt_M_err_trap (nc_rename_att (ncid, z_id, "_FillValue", "_fillValue")); #ifdef DOUBLE_PRECISION_GRID gmt_M_err_trap (nc_put_att_double (ncid, z_id, "_fillValue", z_type, 1U, &header->nan_value)); #else gmt_M_err_trap (nc_put_att_float (ncid, z_id, "_fillValue", z_type, 1U, &header->nan_value)); #endif gmt_M_err_trap (nc_rename_att (ncid, z_id, "_fillValue", "_FillValue")); } else { #ifdef DOUBLE_PRECISION_GRID gmt_M_err_trap (nc_put_att_double (ncid, z_id, "_FillValue", z_type, 1U, &header->nan_value)); #else gmt_M_err_trap (nc_put_att_float (ncid, z_id, "_FillValue", z_type, 1U, &header->nan_value)); #endif } /* Limits need to be stored in actual, not internal grid, units */ if (header->z_min <= header->z_max) { dummy[0] = header->z_min * header->z_scale_factor + header->z_add_offset; dummy[1] = header->z_max * header->z_scale_factor + header->z_add_offset; } else dummy[0] = 0.0, dummy[1] = 0.0; gmt_M_err_trap (nc_put_att_double (ncid, z_id, "actual_range", NC_DOUBLE, 2U, dummy)); /* Store values along x and y axes */ nc_inq_nvars (ncid, &nvars); for (j = 0; j < nvars; j++) { gmt_M_err_trap (nc_inq_varndims (ncid, j, &ndims)); } gmt_M_err_trap (nc_enddef (ncid)); xy = gmt_M_memory (GMT, NULL, MAX (header->n_columns,header->n_rows), double); for (col = 0; col < header->n_columns; col++) xy[col] = gmt_M_grd_col_to_x (GMT, col, header); gmt_M_err_trap (nc_put_var_double (ncid, ids[HH->xy_dim[0]], xy)); /* Depending on row_order, write y-coordinate array bottom-to-top or top-to-bottom */ if (HH->row_order == k_nc_start_south) { for (row = 0; row < header->n_rows; row++) xy[row] = (double) gmt_M_col_to_x (GMT, row, header->wesn[YLO], header->wesn[YHI], header->inc[GMT_Y], 0.5 * header->registration, header->n_rows); } else { for (row = 0; row < header->n_rows; row++) xy[row] = (double) gmt_M_row_to_y (GMT, row, header->wesn[YLO], header->wesn[YHI], header->inc[GMT_Y], 0.5 * header->registration, header->n_rows); } gmt_M_err_trap (nc_put_var_double (ncid, ids[HH->xy_dim[1]], xy)); gmt_M_free (GMT, xy); } /* Close NetCDF file, unless job == 'W' */ if (job != 'W') gmt_M_err_trap (nc_close (ncid)); return (GMT_NOERROR); } /* Shift columns in a grid to the right (n_shift < 0) or to the left (n_shift < 0) */ GMT_LOCAL void gmtnc_right_shift_grid (void *gridp, const unsigned n_cols, const unsigned n_rows, int n_shift, size_t cell_size) { char *tmp = NULL, *grid = (char*)gridp; size_t row, n_shift_abs = abs (n_shift), nm; size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows; assert (n_shift_abs != 0 && n_cols_t > n_shift_abs && n_cols_t > 0 && n_rows_t > 0); tmp = calloc (n_shift_abs, cell_size); if (n_shift > 0) { /* right shift */ for (row = 0; row < n_rows_t; ++row) { nm = row * n_cols_t; /* copy last n_shift_abs cols into tmp buffer */ memcpy (tmp, grid + (nm + n_cols_t - n_shift_abs) * cell_size, n_shift_abs * cell_size); /* right shift row */ memmove (grid + (nm + n_shift_abs) * cell_size, grid + nm * cell_size, (n_cols_t - n_shift_abs) * cell_size); /* prepend tmp buffer */ memcpy (grid + nm * cell_size, tmp, n_shift_abs * cell_size); } } else { /* n_shift_abs < 0 */ for (row = 0; row < n_rows_t; ++row) { nm = row * n_cols_t; /* copy first n_shift_abs cols into tmp buffer */ memcpy (tmp, grid + nm * cell_size, n_shift_abs * cell_size); /* left shift row */ memmove (grid + nm * cell_size, grid + (nm + n_shift_abs) * cell_size, (n_cols_t - n_shift_abs) * cell_size); /* append tmp buffer */ memcpy (grid + (nm + n_cols_t - n_shift_abs) * cell_size, tmp, n_shift_abs * cell_size); } } gmt_M_str_free (tmp); } /* Fill padding by replicating the border cells or wrapping around a * row if columns are periodic */ GMT_LOCAL void gmtnc_padding_copy (void *gridp, const unsigned n_cols, const unsigned n_rows, const unsigned *n_pad, size_t cell_size, bool periodic_cols) { /* n_cols and n_rows are dimensions of the padded grid */ char *grid = (char*)gridp; size_t row, cell, nm, nm2; size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows; size_t pad_t[4]; assert (n_cols > n_pad[XLO] + n_pad[XHI] && n_rows > n_pad[YLO] + n_pad[YHI] && n_pad[XLO] + n_pad[XHI] + n_pad[YLO] + n_pad[YHI] > 0 && cell_size > 0); for (row = 0; row < 4; row++) pad_t[row] = (size_t)n_pad[row]; if (periodic_cols) { /* A periodic grid wraps around */ for (row = pad_t[YHI]; (row + pad_t[YLO]) < n_rows_t; ++row) { nm = row * n_cols_t; /* Iterate over rows that contain data */ for (cell = 0; cell < pad_t[XLO]; ++cell) { /* Copy end of this row into first pad_t[XLO] columns: * X X 0 1 2 3 4 5 X X -> 4 5 0 1 2 3 4 5 X X */ memcpy (grid + (nm + cell) * cell_size, grid + (nm + n_cols_t + cell - pad_t[XLO] - pad_t[XHI]) * cell_size, cell_size); } for (cell = 0; cell < pad_t[XHI]; ++cell) { /* Copy start of this row into last n_pad[XHI] columns: * 4 5 0 1 2 3 4 5 X X -> 4 5 0 1 2 3 4 5 0 1 */ memcpy (grid + (nm + n_cols_t - cell - 1) * cell_size, grid + (nm + pad_t[XLO] + pad_t[XHI] - cell - 1) * cell_size, cell_size); } } } else { /* !periodic_cols */ for (row = pad_t[YHI]; (row + pad_t[YLO]) < n_rows_t; ++row) { nm = row * n_cols_t; /* Iterate over rows that contain data */ for (cell = 0; cell < pad_t[XLO]; ++cell) { /* Duplicate first pad_t[XLO] columns in this row: * 4 5 0 1 2 3 4 5 X X -> 0 0 0 1 2 3 4 5 X X */ memcpy (grid + (nm + cell) * cell_size, grid + (nm + pad_t[XLO]) * cell_size, cell_size); } for (cell = 0; cell < pad_t[XHI]; ++cell) { /* Duplicate last pad_t[XHI] columns in this row: * 0 0 0 1 2 3 4 5 X X -> 0 0 0 1 2 3 4 5 5 5 */ memcpy (grid + (nm + n_cols_t - cell - 1) * cell_size, grid + (nm + n_cols_t - pad_t[XHI] - 1) * cell_size, cell_size); } } } for (cell = 0; cell < pad_t[YHI]; ++cell) { nm = cell * n_cols_t; /* Duplicate pad_t[YHI] rows in the beginning */ memcpy (grid + nm * cell_size, grid + pad_t[YHI] * n_cols_t * cell_size, n_cols_t * cell_size); } nm2 = (n_rows_t - pad_t[YLO] - 1) * n_cols_t; for (cell = 0; cell < pad_t[YLO]; ++cell) { nm = (n_rows_t - cell - 1) * n_cols_t; /* Duplicate last pad_t[YLO] rows */ memcpy (grid + nm * cell_size, grid + nm2 * cell_size, n_cols_t * cell_size); } } /* Fill padding with zeros */ GMT_LOCAL void gmtnc_padding_zero (void *gridp, const unsigned n_cols, const unsigned n_rows, const unsigned *n_pad, size_t cell_size) { /* n_cols and n_rows are dimensions of the padded grid */ char *grid = (char*)gridp; size_t row, nm; size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows; size_t pad_t[4]; assert (n_cols > n_pad[XLO] + n_pad[XHI] && n_rows > n_pad[YLO] + n_pad[YHI] && n_pad[XLO] + n_pad[XHI] + n_pad[YLO] + n_pad[YHI] > 0 && cell_size > 0); for (row = 0; row < 4; row++) pad_t[row] = (size_t)n_pad[row]; /* Iterate over rows that contain data */ for (row = pad_t[YHI]; (row + pad_t[YLO]) < n_rows_t; ++row) { nm = row * n_cols_t; /* Zero n cells at beginning of row */ memset (grid + nm * cell_size, 0, pad_t[XLO] * cell_size); /* Zero n cells at end of row */ memset (grid + (nm + n_cols_t - pad_t[XHI]) * cell_size, 0, pad_t[XHI] * cell_size); } /* Zero pad_t[YHI] rows in the beginning */ memset (grid, 0, pad_t[YHI] * n_cols_t * cell_size); /* Zero last pad_t[YLO] rows */ memset(grid + (n_rows_t-pad_t[YLO]) * n_cols_t * cell_size, 0, pad_t[YLO] * n_cols_t * cell_size); } /* Fill mode for grid padding */ enum Grid_padding_mode { k_pad_fill_none = 0, /* Leave padded cells untouched */ k_pad_fill_zero, /* Fill padded grid cells with zeros */ k_pad_fill_copy, /* Padded cells get the value of their nearest neighbor */ k_pad_fill_copy_wrap /* Padded cells get wrapped values from the other side of the row (gridswith periodic columns) */ }; /* Add padding to a matrix/grid and reshape data */ GMT_LOCAL void gmtnc_pad_grid (void *gridp, const unsigned n_cols, const unsigned n_rows, const unsigned *n_pad, size_t cell_size, unsigned filltype) { /* n_cols and n_rows are dimensions of the grid without padding * cell_size is the size in bytes of each element in grid * n_pad[4] contains the number of cols/rows to pad on each side {W,E,S,N} * * Note: when grid is complex, we pass 2x n_rows */ char *grid = (char*)gridp; size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows; size_t new_row; size_t old_row = n_rows_t-1; size_t n_new_cols = n_cols_t + n_pad[XLO] + n_pad[XHI]; size_t n_new_rows = n_rows_t + n_pad[YLO] + n_pad[YHI]; #ifdef NC4_DEBUG fprintf (stderr, "pad grid w:%u e:%u s:%u n:%u\n", n_pad[XLO], n_pad[XHI], n_pad[YLO], n_pad[YHI]); #endif if (n_pad[XLO] + n_pad[XHI] + n_pad[YLO] + n_pad[YHI] == 0) return; /* nothing to pad */ assert (n_cols > 0 && n_rows > 0 && cell_size > 0); /* Reshape matrix */ if (n_pad[XLO] + n_pad[XHI] + n_pad[YHI] != 0) { /* When padding W, E, and N (not necessary when padding S only). */ for (new_row = n_new_rows - n_pad[YLO] - 1; new_row + 1 > n_pad[YHI]; --new_row, --old_row) { /* Copy original row to new row, bottom upwards */ void *from = grid + old_row * n_cols_t * cell_size; void *to = grid + (new_row * n_new_cols + ((size_t)n_pad[XLO])) * cell_size; if (n_pad[YHI] == 0) /* rows overlap! */ memmove (to, from, n_cols_t * cell_size); else /* no overlap, memcpy is safe */ memcpy (to, from, n_cols_t * cell_size); } } /* Fill padded grid cells */ switch (filltype) { case k_pad_fill_zero: gmtnc_padding_zero (grid, (const unsigned)n_new_cols, (const unsigned)n_new_rows, n_pad, cell_size); break; case k_pad_fill_copy: gmtnc_padding_copy (grid, (const unsigned)n_new_cols, (const unsigned)n_new_rows, n_pad, cell_size, false); break; case k_pad_fill_copy_wrap: gmtnc_padding_copy (grid, (const unsigned)n_new_cols, (const unsigned)n_new_rows, n_pad, cell_size, true); break; } } /* Remove padding from a matrix/grid and reshape data */ GMT_LOCAL void gmtnc_unpad_grid (void *gridp, const unsigned n_cols, const unsigned n_rows, const unsigned *n_pad, size_t cell_size) { /* n_cols and n_rows are dimensions of the grid without padding * cell_size is the size in bytes of each element in grid * n_pad[4] contains the number of cols/rows to pad on each side {W,E,S,N} * * Note: when grid is complex, we pass 2x n_rows */ char *grid = (char*)gridp; size_t n_old_cols = n_cols + n_pad[XLO] + n_pad[XHI]; size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows; size_t row; #ifdef NC4_DEBUG fprintf (stderr, "unpad grid w:%u e:%u s:%u n:%u\n", n_pad[XLO], n_pad[XHI], n_pad[YLO], n_pad[YHI]); #endif if (n_pad[XLO] + n_pad[XHI] + n_pad[YHI] == 0) return; /* nothing to unpad (we just ignore n_pad[YLO]) */ assert (n_cols > 0 && n_rows > 0 && cell_size > 0); /* Reshape matrix */ for (row = 0; row < n_rows_t; ++row) { size_t old_row = row + (size_t)n_pad[YHI]; void *from = grid + (old_row * n_old_cols + ((size_t)n_pad[XLO])) * cell_size; void *to = grid + row * n_cols_t * cell_size; /* Copy original row to new row */ if (n_pad[YHI] == 0) /* rows overlap! */ memmove (to, from, n_cols_t * cell_size); else /* no overlap, memcpy is safe */ memcpy (to, from, n_cols_t * cell_size); } } /* Ensure that repeating columns in geographic gridline registered grids * do not contain conflicting information */ GMT_LOCAL void gmtnc_grid_fix_repeat_col (struct GMT_CTRL *GMT, void *gridp, const unsigned n_cols, const unsigned n_rows, size_t cell_size) { /* Note: when grid is complex, pass 2x n_rows */ char *grid = (char*)gridp; unsigned n_conflicts = 0; size_t row; size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows; for (row = 0; row < n_rows_t; ++row) { char *first = grid + row * n_cols_t * cell_size; /* first element in row */ char *last = grid + (row * n_cols_t + n_cols_t - 1) * cell_size; /* last element in row */ if ( memcmp(last, first, cell_size) ) { /* elements differ: replace value of last element in row with value of first */ memcpy (last, first, cell_size); ++n_conflicts; } } if (n_conflicts) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected %u inconsistent values along east boundary of grid. Values fixed by duplicating west boundary.\n", n_conflicts); } /* Change the default chunk cache settings in the HDF5 library for all variables * in the nc-file. The settings apply for subsequent file opens/creates. */ static inline void gmtnc_setup_chunk_cache (void) { static bool already_setup = false; if (!already_setup) { nc_set_chunk_cache (NC_CACHE_SIZE, NC_CACHE_NELEMS, NC_CACHE_PREEMPTION); already_setup = true; } } GMT_LOCAL int gmtnc_grd_prep_io (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, double wesn[4], unsigned int *width, unsigned int *height, int *n_shift, unsigned origin[2], unsigned dim[2], unsigned origin2[2], unsigned dim2[2]) { /* Determines which rows and columns to extract from a grid, based on * w,e,s,n. This routine first rounds the w,e,s,n boundaries to the nearest * gridlines or pixels, then determines the first and last columns and rows, * and the width and height of the subset (in cells). */ bool is_gridline_reg, is_global, is_global_repeat; unsigned last_row, first_col, last_col, first_row; unsigned first_col2, last_col2; struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header); *n_shift = 0; memset (origin2, 0, 2 * sizeof(unsigned)); memset (dim2, 0, 2 * sizeof(unsigned)); is_global = gmt_grd_is_global (GMT, header); is_global_repeat = HH->grdtype == GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT; is_gridline_reg = header->registration != GMT_GRID_PIXEL_REG; #ifdef NC4_DEBUG GMT_Report (GMT->parent, GMT_MSG_WARNING, " x-region: %g %g, grid: %g %g\n", wesn[XLO], wesn[XHI], header->wesn[XLO], header->wesn[XHI]); GMT_Report (GMT->parent, GMT_MSG_WARNING, " y-region: %g %g, grid: %g %g\n", wesn[YLO], wesn[YHI], header->wesn[YLO], header->wesn[YHI]); #endif if (wesn[XLO] == 0 && wesn[XHI] == 0 && wesn[YLO] == 0 && wesn[YHI] == 0) { /* When -R... matches the exact dimensions of the grid: read whole grid */ *width = header->n_columns; *height = header->n_rows; first_col = first_row = 0; last_col = header->n_columns - 1; last_row = header->n_rows - 1; gmt_M_memcpy (wesn, header->wesn, 4, double); } else { /* Must deal with a subregion */ double x; x = fabs (header->wesn[YLO] - wesn[YLO]); /* if |x| < GMT_CONV4_LIMIT * header->inc[GMT_Y] we set wesn to the grid limit */ if (x > 0.0 && x < GMT_CONV4_LIMIT * header->inc[GMT_Y]) wesn[YLO] = header->wesn[YLO]; /* Avoid snafu */ x = fabs (header->wesn[YHI] - wesn[YHI]); /* if |x| < GMT_CONV4_LIMIT * header->inc[GMT_Y] we set wesn to the grid limit */ if (x > 0.0 && x < GMT_CONV4_LIMIT * header->inc[GMT_Y]) wesn[YHI] = header->wesn[YHI]; /* Avoid snafu */ if (wesn[YLO] < header->wesn[YLO] || wesn[YHI] > header->wesn[YHI]) return (GMT_GRDIO_DOMAIN_VIOLATION); /* Calling program goofed... */ /* Make sure w,e,s,n are proper multiples of x_inc,y_inc away from x_min,y_min */ gmt_M_err_pass (GMT, gmt_adjust_loose_wesn (GMT, wesn, header), HH->name); /* Global grids: ensure that wesn >= header->wesn (w+e only) */ if ( is_global ) { /* Deal with wrapping issues */ while (wesn[XLO] > header->wesn[XHI]) { wesn[XLO] -= 360; wesn[XHI] -= 360; } while (wesn[XLO] < header->wesn[XLO]) { wesn[XLO] += 360; wesn[XHI] += 360; } } else /* Apply a 10^-4 times inc sloppiness in the test [it was zero which was too harsh for sloppy grids] */ assert ((wesn[XLO]+GMT_CONV4_LIMIT*header->inc[GMT_X]) >= header->wesn[XLO] && (wesn[XHI]-GMT_CONV4_LIMIT*header->inc[GMT_X]) <= header->wesn[XHI]); /* Get dimension of subregion */ *width = urint ((wesn[XHI] - wesn[XLO]) * HH->r_inc[GMT_X]) + is_gridline_reg; *height = urint ((wesn[YHI] - wesn[YLO]) * HH->r_inc[GMT_Y]) + is_gridline_reg; /* Get first and last row and column numbers */ first_col = urint ((wesn[XLO] - header->wesn[XLO]) * HH->r_inc[GMT_X]); last_col = urint ((wesn[XHI] - header->wesn[XLO]) * HH->r_inc[GMT_X]) + is_gridline_reg - 1; first_row = urint ((header->wesn[YHI] - wesn[YHI]) * HH->r_inc[GMT_Y]); last_row = urint ((header->wesn[YHI] - wesn[YLO]) * HH->r_inc[GMT_Y]) + is_gridline_reg - 1; /* Adjust first_row */ if (HH->row_order == k_nc_start_south) first_row = header->n_rows - 1 - last_row; /* Global grids: if -R + padding is equal or exceeds grid bounds */ if (is_global && 1 + is_global_repeat + last_col - first_col >= header->n_columns) { /* Number of requested cols >= n_columns: read whole grid and shift */ *n_shift = -(int)first_col; first_col = 0; last_col = header->n_columns - 1; } else if (is_global && last_col + 1 > header->n_columns) { /* Subset of a global grid that wraps around east boundary. This means * we have to read the grid in two parts and stich them together. */ first_col2 = 0; last_col2 = is_global_repeat + last_col - header->n_columns; last_col = header->n_columns - 1; #ifdef NC4_DEBUG GMT_Report (GMT->parent, GMT_MSG_WARNING, "col2: %u %u\n", first_col2, last_col2); #endif origin2[0] = first_row; origin2[1] = first_col2; dim2[0] = *height; assert (first_col2 <= 1 + last_col2); /* check for unsigned overflow */ dim2[1] = 1 + last_col2 - first_col2; } assert (last_col + 1 <= header->n_columns); } /* Do not read last column from global repeating grids */ if (is_global_repeat && last_col + 1 == header->n_columns) --last_col; origin[0] = first_row; origin[1] = first_col; dim[0] = *height; assert (first_col <= 1 + last_col); /* check for unsigned overflow */ dim[1] = 1 + last_col - first_col; #ifdef NC4_DEBUG GMT_Report (GMT->parent, GMT_MSG_WARNING, "-> x-region: %g %g, grid: %g %g\n", wesn[XLO], wesn[XHI], header->wesn[XLO], header->wesn[XHI]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "-> y-region: %g %g, grid: %g %g\n", wesn[YLO], wesn[YHI], header->wesn[YLO], header->wesn[YHI]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "row: %u %u col: %u %u r_shift: %d\n", first_row, last_row, first_col, last_col, *n_shift); GMT_Report (GMT->parent, GMT_MSG_WARNING, "origin : %u,%u dim : %u,%u\n", origin[0], origin[1], dim[0], dim[1]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "origin2: %u,%u dim2: %u,%u\n", origin2[0], origin2[1], dim2[0], dim2[1]); #endif return GMT_NOERROR; } int gmtlib_is_nc_grid (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) { /* Returns GMT_NOERROR if NetCDF grid */ int ncid, z_id = -1, j = 0, nvars, ndims, err, old = false; nc_type z_type; char varname[GMT_GRID_VARNAME_LEN80]; struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header); /* Extract levels name from variable name */ strncpy (varname, HH->varname, GMT_GRID_VARNAME_LEN80-1); if (varname[0]) { j = 0; while (varname[j] && varname[j] != '[' && varname[j] != '(') j++; if (varname[j]) varname[j] = '\0'; } if (!strcmp (HH->name, "=")) return (GMT_GRDIO_NC_NO_PIPE); /* Open the file and look for the required variable */ if (gmt_access (GMT, HH->name, F_OK)) return (GMT_GRDIO_FILE_NOT_FOUND); if ((err = nc_open (HH->name, NC_NOWRITE, &ncid))) return (GMT_GRDIO_OPEN_FAILED); if (!nc_inq_dimid (ncid, "xysize", &z_id)) { /* Old style GMT netCDF grid */ old = true; if (nc_inq_varid (ncid, "z", &z_id)) return (GMT_GRDIO_NO_VAR); } else if (varname[0]) { /* ? used */ if (nc_inq_varid (ncid, varname, &z_id)) return (GMT_GRDIO_NO_VAR); } else { /* Look for first 2D grid, with fallback to first higher-dimension (3-4D) grid if 2D not found */ int ID = -1, dim = 0; nc_inq_nvars (ncid, &nvars); for (j = 0; j < nvars && z_id < 0; j++) { gmt_M_err_trap (nc_inq_varndims (ncid, j, &ndims)); if (ndims == 2) /* Found the first 2-D grid */ z_id = j; else if (ID == -1 && ndims > 2 && ndims < 5) { /* Also look for higher-dim grid in case no 2-D */ ID = j; dim = ndims; } } if (z_id < 0) { /* No 2-D grid found, check if we found a higher dimension cube */ if (ID == -1) return (GMT_GRDIO_NO_2DVAR); /* No we didn't */ z_id = ID; /* Pick the higher dimensioned cube instead, get its name, and warn */ nc_inq_varname (ncid, z_id, varname); GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No 2-D array in file %s. Selecting first 2-D slice in the %d-D array %s\n", HH->name, dim, varname); } } gmt_M_err_trap (nc_inq_vartype (ncid, z_id, &z_type)); switch (z_type) { case NC_UBYTE: header->type = old ? GMT_GRID_IS_CB : GMT_GRID_IS_NB; HH->orig_datatype = GMT_UCHAR; break; case NC_BYTE: header->type = old ? GMT_GRID_IS_CB : GMT_GRID_IS_NB; HH->orig_datatype = GMT_CHAR; break; case NC_SHORT: header->type = old ? GMT_GRID_IS_CS : GMT_GRID_IS_NS; HH->orig_datatype = GMT_SHORT; break; case NC_INT: header->type = old ? GMT_GRID_IS_CI : GMT_GRID_IS_NI; HH->orig_datatype = GMT_INT; break; case NC_FLOAT: header->type = old ? GMT_GRID_IS_CF : GMT_GRID_IS_NF; HH->orig_datatype = GMT_FLOAT; break; case NC_DOUBLE: header->type = old ? GMT_GRID_IS_CD : GMT_GRID_IS_ND; HH->orig_datatype = GMT_DOUBLE; break; default: header->type = k_grd_unknown_fmt; break; } nc_close (ncid); return GMT_NOERROR; } int gmt_nc_read_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) { return (gmtnc_grd_info (GMT, header, 'r')); } int gmt_nc_update_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) { return (gmtnc_grd_info (GMT, header, 'u')); } int gmt_nc_write_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) { return (gmtnc_grd_info (GMT, header, 'w')); } int gmt_nc_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) { /* header: grid structure header * grid: array with final grid * wesn: Sub-region to extract [Use entire file if 0,0,0,0] * padding: # of empty rows/columns to add on w, e, s, n of grid, respectively * complex_mode: &4 | &8 if complex array is to hold real (4) and imaginary (8) parts (otherwise read as real only) * Note: The file has only real values, we simply allow space in the complex array * for real and imaginary parts when processed by grdfft etc. * * Reads a subset of a grdfile and optionally pads the array with extra rows and columns * header values for n_columns and n_rows are reset to reflect the dimensions of the logical array, * not the physical size (i.e., the padding is not counted in n_columns and n_rows) */ bool adj_nan_value; /* if we need to change the fill value */ int err; /* netcdf errors */ int n_shift; int error; unsigned dim[2], dim2[2], origin[2], origin2[2]; /* dimension and origin {y,x} of subset to read from netcdf */ unsigned width = 0, height = 0; size_t width_t, height_t, row, stride_t; uint64_t imag_offset; gmt_grdfloat *pgrid = NULL; struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header); /* Check type: is file in old NetCDF format or not at all? */ if (GMT->session.grdformat[header->type][0] == 'c') return (gmt_cdf_read_grd (GMT, header, grid, wesn, pad, complex_mode)); else if (GMT->session.grdformat[header->type][0] != 'n') return (NC_ENOTNC); if ((error = gmt_M_err_fail (GMT, gmtnc_grd_prep_io (GMT, header, wesn, &width, &height, &n_shift, origin, dim, origin2, dim2), HH->name))) return (error); /* Set stride and offset if complex */ (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */ pgrid = grid + imag_offset; /* Start of this complex component (or start of non-complex grid) */ #ifdef NC4_DEBUG GMT_Report (GMT->parent, GMT_MSG_WARNING, " wesn: %g %g %g %g\n", wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->wesn: %g %g %g %g\n", header->wesn[XLO], header->wesn[XHI], header->wesn[YLO], header->wesn[YHI]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->registration:%u\n", header->registration); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->row_order: %s\n", HH->row_order == k_nc_start_south ? "S->N" : "N->S"); GMT_Report (GMT->parent, GMT_MSG_WARNING, "width: %3d height:%3d\n", width, height); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->n_columns: %3d head->n_rows:%3d\n", header->n_columns, header->n_rows); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->mx: %3d head->my:%3d\n", header->mx, header->my); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->nm: %3d head->size:%3d\n", header->nm, header->size); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->t-index %d,%d,%d\n", HH->t_index[0], HH->t_index[1], HH->t_index[2]); GMT_Report (GMT->parent, GMT_MSG_WARNING, " pad xlo:%u xhi:%u ylo:%u yhi:%u\n", pad[XLO], pad[XHI], pad[YLO], pad[YHI]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->pad xlo:%u xhi:%u ylo:%u yhi:%u\n", header->pad[XLO], header->pad[XHI], header->pad[YLO], header->pad[YHI]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->BC xlo:%u xhi:%u ylo:%u yhi:%u\n", HH->BC[XLO], HH->BC[XHI], HH->BC[YLO], HH->BC[YHI]); GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->grdtype:%u %u\n", HH->grdtype, GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT); GMT_Report (GMT->parent, GMT_MSG_WARNING, "imag_offset: %" PRIu64 "\n", imag_offset); #endif /* Open the NetCDF file */ if (!strcmp (HH->name,"=")) return (GMT_GRDIO_NC_NO_PIPE); gmtnc_setup_chunk_cache(); gmt_M_err_trap (nc_open (HH->name, NC_NOWRITE, &HH->ncid)); /* Read grid */ if (dim2[1] == 0) gmtnc_io_nc_grid (GMT, header, dim, origin, HH->stride, k_get_netcdf, pgrid + HH->data_offset); else { /* Read grid in two parts */ unsigned int stride_or_width = HH->stride != 0 ? HH->stride : width; gmtnc_io_nc_grid (GMT, header, dim, origin, stride_or_width, k_get_netcdf, pgrid + HH->data_offset); gmtnc_io_nc_grid (GMT, header, dim2, origin2, stride_or_width, k_get_netcdf, pgrid + HH->data_offset + dim[1]); } /* If we need to shift grid */ if (n_shift) { #ifdef NC4_DEBUG GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtnc_right_shift_grid: %d\n", n_shift); #endif gmtnc_right_shift_grid (pgrid, dim[1], dim[0], n_shift, sizeof(grid[0])); } /* If dim[1] + dim2[1] was < requested width: wrap-pad east border */ if (gmt_grd_is_global(GMT, header) && width > dim[1] + dim2[1]) { unsigned int fix_pad[4] = {0,0,0,0}; fix_pad[XHI] = width - dim[1] - dim2[1]; gmtnc_pad_grid (pgrid, width - fix_pad[XHI], height, fix_pad, sizeof(grid[0]), k_pad_fill_copy_wrap); } else assert (width == dim[1] + dim2[1]); /* Get stats */ width_t = (size_t)width; height_t = (size_t)height; stride_t = (size_t)HH->stride; header->z_min = DBL_MAX; header->z_max = -DBL_MAX; adj_nan_value = !isnan (header->nan_value); HH->has_NaNs = GMT_GRID_NO_NANS; /* We are about to check for NaNs and if none are found we retain 1, else 2 */ for (row = 0; row < height_t; ++row) { gmt_grdfloat *p_data = pgrid + row * (HH->stride ? stride_t : width_t); unsigned col; for (col = 0; col < width; col ++) { if (adj_nan_value && p_data[col] == header->nan_value) { p_data[col] = (gmt_grdfloat)NAN; HH->has_NaNs = GMT_GRID_HAS_NANS; continue; } else if (!isnan (p_data[col])) { header->z_min = MIN (header->z_min, p_data[col]); header->z_max = MAX (header->z_max, p_data[col]); } else HH->has_NaNs = GMT_GRID_HAS_NANS; } } /* Check limits */ if (header->z_min > header->z_max) { header->z_min = NAN; header->z_max = NAN; GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No valid values in grid [%s]\n", HH->name); } else { /* Report z-range of grid (with scale and offset applied): */ unsigned int level; #ifdef NC4_DEBUG level = GMT_MSG_WARNING; #else level = GMT_MSG_DEBUG; #endif GMT_Report (GMT->parent, level, "packed z-range: [%g,%g]\n", header->z_min, header->z_max); } /* Flip grid upside down */ if (HH->row_order == k_nc_start_south) gmt_grd_flip_vertical (pgrid + HH->data_offset, width, height, HH->stride, sizeof(grid[0])); /* Add padding with border replication */ gmtnc_pad_grid (pgrid, width, height, pad, sizeof(grid[0]), k_pad_fill_zero); #ifdef NC4_DEBUG if (header->size < 160) { unsigned n; unsigned pad_x = pad[XLO] + pad[XHI]; unsigned stride = HH->stride ? HH->stride : width; gmt_grdfloat *p_data = pgrid + HH->data_offset; for (n = 0; n < (stride + pad_x) * (height + pad[YLO] + pad[YHI]); n++) { if (n % (stride + pad_x) == 0) fprintf (stderr, "\n"); fprintf (stderr, "%4.0f", p_data[n]); } fprintf (stderr, "\n"); } #endif /* Adjust header */ gmt_M_memcpy (header->wesn, wesn, 4, double); header->n_columns = width; header->n_rows = height; gmt_M_err_trap (nc_close (HH->ncid)); return GMT_NOERROR; } int gmt_nc_write_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) { /* header: grid structure header * grid: array with final grid * wesn: Sub-region to write out [Use entire file if 0,0,0,0] * padding: # of empty rows/columns to add on w, e, s, n of grid, respectively * complex_mode: &4 | &8 if complex array is to hold real (4) and imaginary (8) parts (otherwise read as real only) * Note: The file has only real values, we simply allow space in the complex array * for real and imaginary parts when processed by grdfft etc. */ int status = NC_NOERR; bool adj_nan_value; /* if we need to change the fill value */ bool do_round = true; /* if we need to round to integral */ unsigned int width, height, *actual_col = NULL; unsigned int dim[2], origin[2]; /* dimension and origin {y,x} of subset to write to netcdf */ int first_col, last_col, first_row, last_row; uint64_t imag_offset; size_t n, nm; size_t width_t, height_t; double limit[2]; /* minmax of z variable */ gmt_grdfloat *pgrid = NULL; struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header); /* Determine the value to be assigned to missing data, if not already done so */ switch (header->type) { case GMT_GRID_IS_NB: if (isnan (header->nan_value)) header->nan_value = NC_MIN_BYTE; break; case GMT_GRID_IS_NS: if (isnan (header->nan_value)) header->nan_value = NC_MIN_SHORT; break; case GMT_GRID_IS_NI: if (isnan (header->nan_value)) header->nan_value = NC_MIN_INT; break; case GMT_GRID_IS_ND: #ifndef DOUBLE_PRECISION_GRID GMT_Report (GMT->parent, GMT_MSG_WARNING, "Precision loss! GMT's internal grid representation is 32-bit float.\n"); #endif /* Intentionally fall through */ default: /* don't round float */ do_round = false; } gmt_M_err_pass (GMT, gmt_grd_prep_io (GMT, header, wesn, &width, &height, &first_col, &last_col, &first_row, &last_row, &actual_col), HH->name); gmt_M_free (GMT, actual_col); /* Adjust header */ gmt_M_memcpy (header->wesn, wesn, 4, double); header->n_columns = width; header->n_rows = height; /* Adjust first_row */ if (HH->row_order == k_nc_start_south) first_row = header->n_rows - 1 - last_row; /* Write grid header without closing file afterwards */ gmtnc_setup_chunk_cache(); status = gmtnc_grd_info (GMT, header, 'W'); if (status != NC_NOERR) goto nc_err; /* Set stride and offset if complex */ (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */ pgrid = grid + imag_offset; /* Beginning of this complex component (or the regular non-complex grid) */ /* Remove padding from grid */ gmtnc_unpad_grid (pgrid, width, height, pad, sizeof(grid[0])); /* Check that repeating columns do not contain conflicting information */ if (HH->grdtype == GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT) gmtnc_grid_fix_repeat_col (GMT, pgrid, width, height, sizeof(grid[0])); /* Flip grid upside down */ if (HH->row_order == k_nc_start_south) gmt_grd_flip_vertical (pgrid, width, height, 0, sizeof(grid[0])); /* Get stats */ header->z_min = DBL_MAX; header->z_max = -DBL_MAX; adj_nan_value = !isnan (header->nan_value); n = 0; width_t = (size_t)width; height_t = (size_t)height; nm = width_t * height_t; while (n < nm) { if (adj_nan_value && isnan (pgrid[n])) pgrid[n] = header->nan_value; else if (!isnan (pgrid[n])) { if (do_round) pgrid[n] = rintf (pgrid[n]); /* round to int */ header->z_min = MIN (header->z_min, pgrid[n]); header->z_max = MAX (header->z_max, pgrid[n]); } n++; } /* Write grid */ dim[0] = height, dim[1] = width; origin[0] = first_row, origin[1] = first_col; status = gmtnc_io_nc_grid (GMT, header, dim, origin, 0, k_put_netcdf, pgrid); if (status != NC_NOERR) goto nc_err; if (header->z_min <= header->z_max) { /* Warn if z-range exceeds the precision of a single precision float: */ static const uint32_t exp2_24 = 0x1000000; /* exp2 (24) */ unsigned int level; if (fabs(header->z_min) >= exp2_24 || fabs(header->z_max) >= exp2_24) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "The z-range, [%g,%g], might exceed the significand's precision of 24 bits; round-off errors may occur.\n", header->z_min, header->z_max); /* Report z-range of grid (with scale and offset applied): */ #ifdef NC4_DEBUG level = GMT_MSG_WARNING; #else level = GMT_MSG_DEBUG; #endif GMT_Report (GMT->parent, level, "packed z-range: [%g,%g]\n", header->z_min, header->z_max); /* Limits need to be written in actual, not internal grid, units: */ limit[0] = header->z_min * header->z_scale_factor + header->z_add_offset; limit[1] = header->z_max * header->z_scale_factor + header->z_add_offset; } else { GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No valid values in grid [%s]\n", HH->name); limit[0] = limit[1] = NAN; /* Set limit to NaN */ } status = nc_put_att_double (HH->ncid, HH->z_id, "actual_range", NC_DOUBLE, 2, limit); if (status != NC_NOERR) goto nc_err; /* Close grid */ status = nc_close (HH->ncid); if (status != NC_NOERR) goto nc_err; return GMT_NOERROR; nc_err: /* exit gracefully */ nc_close(HH->ncid); /* close nc-file */ unlink (HH->name); /* remove nc-file */ if (status == NC_ERANGE) { /* report out of range z variable */ GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cannot write format %s.\n", GMT->session.grdformat[header->type]); GMT_Report (GMT->parent, GMT_MSG_ERROR, "The packed z-range, [%g,%g], exceeds the maximum representable size. Adjust scale and offset parameters or remove out-of-range values.\n", header->z_min, header->z_max); } return status; } /* Examine the netCDF data cube and determine if it is a 3-D cube and return the knots */ int gmt_examine_nc_cube (struct GMT_CTRL *GMT, char *file, uint64_t *nz, double **zarray) { int i, err, ID = -1, dim = 0, ncid, z_id = -1, ids[5] = {-1,-1,-1,-1,-1}, dims[5], nvars; int has_vector,ndims = 0, z_dim, status; uint64_t n_layers = 0; size_t lens[5]; char varname[GMT_GRID_VARNAME_LEN80], dimname[GMT_GRID_UNIT_LEN80], z_units[GMT_GRID_UNIT_LEN80]; double *z = NULL; gmt_M_err_trap (nc_open (file, NC_NOWRITE, &ncid)); gmt_M_err_trap (nc_inq_nvars (ncid, &nvars)); i = 0; while (i < nvars && z_id < 0) { /* Look for first 3D grid, with fallback to first higher-dimension (3-4D) grid if 3D not found */ gmt_M_err_trap (nc_inq_varndims (ncid, i, &ndims)); if (ndims == 3) /* Found the first 3-D grid */ z_id = i; else if (ID == -1 && ndims > 3 && ndims < 5) { /* Also look for higher-dim grid in case no 3-D */ ID = i; dim = ndims; } i++; } if (z_id < 0) { /* No 3-D grid found, check if we found a higher dimension cube */ if (ID == -1) { /* No we didn't */ GMT_Report (GMT->parent, GMT_MSG_ERROR, "No 3-D data cube found in file %s.\n", file); return GMT_GRDIO_NO_2DVAR; } z_id = ID; /* Pick the higher dimensioned cube instead, get its name, and warn */ nc_inq_varname (ncid, z_id, varname); GMT_Report (GMT->parent, GMT_MSG_WARNING, "No 3-D array in file %s. Selecting first 3-D slice in the %d-D array %s\n", file, dim, varname); } gmt_M_err_trap (nc_inq_vardimid (ncid, z_id, dims)); /* Get the ids of the x and y (and depth and time) coordinate variables */ for (i = 0; i < ndims; i++) { gmt_M_err_trap (nc_inq_dim (ncid, dims[i], dimname, &lens[i])); if ((status = nc_inq_varid (ncid, dimname, &ids[i])) != NC_NOERR) GMT_Report (GMT->parent, GMT_MSG_WARNING, "\"%s\", %s\n\tIf something bad happens later, try importing via GDAL.\n", dimname, nc_strerror(status)); } z_dim = ndims-3; n_layers = lens[z_dim]; /* Create enough memory to store the level-coordinate values */ z = gmt_M_memory (GMT, NULL, n_layers, double); /* Get information about z variable */ gmtnc_get_units (GMT, ncid, ids[z_dim], z_units); /* Look for the z-coordinate vector */ if ((has_vector = nc_get_var_double (ncid, ids[z_dim], z))) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "No 3rd-dimension coordinate vector found in %s\n", file); return GMT_GRDIO_NO_2DVAR; } *zarray = z; *nz = n_layers; return GMT_NOERROR; } /* Write a 3-D cube to file; cube is represented internally by a stack of 2-D grids */ #define GMT_WRITE_CUBE_LAYERS 0 int gmt_write_nc_cube (struct GMT_CTRL *GMT, struct GMT_GRID **G, uint64_t nlayers, double *layer, char *file, unsigned int mode) { /* Depending on mode, we either write individual layer grid files or a single 3-D data cube */ uint64_t k; if (mode == GMT_WRITE_CUBE_LAYERS) { char gfile[PATH_MAX] = {""}; for (k = 0; k < nlayers; k++) { sprintf (gfile, file, layer[k]); if (GMT_Write_Data (GMT->parent, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, gfile, G[k]) != GMT_NOERROR) { return (GMT->parent->error); } } return (GMT_NOERROR); } else { /* Here we must write a 3-D netcdf data cube */ int status = NC_NOERR; bool adj_nan_value; /* if we need to change the fill value */ bool do_round = true; /* if we need to round to integral */ unsigned int width, height; unsigned int dim[2], origin[2]; /* dimension and origin {y,x} of subset to write to netcdf */ int first_col, first_row, last_row; // int last_col; size_t n, nm; size_t width_t, height_t; double level_min, level_max, limit[2]; /* minmax of level variable */ gmt_grdfloat *pgrid = NULL; struct GMT_GRID_HEADER *header = G[0]->header; struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header); width = header->n_columns; height = header->n_rows; /* Determine the value to be assigned to missing data, if not already done so */ switch (header->type) { case GMT_GRID_IS_NB: if (isnan (header->nan_value)) header->nan_value = NC_MIN_BYTE; break; case GMT_GRID_IS_NS: if (isnan (header->nan_value)) header->nan_value = NC_MIN_SHORT; break; case GMT_GRID_IS_NI: if (isnan (header->nan_value)) header->nan_value = NC_MIN_INT; break; case GMT_GRID_IS_ND: #ifndef DOUBLE_PRECISION_GRID GMT_Report (GMT->parent, GMT_MSG_WARNING, "Precision loss! GMT's internal grid representation is 32-bit float.\n"); #endif /* Intentionally fall through */ default: /* don't round float */ do_round = false; } first_col = first_row = 0; // last_col = header->n_columns - 1; last_row = header->n_rows - 1; level_min = DBL_MAX; level_max = -DBL_MAX; /* Adjust first_row */ if (HH->row_order == k_nc_start_south) first_row = header->n_rows - 1 - last_row; /* Write grid header without closing file afterwards */ gmtnc_setup_chunk_cache(); status = gmtnc_grd_info (GMT, header, 'W'); if (status != NC_NOERR) goto nc_err; /* Get stats */ adj_nan_value = !isnan (header->nan_value); dim[0] = height, dim[1] = width; for (k = 0; k < nlayers; k++) { /* Write each layer separately */ pgrid = G[k]->data; /* Remove padding from grid */ gmtnc_unpad_grid (pgrid, width, height, header->pad, sizeof(pgrid[0])); /* Check that repeating columns do not contain conflicting information */ if (HH->grdtype == GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT) gmtnc_grid_fix_repeat_col (GMT, pgrid, width, height, sizeof(pgrid[0])); /* Flip grid upside down */ if (HH->row_order == k_nc_start_south) gmt_grd_flip_vertical (pgrid, width, height, 0, sizeof(pgrid[0])); /* Update stats */ n = 0; width_t = (size_t)width; height_t = (size_t)height; nm = width_t * height_t; while (n < nm) { if (adj_nan_value && isnan (pgrid[n])) pgrid[n] = header->nan_value; else if (!isnan (pgrid[n])) { if (do_round) pgrid[n] = rintf (pgrid[n]); /* round to int */ level_min = MIN (level_min, pgrid[n]); level_max = MAX (level_max, pgrid[n]); } n++; } /* Write grid layer */ origin[0] = first_row, origin[1] = first_col; status = gmtnc_io_nc_grid (GMT, header, dim, origin, 0, k_put_netcdf, pgrid); if (status != NC_NOERR) goto nc_err; } if (level_min <= level_max) { /* Warn if level-range exceeds the precision of a single precision float: */ static const uint32_t exp2_24 = 0x1000000; /* exp2 (24) */ unsigned int level; if (fabs(level_min) >= exp2_24 || fabs(level_max) >= exp2_24) GMT_Report (GMT->parent, GMT_MSG_WARNING, "The level-range, [%g,%g], might exceed the significand's precision of 24 bits; round-off errors may occur.\n", level_min, level_max); /* Report level-range of grid layer (with scale and offset applied): */ #ifdef NC4_DEBUG level = GMT_MSG_WARNING; #else level = GMT_MSG_DEBUG; #endif GMT_Report (GMT->parent, level, "packed z-range: [%g,%g]\n", level_min, level_max); /* Limits need to be written in actual, not internal grid, units: */ limit[0] = level_min * header->z_scale_factor + header->z_add_offset; limit[1] = level_max * header->z_scale_factor + header->z_add_offset; } else { GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No valid values in grid [%s]\n", HH->name); limit[0] = limit[1] = NAN; /* Set limit to NaN */ } status = nc_put_att_double (HH->ncid, HH->z_id, "actual_range", NC_DOUBLE, 2, limit); if (status != NC_NOERR) goto nc_err; /* Close grid */ status = nc_close (HH->ncid); if (status != NC_NOERR) goto nc_err; return GMT_NOERROR; nc_err: /* exit gracefully */ nc_close(HH->ncid); /* close nc-file */ unlink (HH->name); /* remove nc-file */ if (status == NC_ERANGE) { /* report out of range z variable */ GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cannot write format %s.\n", GMT->session.grdformat[header->type]); GMT_Report (GMT->parent, GMT_MSG_ERROR, "The packed z-range, [%g,%g], exceeds the maximum representable size. Adjust scale and offset parameters or remove out-of-range values.\n", level_min, level_max); } return status; } }