diff --git a/doc/rst/source/grdcut.rst b/doc/rst/source/grdcut.rst index 28e172c77b5..e3e8194c70d 100644 --- a/doc/rst/source/grdcut.rst +++ b/doc/rst/source/grdcut.rst @@ -32,16 +32,18 @@ Description ----------- **grdcut** will produce a new *outgrid* file which is a subregion of -*ingrid*. The subregion is specified with **-R** as in other programs; +*ingrid*. The subregion may be specified with **-R** as in other programs; the specified range must not exceed the range of *ingrid* (but see **-N**). If in doubt, run :doc:`grdinfo` to check range. Alternatively, define the subregion indirectly via a range check on the node values or via distances from a -given point. Finally, you can use **-J** for oblique projections to determine -the corresponding rectangular **-R** setting that will give a grid that fully -covers the oblique domain. **Note**: If input grid is actually an image, then -only options **-R** and **-G** are supported, i.e., you can cut out a sub-region only. -Complementary to **grdcut** there is :doc:`grdpaste`, which -will join together two grid files along a common edge. +fixed point. Finally, you can use **-J** for oblique projections to determine +the corresponding rectangular **-R** setting that will give a subregion that fully +covers the oblique domain. **Note**: If the input grid is actually an image (gray-scale, +RGB, or RGBA), then options **-N** and **-Z** are unavailable, while for multi-layer +Geotiff files only options **-R**, **-S** and **-G** are supported, i.e., you can cut out +a sub-region only (which we do via gdal_translate if you have multiple bands). +Complementary to **grdcut** there is :doc:`grdpaste`, which will join together +two grid files (not images) along a common edge. Required Arguments ------------------ @@ -179,7 +181,21 @@ To determine what grid region and resolution (in text format) most suitable for that is using an oblique projection to display the remote Earth Relief data grid, try:: gmt grdcut @earth_relief -R270/20/305/25+r -JOc280/25.5/22/69/24c -D+t -V - + +Notes +----- + +If the input file is a geotiff with multiple data bands then the output format will +depend on your selection (if any) of the bands to keep: If you do not specify +any bands (which means we consider all the available bands) or you select more +than one band, then the output file can either be another geotiff (if you give +a .tif[f] extension) or it can be a multiband netCDF file (if you give a .nc or .grd +extension). If you select a single band from the input geotiff then GMT will +normally read that in as a single grid layer and thus write a netCDF grid (unless +you append another grid format specifier). However, if your output filename has +a .tif[f] extension then we will instead write it as a one-band geotiff. +All geotiff output operations are done via GDAL's gdal_translate. + See Also -------- diff --git a/src/gmt_api.c b/src/gmt_api.c index 2e74a8a0065..6d00888560d 100644 --- a/src/gmt_api.c +++ b/src/gmt_api.c @@ -4502,8 +4502,12 @@ GMT_LOCAL struct GMT_IMAGE *gmtapi_import_image (struct GMTAPI_CTRL *API, int ob I_obj->data = gmt_M_memory (GMT, NULL, size * I_obj->header->n_bands, unsigned char); else if (I_obj->type <= GMT_USHORT) I_obj->data = gmt_M_memory (GMT, NULL, size * I_obj->header->n_bands, unsigned short); - else if (I_obj->type <= GMT_UINT) - I_obj->data = gmt_M_memory (GMT, NULL, size * I_obj->header->n_bands, unsigned int); + else if (I_obj->type <= GMT_UINT) + I_obj->data = gmt_M_memory (GMT, NULL, size * I_obj->header->n_bands, unsigned int); + else if (I_obj->type <= GMT_FLOAT) + I_obj->data = gmt_M_memory (GMT, NULL, size * I_obj->header->n_bands, float); + else if (I_obj->type <= GMT_DOUBLE) + I_obj->data = gmt_M_memory (GMT, NULL, size * I_obj->header->n_bands, float); /* Not doing double yet */ else { GMT_Report (API, GMT_MSG_ERROR, "Unsupported image data type %d\n", I_obj->type); return_null (API, GMT_NOT_A_VALID_TYPE); diff --git a/src/gmt_gdalread.c b/src/gmt_gdalread.c index 09ba268f807..8f670961661 100644 --- a/src/gmt_gdalread.c +++ b/src/gmt_gdalread.c @@ -758,7 +758,7 @@ int gmt_gdalread (struct GMT_CTRL *GMT, char *gdal_filename, struct GMT_GDALREAD int *whichBands = NULL; int64_t *rowVec = NULL, *colVec = NULL; int64_t off, i_x_nXYSize, startColPos = 0, indent = 0, col_indent = 0, nXSize_withPad = 0, nYSize_withPad; - int64_t n_alloc, n, m, nn, mm, ij; + int64_t n_alloc, n, m, nn, mm, ij, layer_size, layer_offset = 0; unsigned char *tmp = NULL; double adfMinMax[2]; double adfGeoTransform[6]; @@ -1041,7 +1041,8 @@ int gmt_gdalread (struct GMT_CTRL *GMT, char *gdal_filename, struct GMT_GDALREAD if (nReqBands) nBands = MIN(nBands,nReqBands); /* If a band selection was made */ /* Determine allocation size for the subset (actual allocation depends on whether a pointer is passed in or not) */ - n_alloc = ((size_t)nBands) * ((size_t)(nBufXSize[0] + nBufXSize[1] + MAX(pad_w[0],pad_w[1])) + MAX(pad_e[0],pad_e[1])) * ((size_t)nBufYSize + pad_s + pad_n); + layer_size = ((size_t)(nBufXSize[0] + nBufXSize[1] + MAX(pad_w[0],pad_w[1])) + MAX(pad_e[0],pad_e[1])) * ((size_t)nBufYSize + pad_s + pad_n); + n_alloc = ((size_t)nBands) * layer_size; switch (GDALGetRasterDataType(hBand)) { case GDT_Byte: @@ -1145,7 +1146,6 @@ int gmt_gdalread (struct GMT_CTRL *GMT, char *gdal_filename, struct GMT_GDALREAD rowVec = gmt_M_memory(GMT, NULL, (nRowsPerBlock * nRGBA) * nBlocks, size_t); for (m = 0; m < nY; m++) rowVec[m] = m * nX[piece]; /* Steps in pixels as we go down rows */ colVec = gmt_M_memory(GMT, NULL, (nX[piece]+pad_w[piece]+pad_e[piece]) * nRGBA, size_t); /* For now this will be used only to select BIP ordering */ - for (i = 0; i < nBands; i++) { if (!nReqBands) /* No band selection, read them sequentially */ hBand = GDALGetRasterBand(hDataset, i+1); @@ -1289,7 +1289,7 @@ int gmt_gdalread (struct GMT_CTRL *GMT, char *gdal_filename, struct GMT_GDALREAD /* PW: Special case of jp2 tile with possible NaNs that are not recognized without the .aux.xml file */ float f_NaN = GMT->session.f_NaN; /* Shorthand */ for (m = startRow, mm = 0; m < endRow; m++, mm++) { - nn = (pad_w[piece]+m)*nXSize_withPad + startColPos; + nn = layer_offset + (pad_w[piece]+m)*nXSize_withPad + startColPos; for (n = fliplr ? nXSize[piece]-1 : 0; fliplr ? n >= 0 : n < nXSize[piece]; fliplr ? n-- : n++) if (prhs->f_ptr.active) { int16_t tmpI16; @@ -1302,7 +1302,7 @@ int gmt_gdalread (struct GMT_CTRL *GMT, char *gdal_filename, struct GMT_GDALREAD } else { /* All other short int cases */ for (m = startRow, mm = 0; m < endRow; m++, mm++) { - nn = (pad_w[piece]+m)*nXSize_withPad + startColPos; + nn = layer_offset + (pad_w[piece]+m)*nXSize_withPad + startColPos; for (n = fliplr ? nXSize[piece]-1 : 0; fliplr ? n >= 0 : n < nXSize[piece]; fliplr ? n-- : n++) if (prhs->f_ptr.active) { int16_t tmpI16; @@ -1317,7 +1317,7 @@ int gmt_gdalread (struct GMT_CTRL *GMT, char *gdal_filename, struct GMT_GDALREAD case GDT_UInt16: if (!prhs->f_ptr.active) Ctrl->UInt16.active = true; for (m = startRow, mm = 0; m < endRow; m++, mm++) { - nn = (pad_w[piece]+m)*nXSize_withPad + startColPos; + nn = layer_offset + (pad_w[piece]+m)*nXSize_withPad + startColPos; for (n = fliplr ? nXSize[piece]-1 : 0; fliplr ? n >= 0 : n < nXSize[piece]; fliplr ? n-- : n++) if (prhs->f_ptr.active) { uint16_t tmpUI16; @@ -1331,7 +1331,7 @@ int gmt_gdalread (struct GMT_CTRL *GMT, char *gdal_filename, struct GMT_GDALREAD case GDT_Int32: if (!prhs->f_ptr.active) Ctrl->Int32.active = true; for (m = startRow, mm = 0; m < endRow; m++, mm++) { - nn = (pad_w[piece]+m)*nXSize_withPad + startColPos; + nn = layer_offset + (pad_w[piece]+m)*nXSize_withPad + startColPos; for (n = fliplr ? nXSize[piece]-1 : 0; fliplr ? n >= 0 : n < nXSize[piece]; fliplr ? n-- : n++) if (prhs->f_ptr.active) { int32_t tmpI32; @@ -1345,7 +1345,7 @@ int gmt_gdalread (struct GMT_CTRL *GMT, char *gdal_filename, struct GMT_GDALREAD case GDT_UInt32: if (!prhs->f_ptr.active) Ctrl->UInt32.active = true; for (m = startRow, mm = 0; m < endRow; m++, mm++) { - nn = (pad_w[piece]+m)*nXSize_withPad + startColPos; + nn = layer_offset + (pad_w[piece]+m)*nXSize_withPad + startColPos; for (n = fliplr ? nXSize[piece]-1 : 0; fliplr ? n >= 0 : n < nXSize[piece]; fliplr ? n-- : n++) if (prhs->f_ptr.active) { int32_t tmpUI32; @@ -1359,7 +1359,7 @@ int gmt_gdalread (struct GMT_CTRL *GMT, char *gdal_filename, struct GMT_GDALREAD case GDT_Float32: Ctrl->Float.active = true; for (m = startRow, mm = 0; m < endRow; m++, mm++) { - nn = (pad_w[piece]+m)*nXSize_withPad + startColPos; + nn = layer_offset + (pad_w[piece]+m)*nXSize_withPad + startColPos; for (n = 0; n < nXSize[piece]; n++) { memcpy (&Ctrl->Float.data[nn], &tmp[(rowVec[mm]+n) * sizeof(float)], sizeof(float)); nn += incStep; @@ -1369,7 +1369,7 @@ int gmt_gdalread (struct GMT_CTRL *GMT, char *gdal_filename, struct GMT_GDALREAD case GDT_Float64: /* For now we don't care about doubles */ Ctrl->Float.active = true; for (m = startRow, mm = 0; m < endRow; m++, mm++) { - nn = (pad_w[piece]+m)*nXSize_withPad + startColPos; + nn = layer_offset + (pad_w[piece]+m)*nXSize_withPad + startColPos; for (n = 0; n < nXSize[piece]; n++) { double tmpF64; memcpy (&tmpF64, &tmp[(rowVec[mm]+n) * sizeof(double)], sizeof(double)); @@ -1389,6 +1389,7 @@ int gmt_gdalread (struct GMT_CTRL *GMT, char *gdal_filename, struct GMT_GDALREAD gmt_M_str_free (tmp); col_indent = nXSize[piece]; /* The second time (if there is one) we must step in to pick up western section */ indent = pad_w[piece] + nXSize[piece]; /* The second time (if there is one) we must step in to pick up western section */ + layer_offset += layer_size; /* As of now, this is never used after increment for float images */ } #if 0 /* This code is problematic and commented out for now. PW, 5/15/2016 */ if (Ctrl->Float.active && !isnan(prhs->N.nan_value)) { diff --git a/src/gmt_grdio.c b/src/gmt_grdio.c index 7fec344a14d..8fa096ea8e2 100644 --- a/src/gmt_grdio.c +++ b/src/gmt_grdio.c @@ -3434,13 +3434,19 @@ int gmt_raster_type (struct GMT_CTRL *GMT, char *file, bool extra) { break; } - if (extra && code == GMT_IS_GRID && (data[0] == 'I' || data[0] == 'M')) { /* Need to check more to determine if TIFF is grid or image */ + if (extra && code != GMT_IS_IMAGE && (data[0] == 'I' || data[0] == 'M')) { /* Need to check more to determine if TIFF is grid or image */ /* See if input could be an image of a kind that could also be a grid and we don't yet know what it is. Pass GMT_GRID_IS_IMAGE mode */ - struct GMT_GRID_HEADER_HIDDEN *HH = NULL; struct GMT_IMAGE *I = NULL; if ((I = GMT_Read_Data (GMT->parent, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY | GMT_GRID_IS_IMAGE, NULL, file, NULL)) != NULL) { - HH = gmt_get_H_hidden (I->header); /* Get pointer to hidden structure */ - if (I->header->n_bands > 1 || (HH->orig_datatype == GMT_UCHAR || HH->orig_datatype == GMT_CHAR)) code = GMT_IS_IMAGE; + struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (I->header); /* Get pointer to hidden structure */ + if (HH->pocket && strchr (HH->pocket, ',') == NULL) /* Got a single band request which we return as a grid */ + code = GMT_IS_GRID; + else if (HH->orig_datatype == GMT_UCHAR || HH->orig_datatype == GMT_CHAR) /* Got a gray or RGB image with or without transparency */ + code = GMT_IS_IMAGE; + else if (I->header->n_bands > 1) /* Whatever it is we must return multiband as an image */ + code = GMT_IS_IMAGE; + else /* Here we only have one band so it is a grid */ + code = GMT_IS_GRID; GMT_Destroy_Data (GMT->parent, &I); } } @@ -3631,7 +3637,7 @@ int gmtlib_read_image_info (struct GMT_CTRL *GMT, char *file, bool must_be_image gmt_M_free (GMT, from_gdalread); return (GMT_GRDIO_READ_FAILED); } - if (must_be_image && from_gdalread->band_field_names != NULL) { + if (from_gdalread->band_field_names != NULL) { /* Know the data type */ if (!strcmp(from_gdalread->band_field_names[0].DataType, "Byte")) I->type = GMT_UCHAR; else if (!strcmp(from_gdalread->band_field_names[0].DataType, "Int16")) @@ -3642,7 +3648,11 @@ int gmtlib_read_image_info (struct GMT_CTRL *GMT, char *file, bool must_be_image I->type = GMT_INT; else if (!strcmp(from_gdalread->band_field_names[0].DataType, "UInt32")) I->type = GMT_UINT; - else { + else if (!must_be_image && !strcmp(from_gdalread->band_field_names[0].DataType, "Float32")) + I->type = GMT_FLOAT; + else if (!must_be_image && !strcmp(from_gdalread->band_field_names[0].DataType, "Float64")) + I->type = GMT_FLOAT; /* No doubles for grids or images */ + else { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Using this data type (%s) is not implemented\n", from_gdalread->band_field_names[0].DataType); gmt_M_free (GMT, to_gdalread); @@ -3743,8 +3753,14 @@ int gmtlib_read_image (struct GMT_CTRL *GMT, char *file, struct GMT_IMAGE *I, do /* Tell gmt_gdalread that we already have the memory allocated and send in the *data pointer */ if (I->data) { /* Otherwise (subregion) memory is allocated in gdalread */ - to_gdalread->c_ptr.active = true; - to_gdalread->c_ptr.grd = I->data; + if (I->type == GMT_FLOAT || I->type == GMT_DOUBLE) { + to_gdalread->f_ptr.active = true; + to_gdalread->f_ptr.grd = (gmt_grdfloat *)I->data; + } + else { + to_gdalread->c_ptr.active = true; + to_gdalread->c_ptr.grd = I->data; + } } if (HH->grdtype > GMT_GRID_GEOGRAPHIC_LESS360) to_gdalread->R.periodic = true; /* Let gmt_gdalread know we have a periodic image */ if (gmt_gdalread (GMT, file, to_gdalread, from_gdalread)) { diff --git a/src/grdcut.c b/src/grdcut.c index 422e7183c8e..2445da48519 100644 --- a/src/grdcut.c +++ b/src/grdcut.c @@ -157,7 +157,6 @@ static int parse (struct GMT_CTRL *GMT, struct GRDCUT_CTRL *Ctrl, struct GMT_OPT bool F_or_R_or_J, do_file_check = true; unsigned int n_errors = 0, k, n_files = 0; - int ftype; char za[GMT_LEN64] = {""}, zb[GMT_LEN64] = {""}, zc[GMT_LEN64] = {""}, *c = NULL; struct GMT_OPTION *opt = NULL; struct GMTAPI_CTRL *API = GMT->parent; @@ -290,9 +289,16 @@ static int parse (struct GMT_CTRL *GMT, struct GRDCUT_CTRL *Ctrl, struct GMT_OPT n_errors += gmt_M_check_condition (GMT, !Ctrl->G.file && !Ctrl->D.active, "Option -G: Must specify output grid file\n"); n_errors += gmt_M_check_condition (GMT, n_files != 1, "Must specify one input grid file\n"); if (n_errors == 0) { - if ((ftype = gmt_raster_type (GMT, Ctrl->In.file, true)) == GMT_IS_IMAGE || ftype == GMT_IS_GRID) - Ctrl->In.type = ftype; - else /* Must assume it is a grid */ + int ftype = gmt_raster_type (GMT, Ctrl->In.file, true); + if (ftype == GMT_IS_IMAGE) /* Must read file as an image */ + Ctrl->In.type = GMT_IS_IMAGE; + else if (ftype == GMT_IS_GRID) { /* Check extension in case of special case */ + if (strstr (Ctrl->G.file, ".tif")) /* Want to write a single band (normally written as a grid) to geotiff instead */ + Ctrl->In.type = GMT_IS_IMAGE; + else + Ctrl->In.type = GMT_IS_GRID; + } + else /* Just have to assume it is a grid */ Ctrl->In.type = GMT_IS_GRID; if (Ctrl->In.type == GMT_IS_IMAGE) { n_errors += gmt_M_check_condition (GMT, Ctrl->Z.active, "Option -N: Cannot be used with an image\n"); @@ -409,7 +415,7 @@ EXTERN_MSC int GMT_grdcut (void *V_API, int mode, void *args) { int error = 0; unsigned int nx_old, ny_old, add_mode = 0U, side, extend, type = 0U, def_pad[4], pad[4]; uint64_t node; - bool outside[4] = {false, false, false, false}, all, bail = false, geo_to_cart = false; + bool outside[4] = {false, false, false, false}, all, bail = false, geo_to_cart = false, do_via_gdal = false; char *name[2][4] = {{"left", "right", "bottom", "top"}, {"west", "east", "south", "north"}}; @@ -451,10 +457,17 @@ EXTERN_MSC int GMT_grdcut (void *V_API, int mode, void *args) { if (!Ctrl->Z.active) { /* All of -F, -R, -S selections first needs the header */ if (Ctrl->In.type == GMT_IS_IMAGE) { GMT_Report (API, GMT_MSG_INFORMATION, "Processing input image\n"); - if ((I = GMT_Read_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) { + if ((I = GMT_Read_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY | GMT_GRID_IS_IMAGE, NULL, Ctrl->In.file, NULL)) == NULL) { Return (API->error); /* Get header only */ } h = I->header; + if ((I->type == GMT_FLOAT || h->n_bands > 4) && !gmt_M_file_is_memory (Ctrl->In.file) && !gmt_M_file_is_memory (Ctrl->G.file) && strstr (Ctrl->In.file, ".tif")) { + do_via_gdal = true; /* Use gdal_translate for this multi-layer data subset */ + if (!Ctrl->D.active && strstr (Ctrl->G.file, ".tif") == NULL && strstr (Ctrl->G.file, ".nc") == NULL && strstr (Ctrl->G.file, ".grd") == NULL) { + GMT_Report (API, GMT_MSG_INFORMATION, "Option -G: Must give an output file name with extensions .tiff (geotiff) or .nc or .grd (netCDF) when selecting multiband output from a geotiff file\n"); + Return (GMT_RUNTIME_ERROR); /* Get header only */ + } + } } else { GMT_Report (API, GMT_MSG_INFORMATION, "Processing input grid\n"); @@ -592,7 +605,7 @@ EXTERN_MSC int GMT_grdcut (void *V_API, int mode, void *args) { /* Note: The use of g and gmt_M_grd_row_to_y is correct since lon and lat args are not * coordinates computed from west or south in whole increments of dx dy. */ int row, col; - bool wrap; + bool wrap = gmt_M_360_range (h->wesn[XLO], h->wesn[XHI]); /* true if grid is 360 global */ if (gmt_M_is_cartesian (GMT, GMT_IN)) { GMT_Report (API, GMT_MSG_ERROR, "The -S option requires a geographic grid or image\n"); @@ -602,8 +615,10 @@ EXTERN_MSC int GMT_grdcut (void *V_API, int mode, void *args) { Return (GMT_NOT_A_VALID_TYPE); /* Set w/e to center and adjust in case of -/+ 360 stuff */ wesn_new[XLO] = wesn_new[XHI] = Ctrl->S.lon; - while (wesn_new[XLO] < h->wesn[XLO]) wesn_new[XLO] += 360.0, wesn_new[XHI] += 360.0; - while (wesn_new[XLO] > h->wesn[XHI]) wesn_new[XLO] -= 360.0, wesn_new[XHI] -= 360.0; + if (wrap) { + while (wesn_new[XLO] < h->wesn[XLO]) wesn_new[XLO] += 360.0, wesn_new[XHI] += 360.0; + while (wesn_new[XLO] > h->wesn[XHI]) wesn_new[XLO] -= 360.0, wesn_new[XHI] -= 360.0; + } wesn_new[YLO] = wesn_new[YHI] = Ctrl->S.lat; /* First adjust the S and N boundaries */ if (GMT->current.map.dist[GMT_MAP_DIST].arc) /* Got arc distance */ @@ -647,7 +662,6 @@ EXTERN_MSC int GMT_grdcut (void *V_API, int mode, void *args) { wesn_new[XHI] = h->wesn[XHI]; } else { /* Determine longitude limits */ - wrap = gmt_M_360_range (h->wesn[XLO], h->wesn[XHI]); /* true if grid is 360 global */ radius /= cosd (Ctrl->S.lat); /* Approximate e-w width in degrees longitude at center point */ wesn_new[XLO] -= radius; /* Approximate west limit in degrees */ if (!wrap && wesn_new[XLO] < h->wesn[XLO]) { /* Outside non-periodic grid range */ @@ -677,6 +691,7 @@ EXTERN_MSC int GMT_grdcut (void *V_API, int mode, void *args) { } wesn_new[XHI] = lon - (1.0 - h->xy_off) * h->inc[GMT_X]; /* Go one back since last col was outside */ } + if ((wesn_new[XHI] - wesn_new[XLO]) > 360.0) wesn_new[XHI] -= 360.0; /* One of them got off by 360 */ if (wesn_new[XHI] > 360.0) wesn_new[XLO] -= 360.0, wesn_new[XHI] -= 360.0; if (wesn_new[XHI] < 0.0) wesn_new[XLO] += 360.0, wesn_new[XHI] += 360.0; } @@ -759,6 +774,16 @@ EXTERN_MSC int GMT_grdcut (void *V_API, int mode, void *args) { if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_OFF) != GMT_NOERROR) { /* Enables data output and sets access mode */ Return (API->error); } + if (gmt_M_x_is_lon (GMT, GMT_IN)) { + if (wesn_new[XLO] < 0.0 && wesn_new[XHI] <= 0.0) + GMT->current.io.geo.range = GMT_IS_GIVEN_RANGE; + else if ((wesn_new[XHI] - wesn_new[XLO]) > 180.0) + GMT->current.io.geo.range = GMT_IS_GIVEN_RANGE; + else if (wesn_new[XLO] < 0.0 && wesn_new[XHI] >= 0.0) + GMT->current.io.geo.range = GMT_IS_M180_TO_P180_RANGE; + else + GMT->current.io.geo.range = GMT_IS_0_TO_P360_RANGE; + } if (Ctrl->D.text) { enum gmt_col_enum type = gmt_get_column_type (GMT, GMT_IN, GMT_X); char inc[GMT_LEN64] = {""}; @@ -844,6 +869,10 @@ EXTERN_MSC int GMT_grdcut (void *V_API, int mode, void *args) { GMT_Report (API, GMT_MSG_ERROR, "Old and new y_max do not differ by N * dy\n"); Return (GMT_RUNTIME_ERROR); } + if (!GMT->common.R.active[RSET]) { /* Got a new region indirectly via -S or -Z */ + gmt_M_memcpy (GMT->common.R.wesn, wesn_new, 4, double); + GMT->common.R.active[RSET] = true; + } gmt_grd_init (GMT, h, options, true); @@ -883,7 +912,7 @@ EXTERN_MSC int GMT_grdcut (void *V_API, int mode, void *args) { Return (API->error); } } - else { + else if (!do_via_gdal) { if (GMT_Read_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn_new, Ctrl->In.file, I) == NULL) { /* Get subset */ Return (API->error); } @@ -994,7 +1023,34 @@ EXTERN_MSC int GMT_grdcut (void *V_API, int mode, void *args) { } } else { /* Write an image */ - if (GMT_Write_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, I) != GMT_NOERROR) { + if (do_via_gdal) { /* Special case of multi-band geotiff */ + char driver[GMT_LEN128] = {""}, cmd[GMT_LEN256] = {""},the_band[GMT_LEN32] = {""}, *c = strchr (Ctrl->In.file, '='), *b = strstr (Ctrl->In.file, "+b"); + if (c) c[0] = '\0'; /* Chop off =gd request after that so that we only write the actual file name in the command */ + if (b) b[0] = '\0'; /* Chop offany band request after that so that we only write the actual file name in the command */ + if (wesn_new[XLO] > 180.0) wesn_new[XLO] -= 360.0, wesn_new[XHI] -= 360.0; /* GDAL expects -180/+180 */ + if (strstr (Ctrl->G.file, ".tif")) + sprintf (driver, "-of GTiff -co COMPRESS=DEFLATE"); + else + sprintf (driver, "-of netCDF -co COMPRESS=DEFLATE -co FORMAT=NC4 -co ZLEVEL=%d -a_nodata NaN", GMT->current.setting.io_nc4_deflation_level); + sprintf (cmd, "gdal_translate -projwin %.10lg %.10lg %.10lg %.10lg %s %s %s", wesn_new[XLO], wesn_new[YHI], wesn_new[XHI], wesn_new[YLO], driver, Ctrl->In.file, Ctrl->G.file); + if (c) c[0] = '='; /* Restore full file name */ + if (b) b[0] = '+'; /* Restore band requests */ + if (b) { /* Parse and add specific band request(s) to gdal_translate */ + char p[GMT_LEN64] = {""}; + unsigned int pos = 0, band; + while ((gmt_strtok (&b[2], ",", &pos, p))) { + band = atoi (p) + 1; /* We start counting at 0, GDAL at 1 */ + sprintf (the_band, " -b %d", band); + strcat (cmd, the_band); + } + } + GMT_Report (API, GMT_MSG_INFORMATION, "The gdal_translate command: \n%s\n", cmd); + if (system (cmd)) { /* Execute the gdal_translate command */ + GMT_Report (API, GMT_MSG_ERROR, "Error calling %s", cmd); + Return (GMT_RUNTIME_ERROR); + } + } + else if (GMT_Write_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, I) != GMT_NOERROR) { Return (API->error); } }