Paging @joa-quim on this one: Now that grdcut can return a subset of an image (your typical 1, 3, or 4 band byte data) the issue becomes can grdcut also do the same for more general "images" where there are numerous bands and they are stored as floats. GDAL has a solution for this, e.g.,
gdal_translate -projwin -163.051107059 61.219982021 -163.047286129 61.217258591 -of GTiff sentinel_mosaic.tif sentinel_mosaic_clipped.tif
where the sentinel file has 16 double-precision bands, but it is unaware of periodic boundary conditions. Perhaps there are no large global "images" like this we need to worry about, but if there are it would be nice for GMT to handle that since it now can do periodic boundaries for images. Alternatively, we could implement a call to gdal_translate, like we do in gmt_remote.c to convert JP2 files to netCDF. A limitation would be passing the image around in memory.
There seem to be some unimplemented situations in gmt_gdalread|write.c that would need some updates for that to work. So far these are:
- When a float image is encountered, we successfully learn about it in the image header, but the data type, if not an int or byte, is not stored in the header, so it says GMT_CHAR. This means we will allocate a byte array to hold a stack of double precision layers and this will not work.
- It also seems that gmt_gdalwrite.c makes assumptions of layers being byte so things will go wrong here as well.
Other than that, it seems the machinery is set up to do n_bands just fine. So perhaps overcoming these obstacles is mostly about setting image->type to GMT_DOUBLE (for this data set) and use that when allocating the output array as well as the de-padding step in gmt_gdalwrite.c
Paging @joa-quim on this one: Now that grdcut can return a subset of an image (your typical 1, 3, or 4 band byte data) the issue becomes can grdcut also do the same for more general "images" where there are numerous bands and they are stored as floats. GDAL has a solution for this, e.g.,
gdal_translate -projwin -163.051107059 61.219982021 -163.047286129 61.217258591 -of GTiff sentinel_mosaic.tif sentinel_mosaic_clipped.tifwhere the sentinel file has 16 double-precision bands, but it is unaware of periodic boundary conditions. Perhaps there are no large global "images" like this we need to worry about, but if there are it would be nice for GMT to handle that since it now can do periodic boundaries for images. Alternatively, we could implement a call to gdal_translate, like we do in gmt_remote.c to convert JP2 files to netCDF. A limitation would be passing the image around in memory.
There seem to be some unimplemented situations in gmt_gdalread|write.c that would need some updates for that to work. So far these are:
Other than that, it seems the machinery is set up to do n_bands just fine. So perhaps overcoming these obstacles is mostly about setting image->type to GMT_DOUBLE (for this data set) and use that when allocating the output array as well as the de-padding step in gmt_gdalwrite.c