Visualize Raster Data¶

Visualizing raster data is an effective way to explore new datasets. This example shows how to plot raster data from the Geokit dataset. Raster data is typically provided as GeoTIFF data with the .tif extension, which is used in this example as well.

In [1]:

import geokit.core.raster
import geokit.core.geom
import pathlib
from geokit.core.get_test_data import get_test_data

import geokit.core.raster import geokit.core.geom import pathlib from geokit.core.get_test_data import get_test_data

Simple Example Plot¶

The plotting function uses Matplotlib and returns a tuple containing (Axes, Handles, cbar), which allows for further modification of the plot if necessary.

In [2]:

path_to_first_aachen_raster_data = get_test_data(
    file_name="clc-aachen_clipped.tif",
)
axes_handles = geokit.core.raster.drawRaster(
    path_to_first_aachen_raster_data,
    figsize=(6, 6),
)

path_to_first_aachen_raster_data = get_test_data( file_name="clc-aachen_clipped.tif", ) axes_handles = geokit.core.raster.drawRaster( path_to_first_aachen_raster_data, figsize=(6, 6), )

Change the spatial reference system used for display.¶

To analyze the impact of using a different spatial reference system for the raster, the drawRaster function provides a convenient way to change the spatial reference system during visualization, as shown below.

In [3]:

axes_handles = geokit.core.raster.drawRaster(path_to_first_aachen_raster_data, figsize=(6, 6), srs=25832)

axes_handles = geokit.core.raster.drawRaster(path_to_first_aachen_raster_data, figsize=(6, 6), srs=25832)

Plot Multiple Rasters or Overlay Raster with Vectors¶

GeoKit provides the ability to overlay and plot multiple geometric objects next to each other. This can be useful for comparison or for highlighting individual features. It relies on Matplotlib but provides an interface that handles raster and geometric data.

To demonstrate this ability, two rasters are plotted next to each other and overlaid with vector data.

In [4]:

import matplotlib.pyplot as plt

# Create a subplot with two axes to plot the two rasters next to each other
fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(16, 12))

# Draw the first raster
axes_handle_raster_1 = geokit.core.raster.drawRaster(path_to_first_aachen_raster_data, ax=ax[0])

# Extract information about the spatial reference system
print("Spatial Reference System of the first raster as WKT string:")
print(geokit.core.raster.rasterInfo(path_to_first_aachen_raster_data).srs.ExportToWkt())
# Extract information minimum and maximum coordinates of the raster
print("xMin: ", geokit.core.raster.rasterInfo(path_to_first_aachen_raster_data).xMin)
print("yMin: ", geokit.core.raster.rasterInfo(path_to_first_aachen_raster_data).yMin)
print("xMax: ", geokit.core.raster.rasterInfo(path_to_first_aachen_raster_data).xMax)
print("yMax: ", geokit.core.raster.rasterInfo(path_to_first_aachen_raster_data).yMax)

# Create some examples points based on the dimensions of the raster
# Points are actually vector based geometries which are introduced in the next chapters

aachen_pt_1_from_tuple = geokit.core.geom.point((4012100.0 + 15000, 3031800 + 15000), srs=4258)
aachen_pt_2_from_tuple = geokit.core.geom.point((4012100.0 + 15000, 3031800 + 30000), srs=4258)
aachen_pt_3_from_tuple = geokit.core.geom.point((4012100.0 + 30000, 3031800 + 30000), srs=4258)
# Draw the examples points
axes_handle_geometries_raster_1 = geokit.core.geom.drawGeoms(
    geoms=[aachen_pt_1_from_tuple, aachen_pt_2_from_tuple, aachen_pt_3_from_tuple], ax=ax[0]
)
path_to_second_aachen_raster_data = get_test_data(file_name="aachen_eligibility.tif")

# Draw the second raster
axes_handle_raster_2 = geokit.core.raster.drawRaster(path_to_second_aachen_raster_data, ax=ax[1])

print("Spatial Reference System of the second raster as WKT string:")
print(geokit.core.raster.rasterInfo(path_to_second_aachen_raster_data).srs.ExportToWkt())

# Create lines around some features of the raster
aachen_polygon_1 = geokit.core.geom.line(
    [
        (6.18, 50.5),
        (6.23, 50.5),
        (6.4, 50.58),
        (6.4, 50.65),
        (6.25, 50.65),
        (6.19, 50.55),
        (6.18, 50.5),
    ],
    srs=4326,
)

# Draw the lines.
# Instead of passing the Axes object directly, it is also possible to pass the Axes handle of the axis in which the geometries should be drawn.
axes_handle_geometries_raster_2 = geokit.core.geom.drawGeoms(geoms=[aachen_polygon_1], ax=axes_handle_raster_2)

