The Ultimate Beginner’s Guide to Geospatial Raster Data

Most aerial photographs and imagery from satellites are raster files. This format is often used to represent real-world phenomena. If you are working with geographic data, there is a high chance you have to deal with it.

To use geographical raster files with Python, different theoretical concepts are required. Before jumping to the programmatic part, I highly recommend you follow the introductory sections.

Table of content

1. Introduction: first concepts.
2. Applications: where are rasters used?
3. Colormap: discrete and continuous colormaps to visualize rasters.
4. Georeferencing: CRS and Affine Transformations.
5. Raster’s metadata: all the data associated with rasters.
6. Rasterio: read, save, georeference and visualize raster files in Python.

---

Introduction

A raster consists of a matrix of cells (or pixels) organized into rows and columns where each cell contains a value representing information

Each pixel in a geographical raster is associated with a specific geographical location. This means that if the raster has a 1m/px resolution, __ every pixel covers an area of 1m². More details about this are given in the Georeferencing section.

Furthermore,

A raster contains one or more layers of the same size, called bands

Any type of numerical value can possibly be stored in a cell. According to the context, cells may contain integer or floating point values in different ranges.

JPGs, PNGs, and Bitmaps are all raster files but they aren’t considered in this guide as they are non-geographical. Geographical rasters are usually saved in the TIFF format.

---

Applications

As rasters can be applied in various ways, here are the most common applications:

• Satellite images
• Thematic maps
• Digital elevation models (DEM)

Satellite images

Images from satellites are usually saved on multiband rasters. The electromagnetic spectrum is divided into multiple portions which can be sensed by a satellite. Not all of them belong to the visible spectrum, but often some will be in the infrared and invisible to the human eye.

A Raster file perfectly suits this type of imagery because each electromagnetic spectrum portion sensed by the satellite can be stored in a band.

Sentinel-2, which is one of the most popular satellites, takes photos using 13 spectral bands, one part from the visible spectrum and the other from the infrared. As a result, each output file is a raster with 13 bands (3 of them are Red, Green and Blue).

The following is a photo taken from Sentinel-2 (only the RGB bands are shown):

Thematic map

A thematic map is used to classify a geographical area. Each zone is associated with a particular class sharing some characteristics. For example, we can classify an agricultural area according to the type of plantations. Rasters are perfect for this task because each cell can store the integer value representing the class to which the pixel associated area belongs.

Below is an example of a thematic map from Lombardia, an Italian region. Each pixel stores a value between 0 and 6 depending on the class:

Digital elevation model (DEM)

A digital elevation model is used to represent surface reliefs. DEM is a type of raster whose pixels contain float values: the elevation values of the surface.

What is here represented is a DEM of a Mars surface area:

To allow the visualization of the first file only the visible bands were showed, but the second and third files were not directly visualizable because the value saved in each cell is not a color but a piece of information.

In the next section, the workaround required to visualize this kind of raster files will be focussed on.

---

Colormap

As rasters have no constraints on the type and range of numerical values you can store, it is not always possible to show them visually. For example, the last images I showed above are two single-band rasters: in the former, each pixel is an integer between 0 and 6 while in the latter, a float between -4611 and –3871. Their information doesn’t represent a color.

To visualise these kinds of rasters we use a colormap, i.e.

A function which maps cell values to colors

Thus, when visualizing a raster through a colormap, its values will be replaced by colors.

There are 2 main types of colormaps: continuous and non-continuous.

Non-continuous colormaps

They are made by defining a piecewise function using value-color pairs.

In the thematic map example, I defined 7 value-color pairs: <0, black>, <1, red>, <3, orange>, <4, yellow>, <5, blue>, <6, gray>, <2, green>.

This method is often used when the raster includes a small set of values.

Continuous colormaps

They are made by associating the interval of the raster values with an interval of colors using a continuous function.

Usually, before applying this type of colormap, all the raster values are scaled in the range [0, 1] using the following formula:

In the grayscale colormap, values in the [0,1] range are associated with gray values in the [0,255] range using a linear function:

A colorbar shows the result in a proper way:

In this case the 0 value is represented by the black color, the 1 value by white and all the values between them by different shapes of gray.

To better visualize a raster, we can also define an RGB colormap, in which each raster value will be associated with a red, green and blue value.

The following is a popular colormap in literature known as Turbo:

This colormap starts with blue shades for the lowest values and ends with red shades for the highest:

In the DEM example, I used this type of colormap to convert elevation information into colors. This is the kind of colormap used when mapping a range of values into colors (instead of a small set like in non-continuous colormaps).

---

Georeferencing

