Advanced GLAES workflow examples¶

This notebook introduces several practical GLAES features that build on the basic workflow. It is intended for users who already understand how to create an ExclusionCalculator and want to explore additional exclusion methods, visualization options, and raster handling.

Overview:¶

1. How to inspect region grid properties
2. How to apply exclusion rules from a table using exclusion_set
3. How to restart from a precomputed result
4. How to visualize results with basemaps
5. Why padExtent matters near borders
6. How to use excludeRasterType with value ranges
7. How to use inclusion mode and inversion
8. How to work with DEM mosaics for elevation and slope exclusions
9. Working with excludeVectorType()
10. Visualization methods of ec.draw

In [1]:

import geokit as gk
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import glaes as gl

import geokit as gk import matplotlib.pyplot as plt import numpy as np import pandas as pd import glaes as gl

In [2]:

# Use built-in test data so the examples run without extra downloads
region_path = gl._test_data_["aachenShapefile.shp"]
clc_path = gl._test_data_["clc-aachen_clipped.tif"]

# Create a standard exclusion calculator used in several examples
ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100)

# Use built-in test data so the examples run without extra downloads region_path = gl._test_data_["aachenShapefile.shp"] clc_path = gl._test_data_["clc-aachen_clipped.tif"] # Create a standard exclusion calculator used in several examples ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100)

1. Inspect region grid information¶

Each ExclusionCalculator works on a raster grid defined by the region, spatial reference system, and pixel resolution.

Understanding this grid is useful because it helps you see:

• the spatial extent being processed
• the raster mask used internally
• the size of the calculation grid

In [3]:

# Create a fresh ExclusionCalculator so the grid properties are easy to inspect.
ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100)
# The bounding box around the region, expressed in the chosen coordinate system.
ec.region.extent

# Create a fresh ExclusionCalculator so the grid properties are easy to inspect. ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100) # The bounding box around the region, expressed in the chosen coordinate system. ec.region.extent

Out[3]:

xMin: 4037300.000000
xMax: 4067700.000000
yMin: 3049300.000000
yMax: 3100200.000000
srs: 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"]]

In [4]:

# Region mask of the raster grid covering the region's bounding box.
# True marks pixels inside the study region; False marks pixels outside it.
ec.region.mask

# Region mask of the raster grid covering the region's bounding box. # True marks pixels inside the study region; False marks pixels outside it. ec.region.mask

Out[4]:

array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]], shape=(509, 304))

In [5]:

# The shape of the internal grid as (rows, columns).
ec.region.mask.shape

# The shape of the internal grid as (rows, columns). ec.region.mask.shape

Out[5]:

(509, 304)

2. Applying an Exclusion Set¶

Multiple exclusion rules can be stored in a CSV file and loaded as a pandas.DataFrame. Each row represents one exclusion rule.

In [6]:

exclusion_set = pd.read_csv(gl._test_data_["sample_exclusion_set.csv"])
exclusion_set

exclusion_set = pd.read_csv(gl._test_data_["sample_exclusion_set.csv"]) exclusion_set

Out[6]:

	name	type	value	buffer
0	clc	raster	[12-22]	NaN
1	clc	raster	[1-2]	1000.0
2	osm_roads	vector	type='motorway' OR type='primary' OR type='trunk'	200.0

The name column identifies the dataset used by a rule, while type and value define how the exclusion is applied. Optional columns such as buffer, exclusion_mode, and invert can further control each rule.

The datasets referenced in name are passed to excludeSet as keyword arguments. For example, clc and osm_road in the table correspond to the clc=... and osm_roads=... arguments below.

In [7]:

# Start with a new ExclusionCalculator for this example.
ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100)

# Apply all exclusion rules listed in the table.
# The additional keyword arguments provide the datasets referenced inside the table.
ec.excludeSet(
    exclusion_set=exclusion_set,
    clc=gl._test_data_["clc-aachen_clipped.tif"],
    osm_roads=gl._test_data_["aachenRoads.shp"],
)

ec.draw()

# Start with a new ExclusionCalculator for this example. ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100) # Apply all exclusion rules listed in the table. # The additional keyword arguments provide the datasets referenced inside the table. ec.excludeSet( exclusion_set=exclusion_set, clc=gl._test_data_["clc-aachen_clipped.tif"], osm_roads=gl._test_data_["aachenRoads.shp"], ) ec.draw()

Excluding Raster clc with value [12-22], buffer None, mode exclude, and invert False

Excluding Raster clc with value [1-2], buffer 1000.0, mode exclude, and invert False

Memory usage during calc: 228.87109375 MB
Memory usage during calc: 229.47265625 MB

Excluding Vector osm_roads with where-statement "type='motorway' OR type='primary' OR type='trunk'", buffer 200.0, mode exclude, and invert False 

Done!

Memory usage during calc: 252.9453125 MB

Out[7]:

<Axes: >

3. Reinitialize an ExclusionCalculator from a pre-computed result¶

