This page was generated from notebooks/mclp_gis.ipynb. Interactive online version:

Siting Restaurants with MCLP in San Francisco

Authors: Timothy Ellersiek, James Gaboardi

This notebook implements the maximize covering location problem (MCLP) to site 2 of 10 possible locations for restaurants in San Francisco’s Mission District & Bernal Heights. The objective is to maximize the demand of randomly created households able to reach the restaurants within 3 minutes afoot.

The routing matrix is computed using Open Source Routing Machine’s demo server API.

Prerequisites

To run this jupyter notebook locally, you will have to install the following additional libraries into your python environment.

• folium

• routing-py

• shapely

%config InlineBackend.figure_format = "retina"
%load_ext watermark
%watermark

Last updated: 2023-12-10T13:58:07.138829-05:00

Python implementation: CPython
Python version       : 3.11.6
IPython version      : 8.18.0

Compiler    : Clang 15.0.7
OS          : Darwin
Release     : 23.1.0
Machine     : x86_64
Processor   : i386
CPU cores   : 8
Architecture: 64bit

import geopandas
import folium
import json
import numpy
import pulp
import routingpy
from routingpy import OSRM
import shapely
from shapely.geometry import shape, Point
import spopt
from spopt.locate import MCLP, simulated_geo_points

%watermark -w
%watermark -iv

Watermark: 2.4.3

pulp     : 2.7.0
spopt    : 0.5.1.dev59+g343ef27
shapely  : 2.0.2
routingpy: 1.2.1
folium   : 0.15.0
json     : 2.0.9
geopandas: 0.14.1
numpy    : 1.26.2

---

San Francisco’s Mission District and Bernal Heights merged polygon

SF_MISSION_BERNAL = {
    "type": "Polygon",
    "coordinates": [
        [
            [-122.405115, 37.764635],
            [-122.405212, 37.763469],
            [-122.405832, 37.762199],
            [-122.405859, 37.761815],
            [-122.40604, 37.761865],
            [-122.406365, 37.761199],
            [-122.40645, 37.760114],
            [-122.406167, 37.759315],
            [-122.406015, 37.759409],
            [-122.405604, 37.758851],
            [-122.40405, 37.757619],
            [-122.403428, 37.756871],
            [-122.403585, 37.756818],
            [-122.403328, 37.756082],
            [-122.40315, 37.754488],
            [-122.403524, 37.754463],
            [-122.403437, 37.753199],
            [-122.403364, 37.752366],
            [-122.403022, 37.752336],
            [-122.403052, 37.752048],
            [-122.405091, 37.745281],
            [-122.405614, 37.7442],
            [-122.405493, 37.744084],
            [-122.404865, 37.744385],
            [-122.404524, 37.744285],
            [-122.406652, 37.741241],
            [-122.406936, 37.740554],
            [-122.407045, 37.739587],
            [-122.408136, 37.73964],
            [-122.408009, 37.737734],
            [-122.408218, 37.73765],
            [-122.408284, 37.736846],
            [-122.408613, 37.736074],
            [-122.409287, 37.735258],
            [-122.410052, 37.734675],
            [-122.411978, 37.733732],
            [-122.412695, 37.733197],
            [-122.414554, 37.732371],
            [-122.416255, 37.732034],
            [-122.419881, 37.732016],
            [-122.423718, 37.73155],
            [-122.425944, 37.731698],
            [-122.425256, 37.732125],
            [-122.422274, 37.732383],
            [-122.421889, 37.732533],
            [-122.421787, 37.732774],
            [-122.422472, 37.733967],
            [-122.422901, 37.734388],
            [-122.422243, 37.734846],
            [-122.422133, 37.735162],
            [-122.426983, 37.735455],
            [-122.427849, 37.734718],
            [-122.428578, 37.735082],
            [-122.42801, 37.735441],
            [-122.428454, 37.735882],
            [-122.42551, 37.737774],
            [-122.42452, 37.739868],
            [-122.42427, 37.739867],
            [-122.424394, 37.740405],
            [-122.424168, 37.740805],
            [-122.426453, 37.764634],
            [-122.421885, 37.764908],
            [-122.42173, 37.763304],
            [-122.421031, 37.76335],
            [-122.421248, 37.764947],
            [-122.407553, 37.765798],
            [-122.407431, 37.764497],
            [-122.405115, 37.764635],
        ]
    ],
}

