Skip to content

Let multi-layer geotiff files be cut via gdal_translate #5819

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
PaulWessel merged 11 commits into from
Oct 5, 2021
Merged

Let multi-layer geotiff files be cut via gdal_translate #5819

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
3d4596c
Let multi-layer geotiff files be cut via gdal_translate
PaulWessel Sep 29, 2021
4f96621
Merge branch 'master' into grdcut-gdal
PaulWessel Oct 1, 2021
662bd7a
Add compression
PaulWessel Oct 1, 2021
e6963d8
Handle band requests
PaulWessel Oct 1, 2021
afaf086
Document and add checks
PaulWessel Oct 1, 2021
a0133e8
Fix so -S also works for images
PaulWessel Oct 2, 2021
a8c0ca5
Handle netCDF output via GDAL as well
PaulWessel Oct 2, 2021
cdb134f
Use IO_NC4_DEFLATION_LEVEL to set gdal_translate compression level
PaulWessel Oct 2, 2021
0b0d641
Update grdcut.c
PaulWessel Oct 2, 2021
d91c185
Merge branch 'master' into grdcut-gdal
PaulWessel Oct 4, 2021
8c25cce
Update grdcut.c
PaulWessel Oct 5, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 24 additions & 8 deletions doc/rst/source/grdcut.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
------------------
Expand Down Expand Up @@ -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
--------

Expand Down
8 changes: 6 additions & 2 deletions src/gmt_api.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
21 changes: 11 additions & 10 deletions src/gmt_gdalread.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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));
Expand All @@ -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)) {
Expand Down
32 changes: 24 additions & 8 deletions src/gmt_grdio.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
}
Expand Down Expand Up @@ -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"))
Expand All @@ -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);
Expand Down Expand Up @@ -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)) {
Expand Down
Loading