Geospatial Computing from the Command Line

Installation

The commands on this page use gdal and jq.

Conda

You can install these packages on any operating system using conda/mamba.

mamba create -n my-gdal-env python=3 gdal jq

Then activate the environment every time you want to work with gdal.

mamba activate my-gdal-env

Linux

On Ubuntu these can be installed using:

sudo apt install gdal-bin jq

Get information about a raster

gdalinfo provides information about raster files.

gdalinfo myraster.tif will produce a basic readable output to the screen.

This output can also be written to JSON

gdalinfo -json myraster.tif

Writing to JSON makes it easy to use individual pieces, e.g., to look up the dimensions of the raster (you’ll need to install jq to do this).

gdalinfo -json myraster.tif | jq -r .size

or just the width of the raster

gdalinfo -json myraster.tif | jq -r .size[0]

Splitting rasters using gdal

One way to split a raster into pieces is to use the gdal_retile.py Python script bundled with gdal.

The following command will split myraster.tif into 1500x1500 pixel rasters stored in outputdir. The first number is the width (in pixels) and the second is the height (in pixels) of each chunk.

gdal_retile.py -ps 1500 1500 -targetDir outputdir myraster.tif

Files will be labeled with _row_col and so if the original image was 4500x1500 then the above command would produce three output files:

myraster_1_1.tif
myraster_2_1.tif
myraster_3_1.tif

Representing the top of the original raster (_1_1), the middle of the original raster (_2_1), and the bottom of the original raster (_3_1).

Split raster into horizontal strips

Our most common usage is to split large rasters into horizontal strips with manageable file sizes (< 3 GB). This can be automated by changing myraster.tif to the location of your raster and outputdir to the directory you want the split raster pieces stored in and running the code below:

RASTER=myraster.tif
OUTPUTDIR=outputdir
WIDTH=$(gdalinfo -json $RASTER | jq -r .size[0])
WIDTHPAD=$((WIDTH + 10)) # Padding prevents periodic inclusion of single pixel strip 
HEIGHT=$(expr 1000000000 / $WIDTH)
gdal_retile.py -ps $WIDTHPAD $HEIGHT -targetDir $OUTPUTDIR $RASTER

Combining/merging rasters using gdal

One way to combine rasters is to use the gdal_merge.py Python script bundled with gdal. We use LZW compression to reduce file sizes while ensuring that the resulting GeoTIFF can be used in all geospatial computing systems.

gdal_merge.py -o output_file.tif input_file_1.tif input_file_2.tif input_file_3.tif
gdal_translate -co COMPRESS=LZW -co PREDICTOR=2 -co BIGTIFF=YES output_file.tif compressed_file.tif

We use this two command approach instead of including compression in the merge command because gdal_merge.py doesn’t currently support BigTIFF creation correctly and many of our combined files are greater than the 4GB maximum for regular TIFFs. PREDICTOR=2 produces more efficient compression in the presence of spatial autocorrelation, which we often have.

Removing alpha channels

All of our code works with 3 band RGB rasters. Occasionally we accidentally produce a raster that contains a 4th alpha channel. This can be removed using GDAL.

First check to make sure the bands you want are the first three bands (they pretty much always are):

gdalinfo four_band_ortho.tif

This should show something like the following info about channels:

Band 1 Block=1501x1 Type=Byte, ColorInterp=Red
Band 2 Block=1501x1 Type=Byte, ColorInterp=Green
Band 3 Block=1501x1 Type=Byte, ColorInterp=Blue
Band 4 Block=1501x1 Type=Byte, ColorInterp=Alpha

The use gdal_translate.py to just keep the first three bands:

gdal_translate -b 1 -b 2 -b 3 four_band_ortho.tif three_band_ortho.tif

Original source

Previous
Next