If a previous exclusion result has already been saved as a raster, you can use it as the starting point for a new ExclusionCalculator. Note: Reusing a previously saved .tif as initialValue reproduces the previous availability result exactly only if the new ExclusionCalculator uses the same spatial context as the one used to create the file. This includes the same region, SRS, pixel resolution, bounds/extent, and raster alignment. If these settings differ, GLAES will warp the raster to the new calculator context, which provides a similar starting layer but may not be pixel-identical to the original result.

In [8]:

## Set the initial value to those we computed previously. This is particularly useful if the original calculation took a while
ec = gl.ExclusionCalculator(
    region_path,
    srs=3035,
    pixelRes=100,
    initialValue="data/aachens_best_pv_spots.tif",
)

ec.draw()

## Set the initial value to those we computed previously. This is particularly useful if the original calculation took a while ec = gl.ExclusionCalculator( region_path, srs=3035, pixelRes=100, initialValue="data/aachens_best_pv_spots.tif", ) ec.draw()

Memory usage during calc: 261.4296875 MB

Out[8]:

<Axes: >

4. Draw exclusions with a basemap¶

drawWithSmopyBasemap() can be used to display exclusion results on top of a map background. This is useful when the result should be interpreted in its geographic context, for example to compare excluded areas with roads, settlements, or other recognizable landscape features. Keep the zoom level moderate to avoid slow loading.

In [9]:

ec.drawWithSmopyBasemap?

ec.drawWithSmopyBasemap?

In [10]:

ax = ec.drawWithSmopyBasemap(zoom=11, figsize=(40, 20))

ax = ec.drawWithSmopyBasemap(zoom=11, figsize=(40, 20))

5. Considering region boundary effects with padExtent¶

Features outside the study region can still affect the available area inside it. For example, a road just outside the region boundary may have a buffer that extends into the region.

Without additional padding, features outside the region may not be considered during processing. As a result, exclusions close to the boundary can be incomplete.

padExtent extends the area considered during processing beyond the region boundary. The value specifies the additional distance in the units of the calculator's SRS. The final availability, however, is still evaluated only within the original study region.

Below, we compare the processing grid with and without padExtent.

In [11]:

# No padding: the grid covers only the region bounding box.
ec0 = gl.ExclusionCalculator(
    region=region_path,
    srs=3035,
    pixelRes=100,
    padExtent=0,
)

print("Extent (padExtent=0):", ec0.region.extent)
print("Mask shape (padExtent=0):", ec0.region.mask.shape)

# No padding: the grid covers only the region bounding box. ec0 = gl.ExclusionCalculator( region=region_path, srs=3035, pixelRes=100, padExtent=0, ) print("Extent (padExtent=0):", ec0.region.extent) print("Mask shape (padExtent=0):", ec0.region.mask.shape)

Extent (padExtent=0): (4037300.00000,3049300.00000,4067700.00000,3100200.00000)
Mask shape (padExtent=0): (509, 304)

In [12]:

# With padding: the grid is expanded, which helps compute exclusions near the region border correctly.
ec_pad = gl.ExclusionCalculator(
    region=region_path,
    srs=3035,
    pixelRes=100,
    padExtent=2000,  # 2 km padding around the region bounding box
)

print("Extent (padExtent=2000):", ec_pad.region.extent)
print("Mask shape (padExtent=2000):", ec_pad.region.mask.shape)

# With padding: the grid is expanded, which helps compute exclusions near the region border correctly. ec_pad = gl.ExclusionCalculator( region=region_path, srs=3035, pixelRes=100, padExtent=2000, # 2 km padding around the region bounding box ) print("Extent (padExtent=2000):", ec_pad.region.extent) print("Mask shape (padExtent=2000):", ec_pad.region.mask.shape)

Extent (padExtent=2000): (4035300.00000,3047300.00000,4069700.00000,3102200.00000)
Mask shape (padExtent=2000): (549, 344)

Note that

• xMin and yMin become smaller
• xMax and yMax become larger

So the box grows outward in every direction, which is why the grid shape also increases.

6. Using value ranges with excludeRasterType¶

excludeRasterType can exclude raster cells based on their values. For categorical rasters such as CORINE Land Cover (CLC), multiple classes can be selected in a single value expression.

Instead of listing every value separately, you can also use compact range syntax.

Here, three groups of CLC classes are selected:

• [1-2] selects all classes from 1 to 2
• 12 selects class 12
• [21-23] selects all classes from 21 to 23

The comma , combines multiple selections. Square brackets [ ] include the endpoints of a range, while parentheses ( ) exclude them. For example, [21-23) includes 21 and 22, but not 23. And vice versa (21-23] selects 22 and 23, but not 21.

In [13]:

ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100)

ec.excludeRasterType(
    source=clc_path,
    value="[1-2],12,[21-23]",
)

ec.draw()

ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100) ec.excludeRasterType( source=clc_path, value="[1-2],12,[21-23]", ) ec.draw()

