<jupyter_start><jupyter_text>Siting First Aid Stations with LSCP on Toronto's University CampusThis notebook implements the location set covering problem (LSCP) to site first aid stations for Toronto's University campus. LSCP determines locations such that as many demand points as possible are within the maximum service radius. By doing so, the number of locations required to cover all demand points is minimized. In contrast to the maximum covering location problem it does not require a maximum amount of locations as an input parameter. Hence, the use cases for LSCP are specific. One common scenario is finding suitable locations for hospitals in cities where the amount is governed by reachability preconditions to households in a certain amount of time.  In this notebook demand is described as buildings within Toronto's University campus. Possible locations for the first aid stations are determined using [Overpass API](https://wiki.openstreetmap.org/wiki/Overpass_API) which fetches geometries directly from [OpenStreetMap.org](https://openstreetmap.org).The routing matrix is computed using [Open Source Routing Machine's demo server API](http://project-osrm.org). PrerequisitesTo run this jupyter notebook locally, you will have to install the following additional libraries into your python environment.- [shapely](https://github.com/shapely/shapely)- [folium](https://python-visualization.github.io/folium)- [overpy](https://github.com/DinoTools/python-overpy)- [routing-py](https://github.com/gis-ops/routing-py)<jupyter_code>import random
import json

import numpy as np
from spopt.locate.coverage import LSCP
import pulp

from shapely.geometry import shape, Point
import folium
from routingpy import OSRM
import overpy


api = overpy.Overpass()

# First of all we want to fetch all buildings in the University of Toronto boundary
result = api.query(
    """
    [out:json];
    
    area[name="Toronto"]->.b;
    way(area.b)[name="University of Toronto"];
    map_to_area -> .a;
    
    way[building="university"](area.a);
    out center;
    """
)

m = folium.Map(
    location=[43.6624674, -79.3988052], tiles="cartodbpositron", zoom_start=15
)

# Let's add our results to the map
university_buildings = []
university_buildings_fg = folium.FeatureGroup("University Buildings").add_to(m)
for way in result.ways:
    folium.CircleMarker(
        location=[way.center_lat, way.center_lon],
        radius=3,
        fill=True,
        fill_opacity=0.3,
        popup=folium.Popup(way.tags.get("name", "n/a"), show=False),
        color="blue",
    ).add_to(university_buildings_fg)

    university_buildings.append(
        {
            "name": way.tags.get("name", "n/a"),
            "center": [way.center_lon, way.center_lat],
        }
    )

m


# And now we want to fetch some roads in the same area and 
# return their center points which we can use for our 
# potential locations for first aid stations
result_roads = api.query(
    """
    [out:json];
    
    area[name="Toronto"]->.b;
    way(area.b)[name="University of Toronto"];
    map_to_area -> .a;
    
    way["highway"~"secondary|tertiary|unclassified"](area.a);
    
    out center;
    """
)

# Similarly we want to add them to the map
first_aid_stations = []
first_aid_stations_fg = folium.FeatureGroup("First Aid Station Locations").add_to(m)
for way in result_roads.ways:
    folium.CircleMarker(
        location=[way.center_lat, way.center_lon],
        radius=2,
        fill=True,
        fill_opacity=1,
        popup=folium.Popup(way.tags.get("name", "n/a"), show=False),
        color="black",
    ).add_to(first_aid_stations_fg)

    first_aid_stations.append(
        {
            "name": way.tags.get("name", "n/a"),
            "center": [way.center_lon, way.center_lat],
        }
    )

m<jupyter_output><empty_output><jupyter_text>Computing the Distance (Routing) MatrixGiven our potential first aid station locations and the university buildings on the campus we want to compute our distance matrix. Instead of using simple euclidean distance we can make use of proper routes following roads which will yield more realistic travel times. The [routing-py library](https://github.com/gis-ops/routing-py) helps us achieve this which features the ability to use different open source frameworks such as OSRM or Valhalla as well as 3rd party API's such as Google Maps or HERE.<jupyter_code>client = OSRM(base_url="https://router.project-osrm.org")


university_buildings_locations = [bld["center"] for bld in university_buildings]
first_aid_stations_locations = [fab["center"] for fab in first_aid_stations]

university_buildings_locations + first_aid_stations_locations

sources_indices = list(range(0, len(university_buildings_locations)))
targets_indices = list(
    range(
        len(university_buildings_locations),
        len(university_buildings_locations) + len(first_aid_stations_locations),
    )
)

osrm_routing_matrix = client.matrix(
    dry_run=False,
    locations=university_buildings_locations + first_aid_stations_locations,
    sources=sources_indices,
    destinations=targets_indices,
    profile="pedestrian",
)
cost_matrix = np.array(osrm_routing_matrix.durations)<jupyter_output><empty_output><jupyter_text>Solving the LSCPLast but not least we want to run the solver. We choose the Pulp CBC solver and set the distance cutoff to 2 minutes afoot. Finally we loop over the solution and add the sited stations and the corresponding allocated buildings to the map with different radiuses.<jupyter_code>solver = pulp.PULP_CBC_CMD(msg=False)  # see solvers available in pulp reference
MAX_COVERAGE = 120  # maximum acceptable service duration, 2 minutes

lscp_from_cost_matrix = LSCP.from_cost_matrix(cost_matrix, MAX_COVERAGE)
result = lscp_from_cost_matrix.solve(solver)
lscp_from_cost_matrix.facility_client_array()

# we init a dict here holding the mappings restaurant to households
solution = {}
colors = [
    "green",
    "red",
    "blue",
    "orange",
    "purple",
    "pink",
    "black",
    "brown",
    "teal",
    "magenta"
]


sited_first_aid_stations_fg = folium.FeatureGroup("Sited First Aid Stations").add_to(m)
allocated_buildings_demand_fg = folium.FeatureGroup("Allocated Buildings").add_to(m)


for i in range(len(first_aid_stations)):
    if lscp_from_cost_matrix.fac2cli[i]:
        color = random.choice(colors)

        folium.CircleMarker(
            location=[
                first_aid_stations[i]["center"][1],
                first_aid_stations[i]["center"][0],
            ],
            radius=30,
            fill=True,
            popup=folium.Popup(f"Sited First Aid Station (index: {i})", show=True),
            color=color,
        ).add_to(sited_first_aid_stations_fg)
        solution[i] = []
        for j in lscp_from_cost_matrix.fac2cli[i]:
            solution[i].append(university_buildings[j])
            folium.CircleMarker(
                location=[
                    university_buildings[j]["center"][1],
                    university_buildings[j]["center"][0],
                ],
                radius=10,
                fill=True,
                popup=f"Building allocated to First Aid Station with index: {i}",
                color=color,
            ).add_to(allocated_buildings_demand_fg)

folium.LayerControl().add_to(m)

m
solution<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0962087240114146, 'fraction_numerical': 0.005707297187117815, 'mean_word_length': 3.6908859145952837, 'pattern_counts': {'":': 4, '<': 13, '<?xml version=': 0, '>': 17, 'https://': 8, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '21053', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2286, 'n_tokens_neox': 2189, 'n_words': 617, 'path': 'notebooks/lscp_gis.ipynb', 'repo_name': 'gis-ops/spopt', 'size': '7359'}