This workflow demonstrates a full pipeline:
- Search: Query a STAC catalog and store results in a temporary table.
- Variable Setup: Use
SET VARIABLE to create arrays of URLs for specific bands (e.g., nir_urls, red_urls). - Read & Mosaic: Use
RT_Read(getvariable(...)) to read all matching tiles as a single mosaic. - Clip & Join: Join the NIR and RED mosaics by their grid coordinates (
tile_x, tile_y), then use RT_CubeClip to clip the resulting cubes to the Area of Interest (AOI). - Export: Compute NDVI and export the result as a Cloud-Optimized GeoTIFF (COG).
-- Part A: Search and store results
CREATE OR REPLACE TEMP TABLE __search_results AS
-- ... (STAC search logic) ...
SELECT
replace(features->'assets'->'B08'->>'href', 'https://sentinel-cogs.s3.us-west-2.amazonaws.com', '/vsis3/sentinel-cogs') AS nir_url,
replace(features->'assets'->'B04'->>'href', 'https://sentinel-cogs.s3.us-west-2.amazonaws.com', '/vsis3/sentinel-cogs') AS red_url
FROM __features;
SET VARIABLE nir_urls = (SELECT list(nir_url) FROM __search_results);
SET VARIABLE red_urls = (SELECT list(red_url) FROM __search_results);
-- Part B: Read, clip, and compute NDVI
COPY (
WITH __aoi AS (SELECT geom FROM ST_Read('./test/data/CATAST_Pol_Township-PNA.geojson')),
__nir AS (
SELECT tile_x, tile_y, geometry, metadata, databand_1 AS nir
FROM RT_Read(getvariable('nir_urls'), blocksize_x := 1024, blocksize_y := 1024)
WHERE ST_Intersects(geometry, (SELECT geom FROM __aoi))
),
__red AS (
SELECT tile_x, tile_y, geometry, metadata, databand_1 AS red
FROM RT_Read(getvariable('red_urls'), blocksize_x := 1024, blocksize_y := 1024)
WHERE ST_Intersects(geometry, (SELECT geom FROM __aoi))
),
__data AS (
SELECT
n.tile_x, n.tile_y, n.geometry, n.metadata,
RT_CubeClip(n.nir, n.tile_x, n.tile_y, n.metadata, (SELECT geom FROM __aoi), (n.metadata->'bands'->0->'nodata')::DOUBLE) AS nir,
RT_CubeClip(r.red, r.tile_x, r.tile_y, r.metadata, (SELECT geom FROM __aoi), (r.metadata->'bands'->0->'nodata')::DOUBLE) AS red
FROM __nir n
JOIN __red r ON n.tile_x = r.tile_x AND n.tile_y = r.tile_y
)
SELECT geometry, RT_Cube2TypeFloat((nir - red) / (nir + red)) AS ndvi FROM __data
)
TO './test/data/sentinel2_to_ndvi.tiff'
WITH (
FORMAT RASTER, DRIVER 'COG', CREATION_OPTIONS ('COMPRESS=LZW'), RESAMPLING 'nearest', SRS 'EPSG:32630',
COMPUTE_VALID_ENVELOPE true, GEOMETRY_COLUMN 'geometry', DATABAND_COLUMNS ['ndvi']
);