Memory usage during calc: 325.15234375 MB

Out[13]:

<Axes: >

7. Alternative exclusion approaches (include, invert, initialValue)¶

By default, a new ExclusionCalculator starts with all cells inside the region marked as eligible. With initialValue=False, this can be reversed so that all cells start as unavailable.

This allows eligibility to be approached from two directions. For example, suppose we want only the CLC classes [1-2], 12, and [21-23] to remain eligible.

The first approach starts with nothing available and explicitly includes these classes:

In [14]:

# Start with everything unavailable.
ec = gl.ExclusionCalculator(
    region=region_path,
    srs=3035,
    pixelRes=100,
    initialValue=False,
)

# Include only the selected land cover classes.
ec.excludeRasterType(
    source=clc_path,
    value="[1-2],12,[21-23]",
    mode="include",
)

print("Percent available with include mode:", ec.percentAvailable)
ec.draw()

# Start with everything unavailable. ec = gl.ExclusionCalculator( region=region_path, srs=3035, pixelRes=100, initialValue=False, ) # Include only the selected land cover classes. ec.excludeRasterType( source=clc_path, value="[1-2],12,[21-23]", mode="include", ) print("Percent available with include mode:", ec.percentAvailable) ec.draw()

Memory usage during calc: 309.43359375 MB
Percent available with include mode: 48.132329724853406

Out[14]:

<Axes: >

Alternatively, we can start with everything available and use invert=True to exclude everything except the selected classes:

In [15]:

# Start with everything available.
ec = gl.ExclusionCalculator(
    region=region_path,
    srs=3035,
    pixelRes=100,
)

# Invert the selection, then exclude it.
# This excludes all cells NOT in the selected classes,
# so only the same selected classes remain eligible.
ec.excludeRasterType(
    source=clc_path,
    value="[1-2],12,[21-23]",
    invert=True,
)

print("Percent available with invert=True:", ec.percentAvailable)
ec.draw()

# Start with everything available. ec = gl.ExclusionCalculator( region=region_path, srs=3035, pixelRes=100, ) # Invert the selection, then exclude it. # This excludes all cells NOT in the selected classes, # so only the same selected classes remain eligible. ec.excludeRasterType( source=clc_path, value="[1-2],12,[21-23]", invert=True, ) print("Percent available with invert=True:", ec.percentAvailable) ec.draw()

Memory usage during calc: 312.0859375 MB
Percent available with invert=True: 48.132329724853406

Out[15]:

<Axes: >

Why do these look similar?

Both approaches leave the same CLC classes eligible, but they get there differently:

• initialValue=False + mode="include": start with nothing and add the selected classes.
• invert=True: start with everything and remove everything except the selected classes.

8. Working with Digital Elevation Model (DEM) mosaics for elevation and slope exclusions¶

Terrain can affect whether an area is suitable for a specific use. A Digital Elevation Model (DEM) provides the elevation of the terrain and can also be used to derive additional information, such as slope and orientation.

In this example, we use DEM data to create three terrain-based exclusion criteria:

1. First, multiple DEM tiles are combined into one elevation raster covering the full study area with rasterMosaic().
2. The elevation raster is used directly to exclude areas above 500 m.
3. The DEM mosaic is used to calculate the slope, and areas steeper than 15° are excluded.
4. Finally, the north–south gradient is calculated to identify and exclude north-facing slopes.

This demonstrates how a single DEM can be used both directly as an exclusion raster and as a basis for deriving additional raster layers for an eligibility analysis.

DEM data for other regions can be downloaded from:

• https://dwtkns.com/srtm30m

8.1 Creating the DEM mosaic¶

In [16]:

import os

dem_dir = os.path.dirname(gl._test_data_["aachenShapefile.shp"])
dem_raster = ec.region.extent.rasterMosaic(os.path.join(dem_dir, "*.hgt"))  # .hgt files are used as DEM inputs
print(type(dem_raster))  # multiple DEM files merged into a single raster mosaic
gk.drawRaster(
    dem_raster, cbarTitle="Elevation [m]"
)  # the mosaic is displayed with `gk.drawRaster()` as an elevation heatmap

import os dem_dir = os.path.dirname(gl._test_data_["aachenShapefile.shp"]) dem_raster = ec.region.extent.rasterMosaic(os.path.join(dem_dir, "*.hgt")) # .hgt files are used as DEM inputs print(type(dem_raster)) # multiple DEM files merged into a single raster mosaic gk.drawRaster( dem_raster, cbarTitle="Elevation [m]" ) # the mosaic is displayed with `gk.drawRaster()` as an elevation heatmap

<class 'osgeo.gdal.Dataset'>

Out[16]:

AxHands(ax=<Axes: >, handles=<matplotlib.image.AxesImage object at 0x7a3968835310>, cbar=<matplotlib.colorbar.Colorbar object at 0x7a396869e270>)