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: 241.69140625 MB
Memory usage during calc: 242.14453125 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: 265.5703125 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: 274.7265625 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: 332.36328125 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: 333.59765625 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: 336.8203125 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 0x7a4b5705ae90>, cbar=<matplotlib.colorbar.Colorbar object at 0x7a4b56fa06e0>)

8.2 Excluding areas above 500 m¶

In [17]:

# exclude areas above 500m elevation. Note, that here the hilly Eifel region is mainly exluded in the region's south.
ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100)
ec.excludeRasterType(dem_raster, value="(500-]", resampleAlg="mode")
ec.draw()

# exclude areas above 500m elevation. Note, that here the hilly Eifel region is mainly exluded in the region's south. ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100) ec.excludeRasterType(dem_raster, value="(500-]", resampleAlg="mode") ec.draw()

Memory usage during calc: 405.265625 MB

Out[17]:

<Axes: >

8.3 Excluding slopes above 15°¶

In [18]:

# Exclude Raster - DEM Gradient Mosaic. Exlude slopes above 15 degrees.

dem_gradient_raster = gk.raster.mutateRaster(
    source=gk.raster.gradient(  # computes the gradient (slope) of a raster
        source=dem_raster, factor="latlonToM"
    ),  # factor to convert lat/lon degrees to meters
    processor=lambda x: np.degrees(np.arctan(x)),  # convert gradient to degrees
)

ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100)
ec.excludeRasterType(dem_gradient_raster, value=(15, None), resampleAlg="mode")

ec.draw()

# Exclude Raster - DEM Gradient Mosaic. Exlude slopes above 15 degrees. dem_gradient_raster = gk.raster.mutateRaster( source=gk.raster.gradient( # computes the gradient (slope) of a raster source=dem_raster, factor="latlonToM" ), # factor to convert lat/lon degrees to meters processor=lambda x: np.degrees(np.arctan(x)), # convert gradient to degrees ) ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100) ec.excludeRasterType(dem_gradient_raster, value=(15, None), resampleAlg="mode") ec.draw()

Memory usage during calc: 377.03125 MB

Out[18]:

<Axes: >

8.4 Excluding north-facing slopes¶

In [19]:

# Exlude slopes with north orientation, as north facing slopes receive less radiation.

dem_ns_gradient_raster = gk.raster.mutateRaster(
    source=gk.raster.gradient(source=dem_raster, factor="latlonToM", mode="north-south"),
    processor=lambda matrix: np.degrees(np.arctan(matrix)),
)

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

# exclude all north facing slopes (values > 0 degrees)
ec.excludeRasterType(dem_ns_gradient_raster, value=(0, None), resampleAlg="mode")

ec.draw()

# Exlude slopes with north orientation, as north facing slopes receive less radiation. dem_ns_gradient_raster = gk.raster.mutateRaster( source=gk.raster.gradient(source=dem_raster, factor="latlonToM", mode="north-south"), processor=lambda matrix: np.degrees(np.arctan(matrix)), ) ec = gl.ExclusionCalculator(region_path, srs=3035, pixelRes=100) # exclude all north facing slopes (values > 0 degrees) ec.excludeRasterType(dem_ns_gradient_raster, value=(0, None), resampleAlg="mode") ec.draw()

Memory usage during calc: 398.828125 MB

Out[19]:

<Axes: >

9. Exclude vector type¶

Vector datasets often contain multiple features and attribute information that can be used to define more specific exclusion criteria. The excludeVectorType() method provides several options to control which vector features are considered and how they are processed during the exclusion.

In the following examples, we explore several options of excludeVectorType():

1. The where argument is used to select only features that match a specific attribute condition.

2. resolutionDiv is used to rasterize vector features at a finer internal resolution.

3. The intermediate argument is used to store and reuse the rasterized result of an exclusion operation. And ec.save exports the combined results of the ExclusionCalculator.

These examples demonstrate how excludeVectorType() can be adapted when more control over the selection and rasterization of vector features is required.

In [20]:

ec.excludeVectorType?

ec.excludeVectorType?

9.1 Exclude vector - dsmple source with where clause¶

The where argument allows the exclusion to be restricted to features matching specific attribute values. It filters the features in the input dataset before applying the exclusion.

In this example, the source WDPA_aachen_bbox.shp contains protected areas from the World Database on Protected Areas (WDPA). Only features with IUCN_CAT='V', corresponding to IUCN protected area management category V (Protected Landscape/Seascape), are excluded.

In [21]:

# Exclude vector - where argument

ec = gl.ExclusionCalculator(region=gk._test_data_["aachenShapefile.shp"], srs="LAEA", pixelRes=50, padExtent=5000)

ec.excludeVectorType(
    source=r"data/WDPA_aachen_bbox.shp",
    where="IUCN_CAT='V'",
)

ec.draw()

# Exclude vector - where argument ec = gl.ExclusionCalculator(region=gk._test_data_["aachenShapefile.shp"], srs="LAEA", pixelRes=50, padExtent=5000) ec.excludeVectorType( source=r"data/WDPA_aachen_bbox.shp", where="IUCN_CAT='V'", ) ec.draw()

Out[21]:

<Axes: >

9.2 Exclude vector - representation of fine details¶

The resolutionDiv parameter can be used to rasterize vector features at a finer internal resolution before the result is transferred to the regular exclusion grid. This can improve the representation of vector boundaries and buffered geometries.

