From 228999ef4f3c4b73def8bd272de2ebd7617e8030 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Thu, 30 Mar 2023 19:40:19 +0200
Subject: [PATCH] COG: add NBITS creation option (fixes #7361)

---
 autotest/gcore/cog.py             | 21 +++++++++++++++++++++
 doc/source/drivers/raster/cog.rst |  6 ++++++
 frmts/gtiff/cogdriver.cpp         |  4 ++++
 3 files changed, 31 insertions(+)

diff --git a/autotest/gcore/cog.py b/autotest/gcore/cog.py
index 8bc1ffb3daaa..90a2d9e382c0 100755
--- a/autotest/gcore/cog.py
+++ b/autotest/gcore/cog.py
@@ -1702,3 +1702,24 @@ def test_cog_write_jpegxl_alpha_distance_zero():
     ds = None
 
     gdal.Unlink(filename)
+
+
+###############################################################################
+# Test NBITS creation option
+
+
+def test_cog_NBITS():
+
+    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1)
+    drv = gdal.GetDriverByName("COG")
+    filename = "/vsimem/test_cog_NBITS.tif"
+    drv.CreateCopy(
+        filename,
+        src_ds,
+        options=["NBITS=7"],
+    )
+    ds = gdal.Open(filename)
+    assert ds.GetRasterBand(1).GetMetadataItem("NBITS", "IMAGE_STRUCTURE") == "7"
+    ds = None
+
+    gdal.Unlink(filename)
diff --git a/doc/source/drivers/raster/cog.rst b/doc/source/drivers/raster/cog.rst
index 73c138ae2635..662b31dfd0e6 100644
--- a/doc/source/drivers/raster/cog.rst
+++ b/doc/source/drivers/raster/cog.rst
@@ -113,6 +113,12 @@ General creation options
    or TARGET_SRS creation options. (Overview generation is also multithreaded since
    GDAL 3.2)
 
+-  **NBITS=n**: (GDAL >= 3.7) Create a file with less than 8 bits per sample by
+   passing a value from 1 to 7. The apparent pixel type should be Byte.
+   Values of n=9...15 (UInt16 type) and n=17...31
+   (UInt32 type) are also accepted. From GDAL 2.2, n=16 is accepted for
+   Float32 type to generate half-precision floating point values.
+
 -  **PREDICTOR=[YES/NO/STANDARD/FLOATING_POINT]**: Set the predictor for LZW,
    DEFLATE and ZSTD compression. The default is NO. If YES is specified, then
    standard predictor (Predictor=2) is used for integer data type,
diff --git a/frmts/gtiff/cogdriver.cpp b/frmts/gtiff/cogdriver.cpp
index e501832b126c..fd11920bb1d3 100644
--- a/frmts/gtiff/cogdriver.cpp
+++ b/frmts/gtiff/cogdriver.cpp
@@ -1136,6 +1136,7 @@ GDALDataset *GDALCOGCreator::Create(const char *pszFilename,
                             CSLFetchNameValue(papszOptions, "GEOTIFF_VERSION"));
     aosOptions.SetNameValue("SPARSE_OK",
                             CSLFetchNameValue(papszOptions, "SPARSE_OK"));
+    aosOptions.SetNameValue("NBITS", CSLFetchNameValue(papszOptions, "NBITS"));
 
     if (EQUAL(osOverviews, "NONE"))
     {
@@ -1390,6 +1391,9 @@ void GDALCOGDriver::InitializeCreationOptionList()
         "   <Option name='NUM_THREADS' type='string' "
         "description='Number of worker threads for compression. "
         "Can be set to ALL_CPUS' default='1'/>"
+        "   <Option name='NBITS' type='int' description='BITS for sub-byte "
+        "files (1-7), sub-uint16_t (9-15), sub-uint32_t (17-31), or float32 "
+        "(16)'/>"
         "   <Option name='BLOCKSIZE' type='int' "
         "description='Tile size in pixels' min='128' default='512'/>"
         "   <Option name='BIGTIFF' type='string-select' description='"