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

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

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]:

# Mask: Pixel grid inside bounding box, according to pixelRes (True = inside regiion, False = outside region)
ec.region.mask

# Mask: Pixel grid inside bounding box, according to pixelRes (True = inside regiion, False = outside region) 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)

In [6]:

# 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[6]:

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 [7]:

# Mask: Pixel grid inside bounding box, according to pixelRes (True = inside regiion, False = outside region)
ec.region.mask

# Mask: Pixel grid inside bounding box, according to pixelRes (True = inside regiion, False = outside region) ec.region.mask

Out[7]:

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 [8]:

# 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[8]:

(509, 304)

2. Working with Exclusion Sets¶

Instead of calling several exclusion functions one by one, you can define exclusion rules in a table and apply them in a single step. In this example, a sample exclusion table is loaded from the built-in GLAES test data.

In [9]:

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[9]:

	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

In [10]:

# 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: 244.5859375 MB
Memory usage during calc: 244.890625 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: 266.54296875 MB

Out[10]:

<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 [11]:

## 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: 275.99609375 MB

Out[11]:

<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 [12]:

ec.drawWithSmopyBasemap?

ec.drawWithSmopyBasemap?

In [13]:

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

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

5. Considering region boundary effects with padExtent¶

Some exclusions (for example buffers around roads or settlements) depend on features that may lie just outside the region boundary.

If you do not include a padding area around the region, then:

• buffers can get cut off at the boundary
• exclusions near the border may be underestimated

padExtent adds a buffer (in meters, in the calculator's SRS units) around the region's bounding box for processing, while the final availability is still evaluated for the region.

Below we compare the region grid with and without padding.

In [14]:

# 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 [15]:

# 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. excludeRasterType with ranges for clc classes¶

excludeRasterType can exclude raster cells based on selected values. For categorical rasters such as CORINE Land Cover (CLC), this lets you exclude several classes at once.

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

In the example below:

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

To the value applies the following:

• , separates several selectors, so "[1-2],12,[21-23]" combines three selections in one string.
• Square brackets include the endpoint values: [21-23] selects 21, 22, and 23.
• Parentheses exclude endpoint values: [21-23) selects 21 and 22, but not 23; (21-23] selects 22 and 23, but not 21.

This is useful when you want to exclude multiple related land cover classes in one step.

In [16]:

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: 336.0 MB

Out[16]:

<Axes: >

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

By default, a new ExclusionCalculator starts with all cells inside the region marked as eligible.

You can change this starting state with initialValue:

• the default behavior starts with all cells available
• initialValue=False starts with all cells unavailable

This is useful when you want to build suitability from the ground up by explicitly including only selected areas.

Two related options in excludeRasterType are:

• mode="include": keep only the selected values as eligible
• invert=True: reverse the selection before applying the exclusion

To make this easier to understand, both examples below use the same land cover classes:

• [1-2]
• 12
• [21-23]

In [17]:

# 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: 327.05859375 MB
Percent available with include mode: 48.132329724853406

Out[17]:

<Axes: >

In [18]:

# 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: 331.078125 MB
Percent available with invert=True: 48.132329724853406

Out[18]:

<Axes: >

Why do these look similar?
Both examples keep the same classes as eligible, but they start from different initial states:

• In the mode="include" example, initialValue=False means all cells start as unavailable, and only the selected classes are added as eligible.
• In the invert=True example, the default initial state means all cells start as available, and then everything except the selected classes is excluded.

8. Working with Digital Elevation Model (DEM) to exclude Raster¶

By default, DEM data is often distributed as separate raster tiles. A Digital Elevation Model (DEM) is a raster dataset in which each pixel stores the ground elevation at that location.

If a study area covers multiple DEM tiles, they can be merged into a single raster with rasterMosaic().

This is useful when you want to create one continuous elevation layer for the full region.

In the example below:

• .hgt files are used as DEM inputs
• all matching files are merged into one raster mosaic
• the mosaic is displayed with gk.drawRaster() as an elevation heatmap

DEM data for other regions can be downloaded from:

• https://dwtkns.com/srtm30m

In [19]:

import os

dem_dir = os.path.dirname(gl._test_data_["aachenShapefile.shp"])
dem_raster = ec.region.extent.rasterMosaic(os.path.join(dem_dir, "*.hgt"))
print(type(dem_raster))  # multiple DEM files merged into a single raster mosaic
gk.drawRaster(dem_raster, cbarTitle="Elevation [m]")

import os dem_dir = os.path.dirname(gl._test_data_["aachenShapefile.shp"]) dem_raster = ec.region.extent.rasterMosaic(os.path.join(dem_dir, "*.hgt")) print(type(dem_raster)) # multiple DEM files merged into a single raster mosaic gk.drawRaster(dem_raster, cbarTitle="Elevation [m]")

<class 'osgeo.gdal.Dataset'>

Out[19]:

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

In [20]:

# exclude areas above 500m elevation. Note, that here the hilly Eiffel 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 Eiffel 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: 400.5546875 MB

Out[20]:

<Axes: >

In [21]:

# 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: 372.32421875 MB

Out[21]:

<Axes: >

In [22]:

# 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: 459.36328125 MB

Out[22]:

<Axes: >

9. ExcludeVectorType¶

This section shows how excludeVectorType() can be used to exclude areas from vector datasets.

Use where clause to exclude IUCN_CAT='V' (IUCN protected area management category, category V: Protected Landscape):

In [23]:

#  Exclude Vector - Simple source w/ where clause

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 - Simple source w/ where clause 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[23]:

<Axes: >

The following example uses resolutionDiv to rasterize vector features at a finer internal resolution.

In [24]:

# 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,
    intermediate="tmp/wdpa_V.tif",  # stores the rasterized exclusion result for these exact arguments, so the same exclusion can be reused later instead of recalculated
)
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, intermediate="tmp/wdpa_V.tif", # stores the rasterized exclusion result for these exact arguments, so the same exclusion can be reused later instead of recalculated ) ec.draw()

Computing intermediate exclusion file: tmp/wdpa_V.tif

Memory usage during calc: 438.11328125 MB

Out[24]:

<Axes: >

ec.save() exports the current availability result of the ExclusionCalculator as a GeoTIFF raster

In [25]:

ec.save?

ec.save?

In [26]:

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

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

Out[26]:

'tmp/test_resolution_div_8.tif'

10 Visualization Methods¶

In [27]:

# 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[27]:

<Axes: >

Additional spatial layers can also be drawn on top of the exclusion result. In the example below, major roads are added as a vector overlay.

In [28]:

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 drawing the roads, we inspect the available road categories in the type column. These values can then be used to build the where filter for selecting only major roads.

In [29]:

roads.type.unique()

roads.type.unique()

Out[29]:

<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

In [30]:

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()

drawWithSmopyBasemap() (cf. section 3) can be used to visualize roads in its geographic context.

In [31]:

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()