48 Georeferencing Rasters using Rasterio in GemGIS

Rasters with no spatial reference can be georeferenced so that they are located properly in space. This can be done using rasterio.

Source: Powell, D. (1995): Interpretation geologischer Strukturen durch Karten - Eine praktische Anleitung mit Aufgaben und Lösungen, page 15, figure 10 A, Springer Verlag Berlin, Heidelberg, New York, ISBN: 978-3-540-58607-4.

Set File Paths and download Tutorial Data

If you downloaded the latest GemGIS version from the Github repository, append the path so that the package can be imported successfully. Otherwise, it is recommended to install GemGIS via pip install gemgis and import GemGIS using import gemgis as gg. In addition, the file path to the folder where the data is being stored is set. The tutorial data is downloaded using Pooch (https://www.fatiando.org/pooch/latest/index.html) and stored in the specified folder. Use pip install pooch if Pooch is not installed on your system yet.

import gemgis as gg

file_path ='data/48_georeferencing_rasters_using_rasterio/'

WARNING (theano.configdefaults): g++ not available, if using conda: `conda install m2w64-toolchain`
C:\Users\ale93371\Anaconda3\envs\test_gempy\lib\site-packages\theano\configdefaults.py:560: UserWarning: DeprecationWarning: there is no c++ compiler.This is deprecated and with Theano 0.11 a c++ compiler will be mandatory
  warnings.warn("DeprecationWarning: there is no c++ compiler."
WARNING (theano.configdefaults): g++ not detected ! Theano will be unable to execute optimized C-implementations (for both CPU and GPU) and will default to Python implementations. Performance will be severely degraded. To remove this warning, set Theano flags cxx to an empty string.
WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

gg.download_gemgis_data.download_tutorial_data(filename="48_georeferencing_rasters_using_rasterio.zip", dirpath=file_path)

Downloading file '48_georeferencing_rasters_using_rasterio.zip' from 'https://rwth-aachen.sciebo.de/s/AfXRsZywYDbUF34/download?path=%2F48_georeferencing_rasters_using_rasterio.zip' to 'C:\Users\ale93371\Documents\gemgis\docs\getting_started\tutorial\data\48_georeferencing_rasters_using_rasterio'.

Loading Data

In order to demonstrate the capabilities of rasterio in georeferencing rasters, we are using the the raster used in the first tutorials.

import rasterio
import matplotlib.pyplot as plt

raster = rasterio.open(file_path + 'task1.tif')
raster

WARNING (theano.configdefaults): g++ not available, if using conda: `conda install m2w64-toolchain`
C:\Users\ale93371\Anaconda3\envs\test_gempy\lib\site-packages\theano\configdefaults.py:560: UserWarning: DeprecationWarning: there is no c++ compiler.This is deprecated and with Theano 0.11 a c++ compiler will be mandatory
  warnings.warn("DeprecationWarning: there is no c++ compiler."
WARNING (theano.configdefaults): g++ not detected ! Theano will be unable to execute optimized C-implementations (for both CPU and GPU) and will default to Python implementations. Performance will be severely degraded. To remove this warning, set Theano flags cxx to an empty string.
WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

<open DatasetReader name='../../../../gemgis_data/data/48_georeferencing_rasters_using_rasterio/task1.tif' mode='r'>

raster.read(1)

array([[217, 229, 227, ..., 227, 233, 204],
       [245, 254, 254, ..., 254, 254, 246],
       [244, 254, 254, ..., 254, 254, 240],
       ...,
       [244, 254, 254, ..., 254, 254, 230],
       [245, 254, 254, ..., 254, 254, 246],
       [244, 254, 254, ..., 254, 254, 246]], dtype=uint8)

The origin of the raster is located at 0,0.

from rasterio.plot import show

show(raster)

<AxesSubplot:>

Defining GroundControlPoints

Ground control points are used by Rasterio to map a row and column of an image to an x and y (and z) value. We now set the origin to 100,100 and add additional 100 meters to the margins of the raster.

raster.read(1).shape

(1069, 972)

point1 = rasterio.control.GroundControlPoint(row=0, col=0, x=100, y=1169)
point1

GroundControlPoint(row=0, col=0, x=100, y=1169, id='8281849d-0eb4-464b-b99e-536d3d4edd56')

point2 = rasterio.control.GroundControlPoint(row=0, col=972, x=1072, y=1169)
point2

GroundControlPoint(row=0, col=972, x=1072, y=1169, id='6ec1385c-2d33-4614-adc4-07f024878aff')

point3 = rasterio.control.GroundControlPoint(row=1069, col=0, x=100, y=100)
point3

GroundControlPoint(row=1069, col=0, x=100, y=100, id='e3dba004-87c6-4f3b-b30e-4e953d49400c')

point4 = rasterio.control.GroundControlPoint(row=1069, col=972, x=1072, y=100)
point4

GroundControlPoint(row=1069, col=972, x=1072, y=100, id='abd9e2e4-782c-421e-b794-614ab58c74b0')

A list of ground control points is created

points = [point1, point2, point3, point4]
points

[GroundControlPoint(row=0, col=0, x=100, y=1169, id='8281849d-0eb4-464b-b99e-536d3d4edd56'),
 GroundControlPoint(row=0, col=972, x=1072, y=1169, id='6ec1385c-2d33-4614-adc4-07f024878aff'),
 GroundControlPoint(row=1069, col=0, x=100, y=100, id='e3dba004-87c6-4f3b-b30e-4e953d49400c'),
 GroundControlPoint(row=1069, col=972, x=1072, y=100, id='abd9e2e4-782c-421e-b794-614ab58c74b0')]

Creating Affine Transformation

An affine transformation is created. The data array can then be saved in the next step.

transformation = rasterio.transform.from_gcps(points)
transformation

Affine(1.0, 0.0, 100.0,
       0.0, -1.0000000000000002, 1169.0000000000005)

Saving the georeferenced raster

The raster is now saved as new tif file with the dimensions of the data array, the projection, the transformation and the actual data.

with rasterio.open(
    file_path + 'task1_new.tif',
    'w',
    driver='GTiff',
    height=raster.read(1).shape[0],
    width=raster.read(1).shape[1],
    count=1,
    dtype=raster.read(1).dtype,
    crs='+proj=latlong',
    transform=transformation,
) as dst:
    dst.write(raster.read(1), 1)

Opening new raster and plotting

The new raster can now be opened as usual.

raster_new = rasterio.open(file_path + 'task1_new.tif')
raster_new

<open DatasetReader name='../../../../gemgis_data/data/48_georeferencing_rasters_using_rasterio/task1_new.tif' mode='r'>

When plotting the raster, it can be seen that the origin is now correctly placed at 100,100 instead of 0,0 as for the original raster. The georeferencing was successful.

from rasterio.plot import show

show(raster_new)

<AxesSubplot:>