INSTALL raster FROM community;
LOAD raster;

-- List available GDAL raster drivers and formats:
SELECT short_name, long_name, help_url FROM RT_Drivers();

┌────────────────┬──────────────────────────────────────────────────────────┬─────────────────────────────────────────────────────┐
│   short_name   │                        long_name                         │                      help_url                       │
│    varchar     │                         varchar                          │                       varchar                       │
├────────────────┼──────────────────────────────────────────────────────────┼─────────────────────────────────────────────────────┤
│ VRT            │ Virtual Raster                                           │ https://gdal.org/drivers/raster/vrt.html            │
│ GTiff          │ GeoTIFF                                                  │ https://gdal.org/drivers/raster/gtiff.html          │
│ COG            │ Cloud optimized GeoTIFF generator                        │ https://gdal.org/drivers/raster/cog.html            │
│  ·             │          ·                                               │                    ·                                │
│  ·             │          ·                                               │                    ·                                │
│  ·             │          ·                                               │                    ·                                │
│ ENVI           │ ENVI .hdr Labelled                                       │ https://gdal.org/drivers/raster/envi.html           │
│ Zarr           │ Zarr                                                     │ NULL                                                │
│ HTTP           │ HTTP Fetching Wrapper                                    │ NULL                                                │
└────────────────┴──────────────────────────────────────────────────────────┴─────────────────────────────────────────────────────┘

-- Read a raster file — one row per tile, one datacube column per band:
SELECT * FROM RT_Read('path/to/raster/file.tif');

┌───────┬───────────┬────────────┬────────────────────────────────┬─────────────────────────┬───────┬────────┬────────┬───────┬───────┬────────────┬────────────┐
│  id   │     x     │     y      │              bbox              │        geometry         │ level │ tile_x │ tile_y │ cols  │ rows  │  metadata  │ databand_1 │
│ int64 │  double   │   double   │ struct(xmin, ymin, xmax, ymax) │ geometry('epsg:25830')  │ int32 │ int32  │ int32  │ int32 │ int32 │    JSON    │    BLOB    │
├───────┼───────────┼────────────┼────────────────────────────────┼─────────────────────────┼───────┼────────┼────────┼───────┼───────┼────────────┼────────────┤
│     0 │ 545619.75 │ 4724508.25 │ {xmin: 545539.75, ...}         │ POLYGON ((...))         │     0 │      0 │      0 │   320 │     8 │ {...}      │ ...        │
│     1 │ 545619.75 │ 4724504.25 │ {xmin: 545539.75, ...}         │ POLYGON ((...))         │     0 │      0 │      1 │   320 │     8 │ {...}      │ ...        │
└───────┴───────────┴────────────┴────────────────────────────────┴─────────────────────────┴───────┴────────┴────────┴───────┴───────┴────────────┴────────────┘

-- Compute NDVI directly in SQL using band algebra:
WITH __input AS (
    SELECT
        databand_1 AS red,
        databand_2 AS nir
    FROM
        RT_Read('path/to/raster/file.tif', blocksize_x := 512, blocksize_y := 512)
)
SELECT
    RT_Cube2TypeFloat((nir - red) / (nir + red)) AS ndvi
FROM
    __input
;

-- Extract pixel values into SQL arrays for aggregate analysis:
WITH __input AS (
    SELECT RT_Cube2ArrayFloat(databand_1, true) AS band
    FROM RT_Read('path/to/raster/file.tif', blocksize_x := 512, blocksize_y := 512)
)
SELECT
    list_min(band.values)        AS band_min,
    list_stddev_pop(band.values) AS band_stddev,
    list_max(band.values)        AS band_max
FROM __input
;

-- Write a raster file using COPY ... FORMAT RASTER:
COPY (
    SELECT geometry, databand_1, databand_2, databand_3
    FROM RT_Read('path/to/raster/file.tif')
)
TO 'path/to/output/file.tif'
WITH (
    FORMAT RASTER,
    DRIVER 'COG',
    CREATION_OPTIONS ('COMPRESS=LZW'),
    RESAMPLING 'nearest',
    COMPUTE_VALID_ENVELOPE true,
    SRS 'EPSG:25830',
    GEOMETRY_COLUMN 'geometry',
    DATABAND_COLUMNS ['databand_3', 'databand_2', 'databand_1']
);

-- Clip raster tiles to a polygon geometry and write to a new file:
LOAD spatial;
LOAD json;

COPY (
    WITH __clip_layer AS (
        SELECT geom FROM ST_Read('./test/data/CATAST_Pol_Township-PNA.geojson')
    ),
    __dataset AS (
        SELECT id, x, y, tile_x, tile_y, cols, rows, geometry, databand_1, metadata
        FROM
            __clip_layer,
            RT_Read([
                './test/data/mosaic/SCL.tif-land-clip00.tiff',
                './test/data/mosaic/SCL.tif-land-clip01.tiff',
                './test/data/mosaic/SCL.tif-land-clip10.tiff',
                './test/data/mosaic/SCL.tif-land-clip11.tiff'
            ])
        WHERE ST_Intersects(__clip_layer.geom, geometry)
    ),
    __clipped AS (
        SELECT
            RT_CubeClip(databand_1,
                        tile_x,
                        tile_y,
                       (metadata->'blocksize_x')::INTEGER,
                       (metadata->'blocksize_y')::INTEGER,
                       (metadata->'transform')::DOUBLE[],
                        __clip_layer.geom,
                       (metadata->'bands'->0->'nodata')::DOUBLE
            ) AS databand,
            * EXCLUDE(databand_1)
        FROM __clip_layer, __dataset
    )
    SELECT geometry, databand FROM __clipped
)
TO 'path/to/output/clipped.tiff'
WITH (
    FORMAT RASTER,
    DRIVER 'GTiff',
    CREATION_OPTIONS ('COMPRESS=LZW'),
    RESAMPLING 'nearest',
    SRS 'EPSG:32630',
    GEOMETRY_COLUMN 'geometry',
    DATABAND_COLUMNS ['databand']
);