SF_MISSION_BERNAL = str(SF_MISSION_BERNAL).replace("'", '"')
SF_MISSION_BERNAL

'{"type": "Polygon", "coordinates": [[[-122.405115, 37.764635], [-122.405212, 37.763469], [-122.405832, 37.762199], [-122.405859, 37.761815], [-122.40604, 37.761865], [-122.406365, 37.761199], [-122.40645, 37.760114], [-122.406167, 37.759315], [-122.406015, 37.759409], [-122.405604, 37.758851], [-122.40405, 37.757619], [-122.403428, 37.756871], [-122.403585, 37.756818], [-122.403328, 37.756082], [-122.40315, 37.754488], [-122.403524, 37.754463], [-122.403437, 37.753199], [-122.403364, 37.752366], [-122.403022, 37.752336], [-122.403052, 37.752048], [-122.405091, 37.745281], [-122.405614, 37.7442], [-122.405493, 37.744084], [-122.404865, 37.744385], [-122.404524, 37.744285], [-122.406652, 37.741241], [-122.406936, 37.740554], [-122.407045, 37.739587], [-122.408136, 37.73964], [-122.408009, 37.737734], [-122.408218, 37.73765], [-122.408284, 37.736846], [-122.408613, 37.736074], [-122.409287, 37.735258], [-122.410052, 37.734675], [-122.411978, 37.733732], [-122.412695, 37.733197], [-122.414554, 37.732371], [-122.416255, 37.732034], [-122.419881, 37.732016], [-122.423718, 37.73155], [-122.425944, 37.731698], [-122.425256, 37.732125], [-122.422274, 37.732383], [-122.421889, 37.732533], [-122.421787, 37.732774], [-122.422472, 37.733967], [-122.422901, 37.734388], [-122.422243, 37.734846], [-122.422133, 37.735162], [-122.426983, 37.735455], [-122.427849, 37.734718], [-122.428578, 37.735082], [-122.42801, 37.735441], [-122.428454, 37.735882], [-122.42551, 37.737774], [-122.42452, 37.739868], [-122.42427, 37.739867], [-122.424394, 37.740405], [-122.424168, 37.740805], [-122.426453, 37.764634], [-122.421885, 37.764908], [-122.42173, 37.763304], [-122.421031, 37.76335], [-122.421248, 37.764947], [-122.407553, 37.765798], [-122.407431, 37.764497], [-122.405115, 37.764635]]]}'

Convert SF_MISSION_BERNAL to a shapely.Polygon then a geopandas.GeoDataFrame

misson_bernal = shape(json.loads(SF_MISSION_BERNAL))
misson_bernal

misson_bernal = geopandas.GeoDataFrame(geometry=[misson_bernal])
misson_bernal

	geometry
