The Power of GDAL

The role of a remote sensing scientist is primarily to capture and extract useful information from spatial data.

Taking spatial data and extracting useful, actionable information is an on-going challenge. There are a wide range of commercial software packages to assist in this task that are dedicated to analysing spatial data. However, out-of-box solutions for specific tasks do not always exist - especially when limited resources place specialised commercial software out of reach.

In these circumstances, an analyst is left to develop and implement their own solution. But what options are there for analysts to achieve these goals in both a low-cost and efficient manner? One option is the Geospatial Data Abstraction Layer (GDAL). 

GDAL is a freely available, open source library for accessing, analysing, and managing spatial data. At its core, it provides abstraction layers for both raster and vector data that can support a wide range of file formats across different storage systems. Analysts can access GDAL at two different levels: through a high level suite of command line driven utilities or through a lower level programming API. 

Command line functions

Command line functions cover a range of common requirements that can be executed using a single command. These commands include the translation of data between different file formats and structures, the reprojection of data to different spatial reference systems and the conversion between vector and raster data (Table 1).  When organised into batch scripts, GDAL command line functions are capable of performing repeated processing on large volumes of data - however, they are naturally limited to a select number of common geospatial processing tasks. The true flexibility of GDAL is revealed when utilising its API.
 

Table 1: A list of selected, high use GDAL command line functions. For an exhaustive overview of GDAL command line functions please visit https://gdal.org/programs/.

GDAL programming API

The GDAL API provides libraries for direct access and rapid processing of both raster and vector spatial data. The API itself is available across several programming languages, primarily for Python and C++, but also more recently for C# and Java.

The Python programming language has, over the past decade, become the dominant language for computer vision and machine learning. Through the Python GDAL API, a broad range of open source machine learning techniques are available for geospatial analysis. 
 

GDAL in action: Image segmentation through unsupervised image classification

Image segmentation is a technique by which we group neighbouring, similar pixels within raster data together into objects. 

Here we illustrate a relatively simple approach to image segmentation using an unsupervised image classification – a k-means cluster analysis. This approach is performed in Python, and utilises a range of libraries in addition to GDAL. We will be performing image segmentation of an image of Tasmania captured using the Sentinel2 satellite (see Figure 1).

Figure 1: Tasmania as captured by the Sentinel 2 multispectral satellite.

1. We begin by importing the libraries. In addition to the standard GDAL libraries, we will also need numpy to handle the data arrays, sklearn for its capacity to construct classification models, and scipy for additional data filtering.
   
   import ogr, os, osr, gdal
   import numpy as np
   from sklearn.cluster import KMeans
   from scipy import ndimage
    
2. The next step is to import our data and extract the required image band information. We will work with the data as a flattened array where each element represents a single pixel with a red, green and blue value. We will later transform the result back into a 2D image.
   
   ds_in = gdal.Open( “C:/Spatial/S2_Tasmania.tif” )
   im_blue = ds_in.GetRasterBand(1).ReadAsArray().flatten()
   im_green = ds_in.GetRasterBand(2).ReadAsArray().flatten()
   im_red = ds_in.GetRasterBand(3).ReadAsArray().flatten()
   im_stacked = np.stack( [im_blue, im_green, im_red] )
   im_stacked = im_stacked.T
    
3. There is a large amount of data in this image, which would result in a very time-consuming k-means model building phase. Instead, we will slice this data to produce a training subset. Here we select every 1,000 pixel.
   
   im_stacked_subset = im_stacked[::1000,:]
    
4. We can then use this subset to generate a k-means model. Here we will attempt to identify eight clusters within the data.
   
   model = KMeans(n_clusters=8).fit(im_stacked_subset)
    
5. We will then use this model to predict the class of the original data.
   
   prediction = model.predict(im_stacked)
    
6. We can now transform this prediction layer back into a 2D image using the dimensions of the original raster data source.
   
   x_dim, y_dim = ds_in.RasterXSize, ds_in.RasterYSize
   prediction = prediction.reshape(y_dim, x_dim)
    
7. To clean up the image, we can apply a median filter. This will pass a moving kernel (a square radius of 10 pixels) across the image. Each pixel will be reassigned the median value within this moving kernel. This is quite a strong filter, but it helps to generalise the data structure by removing local fluctuations (see Figure 2).
   
   prediction = ndimage.median_filter(prediction, size=10)

   

   Figure 2: Illustration of the effects of a median filter. Areas of local fluctuation are smoothed out in favour of greater spatial generalisation.

8. We are now ready to output our image as a new raster image. First though, we will need a spatial reference system (SRS) and geotransformation information (including pixel size). We can extract these from our original dataset.
   
   prj = ds_in.GetProjectionRef()
   srs = osr.SpatialReference()
   srs.ImportFromWkt( prj )
   g_transform = ds_in.GetGeoTransform()
    
9. We can use these to create a new output raster file.
   
   driver = gdal.GetDriverByName('GTiff')
   ds_out = driver.Create(“C:/Spatial/S2_Tasmania_Prediction.tif”, x_dim, y_dim, 1, gdal.GDT_Byte)
   ds_out.SetGeoTransform( g_transform )
   ds_out.SetProjection( srs.ExportToWkt() )
    
10. We can now write our prediction layer to the output raster file and close it.
    
    ds_out.GetRasterBand(1).WriteArray( prediction )
    del ds_out
     
11. We want to convert our raster file to a vector file so that we can explore our image segmentation. First, we reopen our prediction raster file and select the first image band.
    
    r_ds = gdal.Open(fnr_out)
    r_bnd = r_ds.GetRasterBand(1)
     
12. We need to generate a shapefile to store our vector data in. Once again, we will need a spatial reference system, but we can use the one extracted from the original raster file. After that, we’ll need to create a new layer within the shapefile to store the vectors.
    
    driver = ogr.GetDriverByName("ESRI Shapefile")
    v_ds = driver.CreateDataSource(“C:/Spatial/S2_Tasmania_Prediction_Vectors.shp”)
    srs = osr.SpatialReference()
    srs.ImportFromWkt(r_ds.GetProjection())
    v_lay = v_ds.CreateLayer('Vector', srs = srs )
     
13. We can then use GDAL to convert the raster prediction layer into vectors, and store them in our new shapefile. Finally, we need to close all our raster and vector datasets.
    
    gdal.Polygonize( r_bnd, r_bnd, v_lay, -1, [], callback=None )
    del r_ds
    v_ds.Destroy()
    
    The result is a simple image segmentation of relatively homogeneous regions across Tasmania (see Figure 3). We can use this in further analysis, extracting information to compare and contrast different land cover types.

     

    Figure 3: Illustration of Tasmania Sentinel 2 data following a K-means clustering based image segmentation.

A flexible suite of geospatial command line functionality 

GDAL offers an open source solution for data processing. It is a flexible suite of geospatial command line simply coupled with an API that is capable of handling a wide range of data formats and structures. 

It’s a strong option when performing repetitive tasks over a large volume of spatial data, or when adapting workflows to fit within limited budgets or resource constraints. When coupled with a minimal amount of programming experience, a remote sensing analyst is able to successfully develop and implement a wide range of novel analyses when adopting the GDAL API.

Back To News Stories