The world’s leading publication for data science, AI, and ML professionals.

Reading and Visualizing GeoTiff | Satellite Images with Python

An introduction to python libraries for working with GeoTiff or satellite images.

Unsplash Image Link
Unsplash Image Link

Reading images is super easy in Python, as there are many libraries available for reading, editing, and visualizing different formats of images. Examples of these libraries are Matplotlib, OpenCV, Pillow, etc. These libraries work perfectly well with popular image formats such as .png, .tiff, .jpg, .jpeg, etc. but they don’t seem to work for GeoTiff images.

The reason can be easily guessed if we understand how GeoTiff format is different from other image formats. A simple definition would be something given below.

A GeoTIFF is a standard .tif or image file format that includes additional spatial (georeferencing) information embedded in the .tif file as tags.

By tags, it is meant that it has some additional metadata like Spatial extent, CRS, Resolution, etc. along with the pixel values. It is a popular distribution format for satellite and aerial photography imagery.

This article discusses different ways of reading and visualizing these images with python using a jupyter notebook. The libraries used are GDAL, rasterio, georaster, and Matplotlib(for visualization). These libraries will help us in converting those images to simple numpy array format from there we can also perform other image transformations using numpy or TensorFlow although I have not discussed transformations in this article.

Installation steps are mentioned in the documentation of these libraries, I have provided their links in the references section, please go through them.

Using RasterIO

RasterIO is implemented by Mapbox and it provides python API for reading geo-referenced dataset. Also, the documentation mentions that unlike GDAL python bindings, it doesn’t suffer from dangling C pointer and other pointer issues that can crash your program.

import rasterio
from rasterio.plot import show
fp = r'GeoTiff_Image.tif'
img = rasterio.open(fp)
show(img)

Don’t confuse yourself with the x and y-axis scale values, they are just longitude and latitude values. If you want to read individual bands use the below code.

fp = r'GeoTiff_Image.tif'
img = rasterio.open(fp) 
# mention band no. in read() method starting from 1 not 0
show(img.read(3))

Image related data with RasterIO

# No. of Bands and Image resolution
dataset.count, dataset.height, dataset.width
# Coordinate Reference System
dataset.crs
image related data with RasterIO
image related data with RasterIO

Using georaster

The version installed with pip for georaster is 1.25. This version may give you some issues, if you face any, just go to your georaster installation directory first. The code below will tell you the installation directory path.

import georaster
print(georaster.__file__)

Now, find the georaster.py file there and change its content with content in this georaster.py file.

Below is the code for visualizing the image having more than one band (or channels).

import georaster
import matplotlib.pyplot as plt
# Use SingleBandRaster() if image has only one band
img = georaster.MultiBandRaster('GeoTiff_Image.tif')
# img.r gives the raster in [height, width, band] format 
# band no. starts from 0
plt.imshow(img.r[:,:,2])

Plotting image without matplotlib

We can directly use the plot method instead of selecting and visualizing individual bands.

georaster plot method
georaster plot method

Using GDAL

This is the most popular library when dealing with GeoTiff images but it is sometimes difficult to install and takes lots of effort to be finally able to use it. GDAL has most of its methods and classes implemented in C++, we are using its python bindings here. Most of the libraries like georaster utilize GDAL and provides a nice and simple python interface to it.

from osgeo import gdal
import matplotlib.pyplot as plt
dataset = gdal.Open('GeoTiff_Image.tif', gdal.GA_ReadOnly) 
# Note GetRasterBand() takes band no. starting from 1 not 0
band = dataset.GetRasterBand(1)
arr = band.ReadAsArray()
plt.imshow(arr)

Image related data with GDAL

# For no. of bands and resolution
img.RasterCount, img.RasterXSize, img.RasterYSize
# stats about image
img.GetStatistics( True, True )
Jupyter Notebook implementation for image stats
Jupyter Notebook implementation for image stats

References

Installation

  1. https://mothergeo-py.readthedocs.io/en/latest/development/how-to/gdal-ubuntu-pkg.html
  2. https://github.com/GeoUtils/georaster
  3. https://rasterio.readthedocs.io/en/latest/installation.html
  4. https://pypi.org/project/GDAL/

Docs

  1. https://georaster.readthedocs.io/en/latest/api.html#multibandraster
  2. https://rasterio.readthedocs.io/en/latest/quickstart.html
  3. https://gdal.org/python/index.html

Related Articles