0	POLYGON ((-122.40511 37.76463, -122.40521 37.7...

Generate 200 artificial households within misson_bernal

NUM_HOUSEHOLDS = 200
households = simulated_geo_points(misson_bernal, needed=NUM_HOUSEHOLDS, seed=0)
households = households.geometry.map(lambda pt: [pt.x, pt.y]).to_list()
households[:5]

[[-122.41455252209364, 37.75604380541952],
 [-122.41317377916111, 37.75021115925127],
 [-122.4177510779481, 37.7536705815843],
 [-122.41739502122897, 37.76209144173078],
 [-122.40834467492677, 37.74966359321169]]

Adding a Leaflet Map for Visualisation via Folium

In the following we are using folium to visualize the data on the map. After creating the base map, we will load the merged polygon of Mission District and Bernal Heights boundary onto it and finally the artificial househoulds we created earlier. In the subsequent step we will add some fixed positions for our possible restaurants which are depicted as red circles.

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

folium.GeoJson(
    data=SF_MISSION_BERNAL,
    name="Mission District & Bernal Heights",
    style_function=lambda x: {"color": "darksalmon", "fillColor": "darksalmon"},
).add_to(m)

households_fg = folium.FeatureGroup("Households").add_to(m)
for coord in households:
    folium.CircleMarker(
        location=[coord[1], coord[0]],
        radius=3,
        fill=True,
        fill_opacity=1,
        popup=folium.Popup("Random Household", show=False),
        color="cornflowerblue",
    ).add_to(households_fg)
m

Make this Notebook Trusted to load map: File -> Trust Notebook

Define 10 restaurant locations

restaurants = [
    [-122.41957, 37.76032],
    [-122.40789, 37.75781],
    [-122.41525, 37.74622],
    [-122.42265, 37.75189],
    [-122.41890, 37.74224],
    [-122.41365, 37.73596],
    [-122.40846, 37.76454],
    [-122.41489, 37.75905],
    [-122.42365, 37.74082],
    [-122.40667, 37.74899],
]
n_rests = len(restaurants)

Add them to the map

restaurants_fg = folium.FeatureGroup("Restaurants").add_to(m)
for i, coord in enumerate(restaurants):
    folium.CircleMarker(
        location=[coord[1], coord[0]],
        radius=5,
        fill=True,
        fill_opacity=1,
        popup=folium.Popup(f"Potential Restaurant ({i})", show=False),
        color="crimson",
    ).add_to(restaurants_fg)
m

Make this Notebook Trusted to load map: File -> Trust Notebook

Computing the Distance (Routing) Matrix

Given our 10 restaurants and 200 artificial households 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 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.

all_locations = households + restaurants
source_indices = [i for i in range(NUM_HOUSEHOLDS)]
target_indices = [i for i in range(NUM_HOUSEHOLDS, len(all_locations))]

client = OSRM(base_url="https://router.project-osrm.org")
osrm_routing_matrix = client.matrix(
    dry_run=False,
    locations=all_locations,
    sources=source_indices,
    destinations=target_indices,
    profile="pedestrian",
)
cost_matrix = numpy.array(osrm_routing_matrix.durations)

Solving the MCLP

Last but not least we want to run the solver. We choose the pulp.COIN_CMD solver and set a maximium coverage to 3 minutes. Our weights of the households we set to 1 and treat them equally. Finally we loop over the solution and add the sited restaurants and the corresponding allocated households to the map with different radiuses.

# maximum acceptable service duration, 3
SERVICE_RADIUS = 180

# weights, all set to 1 for demos sake
ai = numpy.ones((1, NUM_HOUSEHOLDS), dtype=int)

# the number of restaurants to site
P_FACILITIES = 2

# set the solver (for available solvers run ``pulp.listSolvers(onlyAvailable=True)``)
solver = pulp.COIN_CMD(msg=False)

mclp = MCLP.from_cost_matrix(cost_matrix, ai, SERVICE_RADIUS, p_facilities=P_FACILITIES)
mclp = mclp.solve(solver)

How did the model perform?

obj_val = int(mclp.problem.objective.value())
minutes = int(SERVICE_RADIUS/60)
obj_perc = mclp.perc_cov
print(
    f"Of the {NUM_HOUSEHOLDS} households, {obj_val} have access to at least 1 "
    f"restaurant within a {minutes}-minute walk along the street network when "
    f"siting {P_FACILITIES} restaurants. This accounts for {obj_perc}% coverage."
)

Of the 200 households, 154 have access to at least 1 restaurant within a 3-minute walk along the street network when siting 2 restaurants. This accounts for 77.0% coverage.

Add the results to the map

sited_restaurants_fg = folium.FeatureGroup("Sited Restaurants").add_to(m)
allocated_households_demand_fg = folium.FeatureGroup("Allocated Households").add_to(m)
solution = {}
colors = [
    "darkcyan",
    "saddlebrown",
    "lavender",
    "lightskyblue",
    "mediumseagreen",
    "coral",
    "blueviolet",
    "darkgoldenrod",
    "mediumvioletred",
    "mediumorchid"
]
for i in range(n_rests):
    if mclp.fac2cli[i]:
        folium.CircleMarker(
            location=[restaurants[i][1], restaurants[i][0]],
            radius=30,
            fill=True,
            popup=folium.Popup(f"Sited Restaurant ({i})", show=True),
            color=colors[i],
        ).add_to(sited_restaurants_fg)
        solution[i] = []
        for j in mclp.fac2cli[i]:
            solution[i].append(households[j])
            folium.CircleMarker(
                location=[households[j][1], households[j][0]],
                radius=10,
                fill=True,
                popup="Household",
                color=colors[i],
            ).add_to(allocated_households_demand_fg)
folium.LayerControl().add_to(m)
m

Make this Notebook Trusted to load map: File -> Trust Notebook

Inspect the solution set

solution

{4: [[-122.41524164428937, 37.74575134212079],
  [-122.41939041609294, 37.746517470353716],
  [-122.4232014632692, 37.735965467842085],
  [-122.42051691306378, 37.74400636648324],
  [-122.4140060513352, 37.74657122463306],
  [-122.42323994562125, 37.737074528368524],
  [-122.41188716363442, 37.74022473080378],
  [-122.4166609618889, 37.73992108767487],
  [-122.42451537332036, 37.735330127834594],
  [-122.42355414116491, 37.74417809964479],
  [-122.4213505852644, 37.735666491828425],
  [-122.42100984111212, 37.73561618691874],
  [-122.42045162186741, 37.74573767903614],
  [-122.4140979332363, 37.74063905928569],
  [-122.41520587274562, 37.73476727461245],
  [-122.42118193788865, 37.73782393776602],
  [-122.41125572058564, 37.740797233065884],
  [-122.41070210763999, 37.74173662069288],
  [-122.41170660507584, 37.74148457789174],
  [-122.41401197674212, 37.751786210327225],
  [-122.41382634618964, 37.73969735334643],
  [-122.40991900141367, 37.74223349220518],
  [-122.41840106253402, 37.7387367287149],
  [-122.42276619176752, 37.74026120078767],
  [-122.41892417399505, 37.73770106674906],
  [-122.41163237344158, 37.740568264758494],
  [-122.42039964169358, 37.744682871447594],
  [-122.41100854960864, 37.738930706925785],
  [-122.42208826738184, 37.73885550860088],
  [-122.42327589475202, 37.7460946279338],
  [-122.41901571198258, 37.74742653113364],
  [-122.41535585953757, 37.73607306850611],
  [-122.41025793398745, 37.74511425270174],
  [-122.41412809295414, 37.73782696783066],
  [-122.41970410021572, 37.736623528205754],
  [-122.41915845180576, 37.73623856050469],
  [-122.41551073208204, 37.73923240960833],
  [-122.42006414249477, 37.73432756040431],
  [-122.41817054462335, 37.73950355490107],
  [-122.41140314623527, 37.73995333620304],
  [-122.42166950005372, 37.73605302291012],
  [-122.42187930857803, 37.74717190213178],
  [-122.42133239235675, 37.74456173838566],
  [-122.42194819694559, 37.74435235234791],
  [-122.41356130372293, 37.740893604514234],
  [-122.4191004858632, 37.7382987149878],
  [-122.41531864893646, 37.742057632289225],
  [-122.42148375412637, 37.73596321664609],
  [-122.41468049513983, 37.747198302368524],
  [-122.41007117187876, 37.745215819217684],
  [-122.41069845878158, 37.742773768312524],
  [-122.42244404189351, 37.737048133592474],
  [-122.41949619370845, 37.743766497582506],
  [-122.42022883226554, 37.73952102093955],
  [-122.42292601164742, 37.740221292002616],
  [-122.41812421353865, 37.73713086318324],
  [-122.41225380452875, 37.7483419775086],
  [-122.41535766353438, 37.746074199675945],
  [-122.41440239836186, 37.741380940454206],
  [-122.4105207768065, 37.74575801806093],
  [-122.42263298457125, 37.74348609136638],
  [-122.42220382143077, 37.7351770739924],
  [-122.42161560269948, 37.744531940576856],
  [-122.42250059810968, 37.73743562493618],
  [-122.41709590262614, 37.74197743401586],
  [-122.41446908756086, 37.74085356398497],
  [-122.41693866931666, 37.745307885159676],
  [-122.42222954548635, 37.74887491191325],
  [-122.42064590761126, 37.74432569801819],
  [-122.42209591165009, 37.74682924903755],
  [-122.42334097025606, 37.743252477037885],
  [-122.41468376365637, 37.741231897645584],
  [-122.41267626306029, 37.74110751909093],
  [-122.4232176302866, 37.735512604334986]],
 7: [[-122.41455252209364, 37.75604380541952],
  [-122.41317377916111, 37.75021115925127],
  [-122.4177510779481, 37.7536705815843],
  [-122.41739502122897, 37.76209144173078],
  [-122.40834467492677, 37.74966359321169],
  [-122.41406105319669, 37.76324983366825],
  [-122.40869142607274, 37.76134617605315],
  [-122.41678443341827, 37.75828156322945],
  [-122.41692062210988, 37.75101772588085],
  [-122.4129352817101, 37.75267875552496],
  [-122.41153810982007, 37.75451800575868],
  [-122.4140060513352, 37.74657122463306],
  [-122.41385911135957, 37.76337653617459],
  [-122.42043665185116, 37.754407470692996],
  [-122.40978938157065, 37.76450303329318],
  [-122.41344777640442, 37.751148483269525],
  [-122.42287692579625, 37.76417974814643],
  [-122.41715126382405, 37.760537804214785],
  [-122.41004340162493, 37.74871935743222],
  [-122.41774596037905, 37.75231775479745],
  [-122.41401197674212, 37.751786210327225],
  [-122.4139005439404, 37.75392082167846],
  [-122.41743903596992, 37.76209659106257],
  [-122.40797490641594, 37.755656776209086],
  [-122.42601660166785, 37.76304044055552],
  [-122.42475869912616, 37.76128158121274],
  [-122.4254136565183, 37.76059258583783],
  [-122.40794615669083, 37.75104056209607],
  [-122.41075451027352, 37.747082929796356],
  [-122.41012514709995, 37.7612218618984],
  [-122.41975219244812, 37.763334928138605],
  [-122.4148893680859, 37.75175323687019],
  [-122.42381965154172, 37.76389286561216],
  [-122.42060974443976, 37.75539837180495],
  [-122.41486518168091, 37.76225919844404],
  [-122.41354296784603, 37.76001174749493],
  [-122.41250374000519, 37.76143653964776],
  [-122.42383388780732, 37.76418120866797],
  [-122.41948996087152, 37.76375791326102],
  [-122.40901934781371, 37.757190231652636],
  [-122.41446616923866, 37.7515671364096],
  [-122.41145761273533, 37.75843991693097],
  [-122.41993600982686, 37.764481854652736],
  [-122.4226566332339, 37.76406227103009],
  [-122.41246627252544, 37.76149261428097],
  [-122.42108957360912, 37.76062461888235],
  [-122.41586666626333, 37.75345065474162],
  [-122.42607748635427, 37.761078335088655],
  [-122.41783069338953, 37.750638766365256],
  [-122.40658575420164, 37.756449811909334],
  [-122.41111606218875, 37.75537378026224],
  [-122.42394850602186, 37.75855610670547],
  [-122.40867766046134, 37.75817465417819],
  [-122.41381831239967, 37.76440867305686],
  [-122.41757863525594, 37.74901705715819],
  [-122.41487544794536, 37.75488633070316],
  [-122.41854278041016, 37.764304983194315],
  [-122.42379568293042, 37.76250964248838],
  [-122.41468049513983, 37.747198302368524],
  [-122.40923776500622, 37.75333381902695],
  [-122.40822541947725, 37.76439953782074],
  [-122.41686980413054, 37.751790025694],
  [-122.40760247546336, 37.762676079668196],
  [-122.42196706998658, 37.760627863987054],
  [-122.41832341219747, 37.76337637245978],
  [-122.42603224084327, 37.763924686918585],
  [-122.41225380452875, 37.7483419775086],
  [-122.42240830574532, 37.75423945656585],
  [-122.42228960924437, 37.75435434305817],
  [-122.41936389765623, 37.759929842009704],
  [-122.40775071665323, 37.76530111041392],
  [-122.4099039982627, 37.761747153832005],
  [-122.41901248673905, 37.75719450024423],
  [-122.41573893502441, 37.763831603122455],
  [-122.41237555483644, 37.761252927558274],
  [-122.40455000117886, 37.75726219498988],
  [-122.41069965976091, 37.76470088472506],
  [-122.41879495724017, 37.7622151757651],
  [-122.41516185537776, 37.75725637834527],
  [-122.42005488320424, 37.76320058943922],
  [-122.40966406763532, 37.754856259950174],
  [-122.41267234354358, 37.75588417673643],
  [-122.41129594712292, 37.761662032164296],
  [-122.41501959094464, 37.7631304733426]]}