import matplotlib.pyplot as plt # Create a subplot with two axes to plot the two rasters next to each other fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(16, 12)) # Draw the first raster axes_handle_raster_1 = geokit.core.raster.drawRaster(path_to_first_aachen_raster_data, ax=ax[0]) # Extract information about the spatial reference system print("Spatial Reference System of the first raster as WKT string:") print(geokit.core.raster.rasterInfo(path_to_first_aachen_raster_data).srs.ExportToWkt()) # Extract information minimum and maximum coordinates of the raster print("xMin: ", geokit.core.raster.rasterInfo(path_to_first_aachen_raster_data).xMin) print("yMin: ", geokit.core.raster.rasterInfo(path_to_first_aachen_raster_data).yMin) print("xMax: ", geokit.core.raster.rasterInfo(path_to_first_aachen_raster_data).xMax) print("yMax: ", geokit.core.raster.rasterInfo(path_to_first_aachen_raster_data).yMax) # Create some examples points based on the dimensions of the raster # Points are actually vector based geometries which are introduced in the next chapters aachen_pt_1_from_tuple = geokit.core.geom.point((4012100.0 + 15000, 3031800 + 15000), srs=4258) aachen_pt_2_from_tuple = geokit.core.geom.point((4012100.0 + 15000, 3031800 + 30000), srs=4258) aachen_pt_3_from_tuple = geokit.core.geom.point((4012100.0 + 30000, 3031800 + 30000), srs=4258) # Draw the examples points axes_handle_geometries_raster_1 = geokit.core.geom.drawGeoms( geoms=[aachen_pt_1_from_tuple, aachen_pt_2_from_tuple, aachen_pt_3_from_tuple], ax=ax[0] ) path_to_second_aachen_raster_data = get_test_data(file_name="aachen_eligibility.tif") # Draw the second raster axes_handle_raster_2 = geokit.core.raster.drawRaster(path_to_second_aachen_raster_data, ax=ax[1]) print("Spatial Reference System of the second raster as WKT string:") print(geokit.core.raster.rasterInfo(path_to_second_aachen_raster_data).srs.ExportToWkt()) # Create lines around some features of the raster aachen_polygon_1 = geokit.core.geom.line( [ (6.18, 50.5), (6.23, 50.5), (6.4, 50.58), (6.4, 50.65), (6.25, 50.65), (6.19, 50.55), (6.18, 50.5), ], srs=4326, ) # Draw the lines. # Instead of passing the Axes object directly, it is also possible to pass the Axes handle of the axis in which the geometries should be drawn. axes_handle_geometries_raster_2 = geokit.core.geom.drawGeoms(geoms=[aachen_polygon_1], ax=axes_handle_raster_2)

Spatial Reference System of the first raster as WKT string:
PROJCS["ETRS89-extended / LAEA Europe",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","3035"]]
xMin:  4012100.0
yMin:  3031800.0
xMax:  4094600.0
yMax:  3111000.0
Spatial Reference System of the second raster as WKT string:
GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

Modify the Plot Appearance with the Matplotlib Interface¶

Although the drawRaster() function provides some customization options, it does not directly expose Matplotlib's full capabilities. Instead, it returns a tuple containing some of the plot's core objects.

1. The first is the axes on which the raster is drawn directly.
2. The handles of the artists drawn in the axes.
3. The color bar, if it was drawn.

Axes object¶

The Axes object provides various methods for changing the appearance of the plot. A few examples are shown below, but the full documentation can be found here. https://matplotlib.org/stable/api/axes_api.html

Colorbar¶

The appearance of the color bar can also be changed using the methods shown at: https://matplotlib.org/stable/api/colorbar_api.html#matplotlib.colorbar.Colorbar

In [5]:

# Create new plot
axes_handles = geokit.core.raster.drawRaster(path_to_first_aachen_raster_data, figsize=(6, 6))

# Modify the axes object
axes_directly = axes_handles[0]
axes_directly.set_ylabel("Latitude")
axes_directly.set_xlabel("Longitude")
axes_directly.set_title("Aachen")

# Colorbar
colorbar_directly = axes_handles[2]
colorbar_directly.minorticks_on()

# Create new plot axes_handles = geokit.core.raster.drawRaster(path_to_first_aachen_raster_data, figsize=(6, 6)) # Modify the axes object axes_directly = axes_handles[0] axes_directly.set_ylabel("Latitude") axes_directly.set_xlabel("Longitude") axes_directly.set_title("Aachen") # Colorbar colorbar_directly = axes_handles[2] colorbar_directly.minorticks_on()