Overview*¶
- about + purpose
- functionality)
- installation + basic use
- resources + tutorials
*focus on raster functionality
About + Purpose¶
Very Brief History¶
- since 1998 --- intitial release in 2000 (Frank Warmerdam)
- taken over by the Open Source Geospatial Foundation after v1.3.2 (ca. 2006)
- previously separate from OGR Simple Features Library, now more tightly coupled
- OGR refers to the vector data model and functionality
- often referred to as GDAL/OGR
- currently v2.3.2 (September 2018)
About¶
- X/MIT License (permissive free software license)
- C/C++ based, bindings for python, R, Java, Perl, C#, Ruby …
- cross-platform geospatial library (MacOS, Linux/Unix, Windows, Android, iOS)
(Intended) Purpose¶
- library to read/write raster and vector geospatial formats
- more than 200 supported geospatial data formats
- presents a single abstract data model for both raster and vector data formats, respectively
- provides various command-line utilities for translation and geospatial processing
Not a Geographic Information System!
Some Supported Raster Data Formats (of over 155)¶
- TerraSAR-X Complex SAR Data Product
- ENVI .hdr Labelled Raster
- GeoPackage
- GRASS Raster Format
- GeoTiff
- JPEG2000
- GDAL VRT
- Geospatial PDF
- Rasdaman
- Sentinel 1 SAFE
- Sentinel 2
- ArcSDE Raster
Some Software Utilising GDAL¶
Functionality¶
Some GDAL Utilities¶
gdalinfo - Report information about a file.
gdal_translate - Copy a raster file, with control of output format.
gdalwarp - Warp an image into a new coordinate system.
gdaldem - Tools to analyze and visualize DEMs.
gdal_merge - Build a quick mosaic from a set of images.
gdal_rasterize - Rasterize vectors into raster file.
gdal_polygonize - Generate polygons from raster.
gdal_pansharpen - Perform a pansharpen operation.
gdalcompare - Compare two images and report on differences.
GDAL Utilities docs here.
Installation¶
Python Bindings¶
- Recommended: the python virtual environment, package and dependency manager
conda- install Anaconda or miniconda (simplified version)
- start and check installation
- manage environments including different versions of python
- search and find different versions of packages
- create an environment with GDAL and whatever else you require
- activate the environment
- Go!
(Tutorial: Getting Started with conda)
Basic Use¶
Set-up necessary python imports¶
(the %matplotlib inline is a jupyter notebook command to show matplotlib plots in the output)
In [1]:
%matplotlib inline
import os
import sys
import numpy
import matplotlib.pyplot as plt
from osgeo import gdal
Check Versions¶
In [2]:
# Check version of python:
print(sys.version)
# We can check which version we're running by printing the "__version__" variable
print("GDAL's version is: " + gdal.__version__)
# This shows where the version of GDAL is located (in this case, in a conda environment)
print(gdal)
# This allows GDAL to throw Python Exceptions
gdal.UseExceptions()
In [4]:
info_options = gdal.InfoOptions(computeMinMax=True, stats=True)
print(gdal.Info(LC08band2_path, options=info_options))
Open and Close a Dataset in Python¶
The gdal.Open() function creates a GDAL Dataset object based on the given file.
In [5]:
# Open the Landsat 8 Band 2 file as a GDAL Dataset object.
in_ds = gdal.Open(LC08band2_path, gdal.GA_ReadOnly)
# Close GDAL Dataset object.
in_ds = None
Access Metadata from Dataset Object in Python¶
In [6]:
in_ds = gdal.Open(LC08band2_path, gdal.GA_ReadOnly)
# Find the number of bands in the dataset object.
print("[ NUMBER OF BANDS ] = ", in_ds.RasterCount)
In [7]:
# Create a GDAL raster band object from the only band in the dataset object.
in_band = in_ds.GetRasterBand(1)
print("[ BAND DATA TYPE ] = ", gdal.GetDataTypeName(in_band.DataType))
print("[ NO DATA VALUE ] = ", in_band.GetNoDataValue())
In [9]:
# Compute statistics if needed.
if in_band.GetMinimum() is None or in_band.GetMaximum()is None:
in_band.ComputeStatistics(0)
print("Statistics computed.")
print("[ OVERALL STATS ] = ", in_band.GetMetadata())
In [10]:
print("[ MIN ] = ", in_band.GetMinimum())
print("[ MAX ] = ", in_band.GetMaximum())
print("[ SCALE ] = ", in_band.GetScale())
# Close everything.
in_ds = None
in_band = None
Create Copies in Different Formats¶
CreateCopy() or Translate() functions in python, gdal_translate in the command-line.
In [11]:
# Check contents of Landsat 8 processed data folder.
print(os.listdir(LC08proc_path))
In [12]:
in_ds = gdal.Open(LC08band2_path, gdal.GA_ReadOnly)
# Define path to save copy and with what driver.
dst_filename = os.path.join(LC08proc_path, "LC08_B2.vrt")
driver = gdal.GetDriverByName("VRT")
dst_ds = driver.CreateCopy(dst_filename, in_ds, strict=0)
print("VRT copy of Band 2 has been created.")
dst_filename = os.path.join(LC08proc_path, "LC08_B2.dat")
driver = gdal.GetDriverByName("ENVI")
dst_ds = driver.CreateCopy(dst_filename, in_ds, strict=0)
print("ENVI copy of Band 2 has been created.")
# Close initial and resulting datasets.
dst_ds = None
in_ds = None
In [13]:
print(os.listdir(LC08proc_path))
Access Sentinel-2 bands as dictionary.¶
Get Sentinel-2 bands in a sorted order accessible via key-value pairs.
In [14]:
# Create dictionary with key-value pairs for easier access to bands.
S2_bands = {}
for entry in os.listdir(S2_path):
if (entry.endswith('.jp2')):
# This makes bands accessible by name (e.g. B01, B8A).
band_name = entry[-7:-4]
S2_bands[band_name]=os.path.join(S2_path, entry)
# Using the dictionary, individual bands can be accessed so:
print("Band 8A is called: " + os.path.basename(S2_bands['B8A']))
del band_name
Open True Colour Band and view with Matplotlib¶
In [15]:
img = gdal.Open(S2_bands['TCI'], gdal.GA_ReadOnly)
# Read GDAL Dataset Object as a Numpy Array
RGB_Image = img.ReadAsArray()
print(RGB_Image.shape)
# Get extent and boundary
ulx, xres, xskew, uly, yskew, yres = img.GetGeoTransform()
lrx = ulx + (img.RasterXSize * xres)
lry = uly + (img.RasterYSize * yres)
print(ulx, lrx, lry, uly)
# Re-order dimensions for matplotlib.
RGB_Image = numpy.dstack(((RGB_Image[0]), (RGB_Image[1]), (RGB_Image[2]))).astype('uint8')
print(RGB_Image.shape)
In [16]:
# Plot image in Matplotlib.
plt.figure(figsize=(10, 10))
plt.imshow(RGB_Image, extent=(ulx, lrx, lry, uly))
Out[16]:
In [17]:
# Optional -- Save figure.
# plt.savefig('TestImage.png', format='png', dpi=300, bbox_inches='tight')
# Clean Up.
plt.close()
img = None
RGB_Image = None
Crop True Colour Image¶
In [18]:
img = gdal.Open(S2_bands['TCI'], gdal.GA_ReadOnly)
# Define cropped area in same CRS as image (EPSG:32637)
ulx_c = 248200
lrx_c = 269800
lry_c = 4014100
uly_c = 4035700
cropped_path = os.path.join(S2proc_path, "cropped_TCI.tif")
translate_options = gdal.TranslateOptions(projWin=[ulx_c, uly_c, lrx_c, lry_c], stats=True)
translated_stack = gdal.Translate(cropped_path, img, options=translate_options)
img = None
in_crs = None
translated_stack = None
Display Cropped TCI Image using Matplotlib¶
In [19]:
img = gdal.Open(cropped_path, gdal.GA_ReadOnly)
RGB_Image_crop = img.ReadAsArray()
print(RGB_Image_crop.shape)
# Get extent and boundary
ulx, xres, xskew, uly, yskew, yres = img.GetGeoTransform()
lrx = ulx + (img.RasterXSize * xres)
lry = uly + (img.RasterYSize * yres)
print(ulx, lrx, lry, uly)
# Re-order dimensions for matplotlib.
RGB_Image_crop = numpy.dstack(((RGB_Image_crop[2]), (RGB_Image_crop[1]), (RGB_Image_crop[0]))).astype('uint8')
print(RGB_Image_crop.shape)
In [20]:
plt.figure(figsize=(10, 10))
plt.imshow(RGB_Image_crop, extent=(ulx, lrx, lry, uly))
Out[20]:
In [21]:
# Optional -- Save figure.
# plt.savefig('TestImage.png', format='png', dpi=300, bbox_inches='tight')
# Clean Up.
plt.close()
img = None
RGB_Image_crop = None
Create Empty output dataset object with GeoTiff driver¶
This creates an empty dataset object intended for a 3 band GeoTIFF stack. It uses Sentinel-2's band 2 as a template to retrieve projection and pixel information to base the empty dataset object on.
In [22]:
img = gdal.Open(S2_bands['B02'], gdal.GA_ReadOnly)
# Get raster georeference and projection info from B02.
projection = img.GetProjection()
transform = img.GetGeoTransform()
print("xOrigin: " + str(transform[0]))
print("yOrigin: " + str(transform[3]))
print("pixelWidth: " + str(transform[1]))
print("pixelHeight: " + str(transform[5]))
img_rows = img.RasterYSize
img_cols = img.RasterXSize
print("Rows: " + str(img_rows))
print("Columns: " + str(img_cols))
img = None
In [23]:
# Get GeoTIFF driver and define output filename and path.
gdal_format = 'GTiff'
driver = gdal.GetDriverByName(gdal_format)
stacked_file = "B348_stack.tif"
stacked_path = os.path.join(S2proc_path, stacked_file)
# Create output dataset object with the right size, 3 bands and as an 8-bit unsigned GeoTiff (driver).
outDs = driver.Create(stacked_path, img_cols, img_rows, 3, gdal.GDT_Byte)
# Georeference the outpute dataset object and set the projection based on B02.
outDs.SetGeoTransform(transform)
outDs.SetProjection(projection)
Out[23]:
Create false colour image in empty Dataset object¶
In [24]:
desired_bands = ['B03', 'B04', 'B08']
In [25]:
for idx, val in enumerate(desired_bands):
# Open the file.
img = gdal.Open(S2_bands[val], gdal.GA_ReadOnly)
# Read the one band as an array with the proper size.
img_band = img.GetRasterBand(1)
img_rows = img.RasterYSize
img_cols = img.RasterXSize
img_array = img_band.ReadAsArray(0, 0, img_cols, img_rows)
# Adjust outliers using numpy (very high reflectance and negative).
outData = img_array / 10000.0
outData = numpy.where((outData > 1), (1), outData)
outData = numpy.where((outData < 0), (0), outData)
img_array = None
# Convert to 8-bit.
outData = ((numpy.absolute(outData) * 255.0) + 0.5).astype(int)
# Write the array to the proper raster band in the output dataset
# object and flush to disk
outBand = outDs.GetRasterBand(idx+1)
outBand.WriteArray(outData, 0, 0)
outBand.FlushCache()
# Calculate statistics and set to band.
stats = outBand.ComputeStatistics(False)
outBand.SetStatistics(stats[0], stats[1], stats[2], stats[3])
print("Added " + val + " as band " + str(idx+1) + " to stack.")
# Clean up band specific variables.
img_band = None
band_id = None
stats = None
img = None
del outData
del outBand
outDs = None
Display False Colour (Nir, R, G) Image using Matplotlib¶
In [27]:
img = gdal.Open(stacked_path, gdal.GA_ReadOnly)
GRNir_Image = img.ReadAsArray()
print(GRNir_Image.shape)
# Get extent and boundary
ulx, xres, xskew, uly, yskew, yres = img.GetGeoTransform()
lrx = ulx + (img.RasterXSize * xres)
lry = uly + (img.RasterYSize * yres)
print(ulx, lrx, lry, uly)
# Re-order dimensions for matplotlib.
NirRG_Image = numpy.dstack(((GRNir_Image[2]), (GRNir_Image[1]), (GRNir_Image[0]))).astype('uint8')
print(NirRG_Image.shape)
In [28]:
plt.figure(figsize=(10, 10))
plt.imshow(NirRG_Image, extent=(ulx, lrx, lry, uly))
Out[28]:
In [29]:
# Optional -- Save figure.
# plt.savefig('TestImage.png', format='png', dpi=300, bbox_inches='tight')
# Clean Up.
plt.close()
img = None
GRNir_Image = None
NirRG_Image = None
Crop False Colour Image¶
In [30]:
img = gdal.Open(stacked_path, gdal.GA_ReadOnly)
# Define cropped area in same CRS as image (EPSG:32637)
ulx_c = 248200
lrx_c = 269800
lry_c = 4014100
uly_c = 4035700
cropped_path = os.path.join(S2proc_path, "false_crop.tif")
translate_options = gdal.TranslateOptions(projWin=[ulx_c, uly_c, lrx_c, lry_c], stats=True)
translated_stack = gdal.Translate(cropped_path, img, options=translate_options)
img = None
in_crs = None
translated_stack = None
In [31]:
print(gdal.Info(cropped_path))
Display Cropped False Colour (Nir, R, G) Image using Matplotlib¶
In [32]:
# Open the file.
img = gdal.Open(cropped_path, gdal.GA_ReadOnly)
GRNir_Image_crop = img.ReadAsArray()
print(GRNir_Image_crop.shape)
# Get extent and boundary
ulx, xres, xskew, uly, yskew, yres = img.GetGeoTransform()
lrx = ulx + (img.RasterXSize * xres)
lry = uly + (img.RasterYSize * yres)
print(ulx, lrx, lry, uly)
# Re-order dimensions for matplotlib.
NirRG_Image_crop = numpy.dstack(((GRNir_Image_crop[2]), (GRNir_Image_crop[1]), (GRNir_Image_crop[0]))).astype('uint8')
print(NirRG_Image_crop.shape)
In [33]:
plt.figure(figsize=(10, 10))
plt.imshow(NirRG_Image_crop, extent=(ulx, lrx, lry, uly))
Out[33]:
In [34]:
# Optional -- Save figure.
# plt.savefig('TestImage.png', format='png', dpi=300, bbox_inches='tight')
# Clean Up.
plt.close()
img = None
GRNir_Image_crop = None
NirRG_Image_crop = None
Further Tutorials and Resources¶
GDAL-related¶
Official GDAL API Tutorial for C++, C and Python
Medium: A Gentle Introduction to GDAL (Parts 1, 2, 3, 4)
Planet Lab: Getting Started with GDAL
University of Edinburgh: Simple Raster Conversion with GDAL
Digital Geography: Reading and Writing Geospatial files using GDAL and python
GeoNode: Introduction and Advanced Processing with GDAL Utilities
Documentation for all GDAL Utilities
GDAL Cheat Sheet by @dwtkns on GitHub
Python GDAL/OGR Cookbook: Raster Layers
Jupyter notebook on Exploring the GDALDataset Class by @ceholden on GitHub
Creating Natural Color RGB Composites from Landsat 8 data in Python
Basic GDAL Commands for Cartographers by Christopher Wesson
Other Python Image Libraries¶
Rasterio: essentially wraps GDAL functionality
rasterstats: allows zonal statistics from raster files or numpy arrays using geometries
scikit-image: for image processing in python
Pillow: Python Imaging Library fork for image processing
scipy.ndimage: functions for multi-dimensional image processing
Thank you for your attention!¶
E-Mail: hannah.augustin@sbg.ac.at
GitHub: @augustinh22
Website: http://hannahaugustin.at