Useful gdal commands

gdal_commands.md

1. Rescale specfic bands and write them to output tiff

gdal_translate -ot Byte -b 5 -b 3 -b 2 -scale_1 0 1970 0 255 -scale_2 0 1860 0 255 -scale_3 0 1174 0 255 WV.TIF WV_RGB.TIF

We take bands 5 (red), 3 (green) and 2(blue), scale them to 8 bit range and write it to an output file. Things to note:

a. Output scale by default is 0 to 255, so there's no real reason to specify. 
b. Once we pick the bands, they are thereafter referred w.r.t output. This is why we specify `scale_1` to rescale band 5 because band 5 in the input is band 1 of output.

2. Converting GeoTiff to PDF

The process is pretty simple as seen below

gdal_translate -of PDF in.tif out.pdf

But conversion to PDF only works if your input tiff is an 8 bit 1/3/4 band image.

3. Merge large number of tiffs

You can use gdalwarp or gdal_merge to merge geotiffs together. Something to note: gdal_merge uses nearest neighbour resampling and it's recommended to use gdalwarp if you want more control over the resampling used. Refer here

For merging a really large number of images or a number of really large images, you should not use gdal_merge since it loads all the images to memory and you might end up with a MemoryError. To circumvent this, use vrt as follows:

1. Create VRT file like this: gdalbuildvrt merged.vrt *.jp2
2. gdal_translate -of GTiff -ot Byte -co PHOTOMETRIC=RGB -co COMPRESS=DEFLATE merged.vrt merged.tif

Ref: https://gis.stackexchange.com/questions/306274/merge-thousands-of-rasters-with-gdal-merge-py

4. Reproject tiffs in src_dir into dst_dir

Shell script for quiet reprojection

for f in src_dir/*.tif*; do gdalwarp -q -t_srs EPSG:4326 "$f" dst_dir/"${f##*/}"; done

Alternative python script

from osgeo import gdal
from pathlib import Path
from tqdm import tqdm

img_dir = "src_dir_path"
out_dir = "dst_dir_path"

img_dir = Path(img_dir)
out_dir = Path(out_dir)

imgs = [im for im in img_dir.iterdir() if im.suffix.lower() in [".tiff", ".tif"]]
epsg = "EPSG:4326"

for img in tqdm(imgs):
    out_file = str(out_dir / img.name)
    input_raster = gdal.Open(str(img))
    gdal.Warp(out_file, input_raster, dstSRS=epsg)

5. Formula for calculating memory size.

If your image resolution is 250m and you have a 32 bit image

memory = pixels * bits per pixel
memory = (area of earth / area of cell) * bits
memory = (510.1 trillion m² / (250 m * 250 m)) * 32 bits = 32.65 GB