In the following example, resolutionDiv=8 is applied to the same category V protected areas, together with a 400 m buffer. The rasterized exclusion is stored using the intermediate argument, allowing it to be reused for the same exclusion parameters without repeating the rasterization step.

In [22]:

# Exclude vector - resolutionDiv
ec = gl.ExclusionCalculator(region=gk._test_data_["aachenShapefile.shp"], srs="LAEA", pixelRes=50, padExtent=5000)

ec.excludeVectorType(
    source=r"data/WDPA_aachen_bbox.shp",
    where="IUCN_CAT='V'",
    buffer=400,
    resolutionDiv=8,  # Region Mask's native resolution divided by factor 8, meaning a 1000 m x 1000 m resolution becomes 125 m x 125 m
    intermediate="tmp/wdpa_V.tif",  # intermediate stores the rasterized result of this particular exclusion
)
ec.draw()

# Exclude vector - resolutionDiv ec = gl.ExclusionCalculator(region=gk._test_data_["aachenShapefile.shp"], srs="LAEA", pixelRes=50, padExtent=5000) ec.excludeVectorType( source=r"data/WDPA_aachen_bbox.shp", where="IUCN_CAT='V'", buffer=400, resolutionDiv=8, # Region Mask's native resolution divided by factor 8, meaning a 1000 m x 1000 m resolution becomes 125 m x 125 m intermediate="tmp/wdpa_V.tif", # intermediate stores the rasterized result of this particular exclusion ) ec.draw()

Computing intermediate exclusion file: tmp/wdpa_V.tif

Memory usage during calc: 429.97265625 MB

Out[22]:

<Axes: >

9.3 Store, save, and export¶

The intermediate argument stores the rasterized result of this particular exclusion in tmp/wdpa_V.tif (example 9.2 above). For the same exclusion arguments, this intermediate result can be reused instead of rasterizing the vector data again.

Note that intermediate and ec.save() serve different purposes:

• intermediate stores the rasterized result of this particular excludeVectorType() operation.
• ec.save() exports the current combined availability result of the ExclusionCalculator as GEOTIFF raster.

In [23]:

ec.save?

ec.save?

In [24]:

ec.save("tmp/test_resolution_div_8.tif")

ec.save("tmp/test_resolution_div_8.tif")

Out[24]:

'tmp/test_resolution_div_8.tif'

10. Visualization methods¶

The ExclusionCalculator provides several options for visualizing the resulting eligibility map. Besides adjusting the appearance of the exclusion result itself, additional spatial data can be added to provide more context.

The following examples show how to customize the colors and legend of ec.draw(), overlay additional vector features, and combine these layers with a basemap.

In [25]:

ec.draw?

ec.draw?

In [26]:

# Regular drawing with custom colors and legend position
ec.draw(
    goodColor="blue",
    excludedColor=(1.0, 0.2, 0.0, 0.6),
    legendargs={"loc": "upper right"},
)

# Regular drawing with custom colors and legend position ec.draw( goodColor="blue", excludedColor=(1.0, 0.2, 0.0, 0.6), legendargs={"loc": "upper right"}, )

Out[26]:

<Axes: >

Additional spatial layers can be drawn on top of the exclusion result.

In the following example, major roads are added as a vector overlay. First, the road features intersecting the study region are extracted from the example dataset.

In [27]:

road_path = gl._test_data_["aachenRoads.shp"]

roads = gk.vector.extractFeatures(
    road_path,
    geom=ec.region.geometry,
)

road_path = gl._test_data_["aachenRoads.shp"] roads = gk.vector.extractFeatures( road_path, geom=ec.region.geometry, )

Before selecting the major roads, we inspect the available road categories in the type column. These values can then be used to construct a where filter.

In [28]:

roads.type.unique()

roads.type.unique()

Out[28]:

<StringArray>
[     'secondary',    'residential',       'tertiary',        'primary',
       'motorway',   'primary_link',  'motorway_link',  'tertiary_link',
 'secondary_link',           'path',           'road',      'rest_area',
        'raceway',          'steps',          'track',       'proposed',
        'service']
Length: 17, dtype: str

Based on the available categories, motorway, primary, and trunk roads are selected and drawn on top of the exclusion result.

ec.draw() returns a Matplotlib Axes object, which can be passed to gk.drawGeoms() to add the selected road geometries to the same plot.

In [29]:

ax = ec.draw()
major_roads = gk.vector.extractFeatures(
    road_path,
    where="type='motorway' OR type='primary' OR type='trunk'",
)
gk.drawGeoms(major_roads, ax=ax, color="k", linestyle="--", srs=ec.region.srs)
plt.show()

ax = ec.draw() major_roads = gk.vector.extractFeatures( road_path, where="type='motorway' OR type='primary' OR type='trunk'", ) gk.drawGeoms(major_roads, ax=ax, color="k", linestyle="--", srs=ec.region.srs) plt.show()

The same vector overlay can also be combined with drawWithSmopyBasemap() to show the exclusion result and the selected roads in their geographic context.

Since the basemap uses Web Mercator, the road geometries are drawn using EPSG:3857.

In [30]:

ax = ec.drawWithSmopyBasemap(
    zoom=9,
    legendargs={"loc": "upper right"},
)
gk.drawGeoms(major_roads, ax=ax, color="k", linestyle="--", srs=gk.srs.EPSG3857)
plt.show()

ax = ec.drawWithSmopyBasemap( zoom=9, legendargs={"loc": "upper right"}, ) gk.drawGeoms(major_roads, ax=ax, color="k", linestyle="--", srs=gk.srs.EPSG3857) plt.show()