Each cell in a geographical raster covers a particular geographical area and its coordinates, represented by the row and the column, can be converted into real-world geographic coordinates. The translation process uses two components: the Coordinate Reference System (CRS) and the Affine Transformations.

Before going forward, it’s important to know that the earth’s shape is approximated by means of a geometrical figure, known as the ellipsoid of revolution or spheroid. As this figure is an approximation, multiple spheroids have been defined over the years using axes with different sizes.

Coordinate Reference System (CRS)

CRS is a framework used to precisely measure locations on the surface of the Earth as coordinates

Each CRS is based on a specific spheroid, thus, if two CRS use different spheroids, the same coordinates refer to two different locations.

CRS can be divided into:

• Geographic coordinate systems, which utilize angular units (degrees). The angular distances are measured from defined origins.
• Projected coordinate systems, based on a geographic coordinate system. It uses spatial projection, which is a set of mathematical calculations performed to flatten the 3D data onto a 2D plane, to project the spheroid. It utilizes linear units (feet, meters, etc.) to measure the distance (on both axes) of the location from the origin of the plane.

One of the most popular geographic coordinate systems is WGS84, ** also known as EPSG:4326. It’s based on a spheroid with a semi-minor axis (known as equatorial radius) equal to 6378137 m and a semi-major axis (known as polar radius) equal to 6356752 m. WGS84 uses Latitude to find out how far north or south a place is from the Equator and Longitud**e to find out how far east or west a place is from the Prime meridian. Both are specified in degrees.

e.g. 40° 43′ 50.1960” N, 73° 56′ 6.8712” W is the New York City position using latitude and longitude.

While one of the most popular projected coordinate systems is UTM/WGS84, also known as EPSG:32632. It projects the WGS84 spheroid into a plane and then coordinates are defined using (x, y) in meters.

The CRS used in the raster file depends on multiple factors such as when the data was collected, the geographic extent of the data and the purpose of the data. Keep in mind that you can convert the coordinates of a CRS to another. There are also CRS used to georeference surfaces outside our earth, such as Moon and Mars.

Affine transformations

Georeferenced rasters use affine transformations to map from image coordinates to real-world coordinates (in the format defined by the CRS).

Affine transformations are used to map pixel positions into the chosen CRS coordinates

An affine transformation is any transformation that preserves collinearity (three or more points are said to be collinear if they all lie on the same straight line) and the ratios of distances between points on a line.

These are all the affine transformations:

In georeferencing, most of the time, only scaling and translation transformations are required. Applying them with the right coefficients is what allows translating raster cell coordinates into real-world coordinates. When reading a geographical raster, these coefficients are already defined inside the metadata.

The relationship used for the conversion is:

If scale_x and scale_y are the pixel_width and pixel_height in CRS unit (degrees, meters, feet, etc.), r is the rotation of the image in real-world, x_origin and y_origin are the coordinates of the top left pixel of the raster in CRS unit, the parameters are:

• A = scale_x ∙ cos(r)
• B = scale_y ∙ sin(r)
• C = x_origin ∙ cos(r) + y_origin ∙ sin(r)
• D = scale_x ∙ sin(r)
• E = scale_y ∙ cos(r)
• F = x_origin ∙ sin(r) + y_origin ∙ cos(r)

Keep in mind that one or more of scale_x, scale_y, x_origin and y_origin can be negative depending on the CRS used.

As most images are north-up, and thus r = 0, parameters can be simplified:

• A = scale_x
• B = 0
• C = x_origin
• E = scale_y
• D = 0
• F = y_origin

A and E define the scaling ratio while C and F the translation from the origin.

---

Metadata

Each geographical raster has metadata associated. Here are the most important fields:

CRS

This field stores the information of the Coordinate Reference System, such as the name, unit of measurement, spheroid axes and the coordinates of the origin.

Tranformation

It stores the coefficients A, B, C, D, E, F, used to map raster pixel coordinates to CRS coordinates.

Data type

It is usually known as dtype. It defines the type of data stored in the raster such as Float32, Float64, Int32, etc.

NoData Value

Each cell of a raster must hold a value and rasters don’t support null values. If for some cells the source that generated the raster couldn’t provide a value, they are filled using the nodata value. If the nodata value is set to 0, this means that when 0 is read it is not a value to consider because it indicates that the source couldn’t provide a correct value. Usually, rasters with Float32 data type, set the nodata value to ≈ -3.4028235 ∙ 10³⁸.

Width, Height and Band count

They are respectively the width of each band, the height of each band and the number of bands of the raster.

Driver

A driver provides more features to raster files. Most of the time, geographical rasters use the GTiff driver, which allows georeferencing information to be integrated into the file.

---

Rasterio

