Create a median composite raster, using Sentinel-2 data, comprised of red, green and blue bands over a specifed date range. The result will be clipped to Bengaluru in India.

This tutorial achieves the same result as Ujaval Gandhi's Python tutorial, but using only GDAL and Bash, with no dependencies. Ujaval has a wealth of information and definitely worth the reading.

The basic premise: if a process is complex to learn as a student, it will be equally difficult to maintain as a professional. Keep it simple. Don't use more code than required. Use the right tools, not the currently popular ones.

To reiterate, this is GDAL and Bash only. Nothing else. Bash is already part of your Linux default environment. If developing for GIS, GDAL is probably installed. Python and it's unwieldy environment and dozens of packages/dependencies is not needed. Why bother with it at all? The following is about 1/3 the size of the Python tutorial.

You probably already have GDAL installed, as almost every proprietary/open-source spatial software uses it. You will need GDAL version 3.13. Let's begin.

# Let's set some Bash variables:
# A prefix for naming resulting TIF, PNG, GPKG files.
f="Bengaluru"

# Date range of Sentinel raster bands. This will result in 8 source tif's for each red, green, blue.
range='2023-04-01T00:00:00Z/2023-05-31T23:59:59Z'

# Source bounding box generally covering Bengaluru to speed up remote queries
srcbbox="61.966,5.051,96.595,33.946"

# where clause for SQL, to get exactly the clipping polygon we need for Bengaluru
where="country = 'IN' and region = 'IN-KA' and \"names.primary\" like 'Bengaluru' and is_land = 1"

# Maximum cloud cover
cloudcover="3"

# Output of clipping polygon
poly="aoi-${f}.gpkg"

# Get the Bengaluru polygon poly from overture
gdal vector pipeline \
   ! read -i /vsis3/overturemaps-us-west-2/release/2026-06-17.0/'theme=divisions/type=division_area' \
   ! filter --bbox ${srcbbox} --where "${where}" \
   ! clip --bbox ${srcbbox} --bbox-crs EPSG:4326 \
   ! write $poly --output-layer aoi --overwrite

# get bbox from poly extent. We'll use bbox throughout
bbox="$(gdal vector info $poly --of json|jq -r '.layers[].geometryFields[].extent|@csv')"

# get sentinel meta data constrained by bbox and date range
gdal vector convert \
   -i "https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a&datetime=${range}&bbox=${bbox}&limit=500" \
   -o ${f}.geojson --overwrite

# SQL template for querying meta data we captured from above, creating a list of all assets intersecting above poly
sql="select group_concat('[ ! read /vsicurl/' || \"assets.%s.href\" || ' ! reproject --bbox %s --bbox-crs EPSG:4326 --output-crs EPSG:4326 ]', ' ') from %s where \"eo:cloud_cover\" < %s"

# download and save our 3 tif's. Each resulting file, ie red.tif, will have 8 bands 
for band in $(echo "red green blue") ; do
   gdal raster pipeline ! stack --resolution highest $(gdal vector sql -i ${f}.geojson --dialect sqlite --lco HEADER=NO --sql "$(printf "$sql" "${band}" "$bbox" $f $cloudcover)" --of CSV -o /vsistdout/|tr -d '"') ! write -o ${band}.tif --co COMPRESS=DEFLATE --overwrite
done

# If you do a 'gdal raster info' on any of these TIF's you note, offset/scale defined, which we use in --calc:
#   Offset: -0.1,   Scale:0.0001

# This is a GDAL pipeline that performs several steps:
# 1. get the median value for all bands in, ie, red.tif. --calc can be any mathematic forumula we want. Median is pre-defined by GDAL in this case.
# 2. stack resulting single median band for each color, resulting in a composited RGB raster
# 3. clip that raster to the Bengaluru polygon
# 4. save as a Cloud Optimized GeoTiff (COG)
gdal raster pipeline \
   ! stack \
      [ ! calc -i "r=red.tif" --calc="avg(r) * .1 + .0001" --ot Byte --nodata 255 --flatten ] \
      [ ! calc -i "r=green.tif" --calc="avg(r) * .1 + .0001" --ot Byte --nodata 255 --flatten ] \
      [ ! calc -i "r=blue.tif" --calc="avg(r) * .1 + .0001" --ot Byte --nodata 255 --flatten ] \
   ! clip --like $poly --allow-bbox-outside-source \
   ! write -o ${f}.tif --co COMPRESS=DEFLATE --of COG --overwrite

# ...and for optimal viewing, create a PNG. (image at top of page)

gdal raster convert -i ${f}.tif -o ${f}.png --overwrite

Here's the entire script without comments:

#!/usr/bin/env bash

f="Bengaluru"
range='2023-04-01T00:00:00Z/2023-05-31T23:59:59Z'
srcbbox="61.966,5.051,96.595,33.946"
where="country = 'IN' and region = 'IN-KA' and \"names.primary\" like 'Bengaluru' and is_land = 1"
cloudcover="3"
poly="aoi-${f}.gpkg"

gdal vector pipeline \
   ! read -i /vsis3/overturemaps-us-west-2/release/2026-06-17.0/'theme=divisions/type=division_area' \
   ! filter --bbox ${srcbbox} --where "${where}" \
   ! clip --bbox ${srcbbox} --bbox-crs EPSG:4326 \
   ! write $poly --output-layer aoi --overwrite

bbox="$(gdal vector info $poly --of json|jq -r '.layers[].geometryFields[].extent|@csv')"

gdal vector convert \
   -i "https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a&datetime=${range}&bbox=${bbox}&limit=500" \
   -o ${f}.geojson --overwrite

sql="select group_concat('[ ! read /vsicurl/' || \"assets.%s.href\" || ' ! reproject --bbox %s --bbox-crs EPSG:4326 --output-crs EPSG:4326 ]', ' ') from %s where \"eo:cloud_cover\" < %s"

for band in $(echo "red green blue") ; do
   gdal raster pipeline ! stack --resolution highest $(gdal vector sql -i ${f}.geojson --dialect sqlite --lco HEADER=NO --sql "$(printf "$sql" "${band}" "$bbox" $f $cloudcover)" --of CSV -o /vsistdout/|tr -d '"') ! write -o ${band}.tif --co COMPRESS=DEFLATE --overwrite
done

gdal raster pipeline \
   ! stack \
      [ ! calc -i "r=red.tif" --calc="avg(r) * .1 + .0001" --ot Byte --nodata 255 --flatten ] \
      [ ! calc -i "r=green.tif" --calc="avg(r) * .1 + .0001" --ot Byte --nodata 255 --flatten ] \
      [ ! calc -i "r=blue.tif" --calc="avg(r) * .1 + .0001" --ot Byte --nodata 255 --flatten ] \
   ! clip --like $poly --allow-bbox-outside-source \
   ! write -o ${f}.tif --co COMPRESS=DEFLATE --overwrite 
   gdal raster convert -i ${f}.tif -o ${f}.png --overwrite

That's it!