[Raster.write] add support for COG driver

Closes #24

Here is a comparison between a COG produced by the changes I implemented (COG_1) and the workaround that was implemented in COPA until now (COG_2):

ref = 's1a-iw-nrb-20210514t171503-20210514t171555-037887-0478b1-32tns-vh-g-lin.tif')
tmp_tif = 'tmp.tif'
cog_1 = 'cog_1.tif'
cog_2 = 'cog_2.tif'

overviews = [2, 4, 8, 16, 32]
options_1 = ['BLOCKSIZE=512', 'COMPRESS=LERC_ZSTD', 'MAX_Z_ERROR=2e-5']
options_2 = ['BLOCKXSIZE=512', 'BLOCKYSIZE=512', 'TILED=YES', 'INTERLEAVE=BAND', 'COMPRESS=LERC_ZSTD', 'MAX_Z_ERROR=2e-5', 'COPY_SRC_OVERVIEWS=YES']

# COG_1
with Raster(ref) as ras:
    out = np.random.rand(ras.cols, ras.rows)
    out[np.isnan(ras.array())] = np.nan
    ras.write(cog_1, format='COG', array=out.astype('float32'), nodata=0, overviews=overviews, options=options_1)

# COG_2
with Raster(ref) as ras:
    out = np.random.rand(ras.cols, ras.rows)
    out[np.isnan(ras.array())] = np.nan
    ras.write(tmp_tif, format='GTiff', array=out.astype('float32'), nodata=0)

    if overviews is not None:
        raster = gdal.Open(tmp_tif)
        raster.BuildOverviews('NEAREST', overviews)
        raster = None

    gdal_translate(tmp_tif, cog_2,
                   options={'format': 'GTiff', 'creationOptions': options_2})

Check COG validity:

>python dev_validate_cog.py 'cog_1.tif' 
cog_1.tif is a valid cloud optimized GeoTIFF

The size of all IFD headers is 6794 bytes

>python dev_validate_cog.py 'cog_2.tif' 
cog_2.tif is a valid cloud optimized GeoTIFF

The size of all IFD headers is 6794 bytes

gdalinfo output only once, because they're identical:

Driver: GTiff/GeoTIFF
Files: cog_1.tif
Size is 10980, 10980
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 32N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State."],
        BBOX[0,6,84,12]],
    ID["EPSG",32632]]
Data axis to CRS axis mapping: 1,2
Origin = (499980.000000000000000,5200020.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LERC_ZSTD
  INTERLEAVE=BAND
  LAYOUT=COG
  LERC_VERSION=2.4
Corner Coordinates:
Upper Left  (  499980.000, 5200020.000) (  8d59'59.05"E, 46d57'13.35"N)
Lower Left  (  499980.000, 5090220.000) (  8d59'59.07"E, 45d57'55.98"N)
Upper Right (  609780.000, 5200020.000) ( 10d26'33.03"E, 46d56'40.64"N)
Lower Right (  609780.000, 5090220.000) ( 10d24'59.96"E, 45d57'24.37"N)
Center      (  554880.000, 5145120.000) (  9d42'52.78"E, 46d27'26.71"N)
Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
  NoData Value=0
  Overviews: 5490x5490, 2745x2745, 1373x1373, 687x687, 344x344

I also checked with gdalcompare, but I'm not sure about the results to be honest:

>gdalcompare.py cog_2.tif cog_1.tif
Files differ at the binary level.
Band 1 checksum difference:
  Golden: 28172
  New:    28157
  Pixels Differing: 120560393
  Maximum Pixel Difference: 0
Band 1 overview 0 checksum difference:
  Golden: 37221
  New:    34881
  Pixels Differing: 30140099
  Maximum Pixel Difference: 0
Band 1 overview 1 checksum difference:
  Golden: 29967
  New:    29484
  Pixels Differing: 7535025
  Maximum Pixel Difference: 0
Band 1 overview 2 checksum difference:
  Golden: 9834
  New:    8066
  Pixels Differing: 1885129
  Maximum Pixel Difference: 0
Band 1 overview 3 checksum difference:
  Golden: 36259
  New:    34378
  Pixels Differing: 471969
  Maximum Pixel Difference: 0
Band 1 overview 4 checksum difference:
  Golden: 40540
  New:    41706
  Pixels Differing: 118336
  Maximum Pixel Difference: 0.9647746048867702
Differences Found: 7

johntruckenbrodt / spatialist

[Raster.write] add support for COG driver #25