The first library ever made for accessing geographical raster files is [Geospatial Data Abstration Library, GDAL](http://GT(1) × GT(5) = 10m × 10m). It was first developed in C and then extended to Python. This way, the Python version provides only little abstraction from the C version. Rasterio which is based on GDAL, tries to solve this problem by providing an easier and higher-level interface.

To install the latest version of Rasterio from PyPI use:

pip install rasterio

These are the required imports:

import rasterio
from rasterio.crs import CRS
from rasterio.enums import Resampling
from matplotlib import pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np
from pprint import pprint

Reading a raster

To open a raster file use:

raster = rasterio.open('raster.tiff')

To print the number of bands use:

print(raster.count)

To read all the raster as a NumPy array use:

raster_array = raster.read()  # shape = (n_bands x H x W)

Alternatively, to read only a specific band use:

first_band = raster.read(1)  # shape = (H x W)

Keep in mind that the band index begins from 1.

To read all the metadata associated with a raster use:

metadata = dataset.meta
pprint(metadata)

Output:

{'count': 1,
 'crs': CRS.from_epsg(32632),
 'driver': 'GTiff',
 'dtype': 'uint8',
 'height': 2496,
 'nodata': 255.0,
 'transform': Affine(10.0, 0.0, 604410.0,
       0.0, -10.0, 5016150.0),
 'width': 3072}

When reading a raster, there might be nodata values, and it is advisable to replace them with NaN values. To do this use:

first_band[first_band == metadata['nodata']] = np.nan

To resize a raster by a given factor, define first the output shape:

out_shape = (raster.count, int(raster.height * 1.5), int(raster.width * 1.5))

Then use:

scaled_raster = raster.read(out_shape=out_shape,
                resampling=Resampling.bilinear)

Keep in mind that after scaling a raster, A and F coefficients of the affine transformation must be changed to the new pixel resolution, otherwise, georeferencing will give the wrong coordinates.

Visualization

To show a raster band with a colormap and a colorbar use:

fig, ax = plt.subplots()
im = ax.imshow(raster.read(1), cmap='viridis')
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.10)
fig.colorbar(im, cax=cax, orientation='vertical')
plt.savefig('cmap_viz')

Output:

In the case of multiple-band satellite rasters, to show the first 3 visible bands use:

rgb_bands = raster.read()[:3]  # shape = (3, H, W)
plt.imshow(rgb_bands)

Usually, Satellite Imagery stores RGB in the first three bands. If this is not the case, the index has to be changed according to the order of the bands.

Georeferencing

Coordinates from the following methods will be returned and must be provided using the current CRS unit.

To find the real coordinates of a pixel (i, j), where i is the row and j the column, use:

x, y = raster.xy(i, j)

To do the opposite use:

i, j = raster.index(x, y)

To show the bounds of the data:

print(raster.bounds)

Output:

BoundingBox(left=604410.0, bottom=4991190.0, right=635130.0, top=5016150.0)

The raster in this example was georeferenced using ESPG:32632 as CRS, therefore the output coordinates are in meters.

Save a raster

Step 1 – Find the EPSG code of the desired CRS, then retrieve its information:

crs = CRS.from_epsg(32632)

In this example, EPSG:32632 is used, which was mentioned in the fourth section of this guide.

Step 2— Define an Affine transformation:

transformation = Affine(10.0, 0.0, 604410.0, 0.0, -10.0, 5016150.0)

The purpose of the coefficients A, B, C, D, E, F was explained in the fourth section of this guide.

Step 3 – Save the raster:

# a NumPy array representing a 13-band raster
array = np.random.rand(13,3000,2000)
with rasterio.open(
    'output.tiff',
    'w',
    driver='GTiff',
    count=array.shape[0],  # number of bands
    height=array.shape[1],
    width=array.shape[2],
    dtype=array.dtype,
    crs=crs,
    transform=transform
) as dst:
    dst.write(array)

In case georeferencing is not required set crs = Noneand transform = None.

There is another syntax of the write method:

dst.write(array_1, 1)
# ...
dst.write(array_13, 13)

But in most cases, it’s easier to write a single 3d array.

---

Conclusion

This guide has shown how in a raster, a set of one or more same-size matrices called bands, each cell contains a piece of information. This information changes according to the tasks, such as satellite imagery, thematic maps, and digital elevation model. Furthermore, depending on the application you might need colormaps to visualize it. Also, you found out that using a coordinate reference system and affine transformations, it’s possible to map each cell position to real-world coordinates. In the end, you saw that Rasterio makes easy to perform read, write, visualization and georeference operations. In case you need a program to open rasters, QGIS and ArcGIS are good options.

Additional resources

• Coordinate reference systems
• Map projections
• Affine transformations library
• Rasterio official documentation

Thanks for reading, I hope you have found this useful.