%load_ext autoreload
%autoreload 2

import warnings
from functools import partial

from pisa.administrative_area import AdministrativeArea
from pisa.facilities import Facilities
from pisa.population import WorldpopPopulation
from pisa.population_served_by_isopolygons import get_population_served_by_isopolygons
from pisa.visualisation import (
    plot_facilities,
    plot_isochrones,
    plot_population,
    plot_population_heatmap,
)

from optimization import jg_opt

INFO:matplotlib.font_manager:Failed to extract font properties from /usr/share/fonts/truetype/noto/NotoColorEmoji.ttf: Can not load face (unknown file format; error code 0x2)

INFO:matplotlib.font_manager:generated new fontManager

import os

from dotenv import load_dotenv

# load environment variables from an `.env` file in the local root directory
load_dotenv()

CBC_SOLVER_PATH = os.getenv("CBC_SOLVER_PATH")  # path to the cbc executable (e.g. /opt/homebrew/bin/cbc)

Pisa Showcase, Example Notebook

Use Case Example notebook

This notebook shows how to use the PISA package to calculate new locations for hospitals where it would have most impact, i.e. reach the highest possible number of people that currently do not have a hospital in their vicinity.

The example is based on the case of Baucau, Timor-Leste, and we are interested in the part of the population that cannot reach a hospital when driving 10 km. WorldPop is used as source to get the population count in Baucau area and the road network from OSM is used to identify their distance to current hospitals and potential locations for new facilities. Depending on the specified budget, the optimization will return a number of locations where new hospitals would reach the highest number of people that currently do not have a hospital in their vicinity.

Define Administrative Area

We need to define the administrative area that we want to work with. This is done in two steps:

1. Create an AdministrativeArea object with the country name and administrative level of interest.

2. Call the get_admin_area_boundaries method on this AdministrativeArea object with the name of the region of interest to get the boundaries of the administrative area. The regions that can be specified here are dependent on which administrative level was specified in the first step.

timor_leste = AdministrativeArea(country_name="Timor-Leste", admin_level=1)

# these are the boundaries of Baucau
# type: Polygon
baucau = timor_leste.get_admin_area_boundaries("Baucau")

INFO:pisa.administrative_area:Validating country name: Timor-Leste

INFO:pisa.administrative_area:Country name 'Timor-Leste' validated successfully

INFO:pisa.administrative_area:Retrieving boundaries of all administrative areas of level 1 for country Timor-Leste

  0%|          | 0.00/6.25M [00:00<?, ?B/s]

  1%|          | 48.1k/6.25M [00:00<00:15, 411kB/s]

  3%|▎         | 201k/6.25M [00:00<00:06, 917kB/s] 

 14%|█▎        | 849k/6.25M [00:00<00:01, 2.99MB/s]

 55%|█████▌    | 3.45M/6.25M [00:00<00:00, 10.4MB/s]

 96%|█████████▋| 6.02M/6.25M [00:00<00:00, 15.3MB/s]

100%|██████████| 6.25M/6.25M [00:00<00:00, 10.7MB/s]

Facilities

Get existing facilities (in our case, hospitals) from OSM

We first identify the existing hospitals in Baucau by calling the get_existing_facilities method on the Facilities class, which takes the administrative area boundaries as input. This method fetches the existing facilities from OpenStreetMap (OSM).

# Suppress user warning about geometry in geographic CRS. Centroid is calculated over a single facility (e.g. a hospital),
# so distances are very small and projection isn't necessary
warnings.filterwarnings(
    "ignore",
    message="Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect",
    category=UserWarning,
)

hospitals_df = Facilities(admin_area_boundaries=baucau).get_existing_facilities()

INFO:pisa.facilities:Retrieving existing facilities with tags {'amenity': 'hospital', 'building': 'hospital'} using OSM.

INFO:pisa.facilities:Successfully retrieved existing facilities from OSM.

plot_facilities(hospitals_df, baucau)

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

Estimate potential locations for new facilities

We then identify potential locations for new hospitals in Baucau by drawing a grid over the area. The spacing parameter specifies the distance between potential facilities. In this case, we use a spacing of 0.05.

potential_hospitals_df = Facilities(admin_area_boundaries=baucau).estimate_potential_facilities(spacing=0.05)

Get population

We get the population data for Baucau from WorldPop by calling the get_population_gdf method on the WorldpopPopulation class, which takes the administrative area boundaries and the ISO3 country code as input parameters. The ISO3 country code is obtained from the AdministrativeArea object.

WorldPop always provides the most recent population data available, which is usually from the previous year. There is also a possibility to use Facebook population data in the PISA package, but the most recent data available is from 2019. In this example, we use WorldPop.

population_gdf = WorldpopPopulation(
    admin_area_boundaries=baucau, iso3_country_code=timor_leste.get_iso3_country_code()
).get_population_gdf()

population_gdf.head()

	longitude	latitude	population	geometry
0	126.14917	-8.67333	0.10	POINT (126.14917 -8.67333)
1	126.14917	-8.67250	0.09	POINT (126.14917 -8.6725)
2	126.15000	-8.67333	0.10	POINT (126.15 -8.67333)
3	126.15000	-8.67250	0.09	POINT (126.15 -8.6725)
4	126.15000	-8.67167	0.09	POINT (126.15 -8.67167)

plot_population_heatmap(population_gdf, baucau)

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

plot_population(population_gdf, baucau, random_sample_n=1000)

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

Calculate isopolygons

We specify the parameters for determining access to existing and potential facilities.

For this we make choices on:

• distance type: in this notebook we use length (i.e. distance in meters) as the distance type, but it could also be time (i.e. travel time in minutes)

• distance values: these are the distances we want to use to determine access to hospitals. In this example, we use 2000, 5000 and 10000 meters (i.e. 2, 5 and 10 km) as the distance values.

• mode of transport: in this example, we use driving as the mode of transport, but it could also be walking or cycling.

Valid values for constants can be found in the script pisa.constants

DISTANCE_TYPE = "length"

DISTANCE_VALUES = [2000, 5000, 10000]

MODE_OF_TRANSPORT = "driving"

Using OSM

In order to determine the population that lives 10 km or less driving from an existing or potential facility, we need to get the road network from OSM. This can be done by calling the get_osm_road_network method on the OsmRoadNetwork class, which takes the administrative area boundaries, mode of transport, and distance type as input parameters. This will return the road network for the specific input parameters that can be used in the following step.

from pisa.osm_road_network import OsmRoadNetwork

road_network = OsmRoadNetwork(
    admin_area_boundaries=baucau,
    mode_of_transport=MODE_OF_TRANSPORT,
    distance_type=DISTANCE_TYPE,
).get_osm_road_network()

INFO:pisa.osm_road_network:OSM road network set with parameters network_type 'drive' and distance_type 'length'

Calculate isopolygons for existing facilities

We calculate the areas (isopolygons) around the existing facilities that are reachable within the specified distances and transport mode. This is done by using the calculate_isopolygons method of the OsmIsopolygonCalculator class, which takes the existing facilities, distance type, distance values, and road network as input parameters.

from pisa.isopolygons import OsmIsopolygonCalculator

with warnings.catch_warnings():
    warnings.simplefilter("ignore", category=UserWarning)
    isopolygons_existing_facilities = OsmIsopolygonCalculator(
        facilities_df=hospitals_df,
        distance_type=DISTANCE_TYPE,
        distance_values=DISTANCE_VALUES,
        road_network=road_network,
    ).calculate_isopolygons()

isopolygons_existing_facilities.head()

	ID_2000	ID_5000	ID_10000
1777614876	POLYGON ((126.59878421714033 -8.54764218472667...	POLYGON ((126.59761276961231 -8.57116537023302...	POLYGON ((126.56068305448515 -8.56751807760911...
1777614896	POLYGON ((126.65907104033573 -8.62868878467725...	POLYGON ((126.6311023795325 -8.646327183432366...	POLYGON ((126.63106042126435 -8.64641589673682...
1777615024	POLYGON ((126.59499222126435 -8.47806089673682...	POLYGON ((126.55950911149382 -8.48560952854904...	POLYGON ((126.57412141532275 -8.51524054033573...
1777615035	POLYGON ((126.68354761081066 -8.63974323564052...	POLYGON ((126.68357328721343 -8.63970069138097...	POLYGON ((126.6311023795325 -8.646327183432366...
1777615036	POLYGON ((126.6803646 -8.578027, 126.680359784...	POLYGON ((126.6650524795325 -8.602129783432366...	POLYGON ((126.65894556961231 -8.62895407023302...

plot_isochrones(isopolygons_existing_facilities, baucau)

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

Calculate isopolygons for potential facilities

We do the same for the potential facilities. This will give us the areas around the potential facilities that are reachable within the specified distances and transport mode.

with warnings.catch_warnings():
    warnings.simplefilter("ignore", category=UserWarning)
    isopolygons_potential_facilities = OsmIsopolygonCalculator(
        facilities_df=potential_hospitals_df,
        distance_type=DISTANCE_TYPE,
        distance_values=DISTANCE_VALUES,
        road_network=road_network,
    ).calculate_isopolygons()

Get population served by isopolygons

We can now calculate the population that lives within the isopolygons for both existing and potential facilities. This is done by using the get_population_served_by_isopolygons function, which takes the population data and the isopolygons as input parameters. We save this value in a dictionary with the distance type as key.

population_served_current = get_population_served_by_isopolygons(
    grouped_population=population_gdf, isopolygons=isopolygons_existing_facilities
)

current = {DISTANCE_TYPE: population_served_current}

population_served_potential = get_population_served_by_isopolygons(
    grouped_population=population_gdf, isopolygons=isopolygons_potential_facilities
)

potential = {DISTANCE_TYPE: population_served_potential}

Prepare optimization data

Preparing the variables that go into the optimization:

1. We define the population count, which is the total population in the area of interest. This is used to calculate the number of people that can be served by the existing and potential facilities.

2. We define the budget for the optimization, which is the number of locations that can be built. In this example, we use a budget of 5, 20 and 50 locations.

3. We set the optimization method to use the CBC solver. The jg_opt.OpenOptimize function is used to create an optimization object with the specified solver path (defined in a .env file).

population_count = population_gdf.population.values

BUDGET = [
    5,
    20,
    50,
]  # budget for the optimization in terms of how many locations can be built

cbc_optimize = partial(jg_opt.OpenOptimize, solver_path=CBC_SOLVER_PATH)

Optimization

We run the optimization using the jg_opt.Solve function, which takes the population count, current population served by existing facilities, population that could be served by the potential facilities, distance type, budget, and optimization method as input parameters. The optimization will return the values (percentage of population that can be reached with the new facilities) and solutions (the best locations for new facilities) for the specified budgets.

values, solutions = jg_opt.Solve(
    household=population_count,
    current=current,
    potential=potential,
    accessibility=DISTANCE_TYPE,
    budgets=BUDGET,
    optimize=cbc_optimize,
    type="ID",
)

WARNING:pyomo.opt:Failed to create solver with name '_cbc_shell':
Failed to set executable for solver cbc. File with name=D:\EiriniK\Downloads\amplbundle.mswin64\ampl.mswin64\cbc.exe either does not exist or it is not executable. To skip this validation, call set_executable with validate=False.
Traceback (most recent call last):
  File "/home/runner/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/pyomo/opt/base/solvers.py", line 148, in __call__
    opt = self._cls[_name](**kwds)
  File "/home/runner/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/pyomo/solvers/plugins/solvers/CBCplugin.py", line 88, in __init__
    super(CBCSHELL, self).__init__(**kwds)
  File "/home/runner/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/pyomo/opt/solver/shellcmd.py", line 66, in __init__
    self.set_executable(name=executable, validate=validate)
  File "/home/runner/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/pyomo/opt/solver/shellcmd.py", line 115, in set_executable
    raise ValueError(
ValueError: Failed to set executable for solver cbc. File with name=D:\EiriniK\Downloads\amplbundle.mswin64\ampl.mswin64\cbc.exe either does not exist or it is not executable. To skip this validation, call set_executable with validate=False.

WARNING:pyomo.opt:Failed to create solver with name 'cbc':
Attempting to use an unavailable solver.

The SolverFactory was unable to create the solver "_cbc_shell"
and returned an UnknownSolver object.  This error is raised at the point
where the UnknownSolver object was used as if it were valid (by calling
method "set_problem_format").

The original solver was created with the following parameters:
	executable: D:\EiriniK\Downloads\amplbundle.mswin64\ampl.mswin64\cbc.exe
	type: _cbc_shell
	_args: ()
	options: {}
Traceback (most recent call last):
  File "/home/runner/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/pyomo/opt/base/solvers.py", line 148, in __call__
    opt = self._cls[_name](**kwds)
  File "/home/runner/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/pyomo/solvers/plugins/solvers/CBCplugin.py", line 50, in __new__
    opt.set_problem_format(ProblemFormat.cpxlp)
  File "/home/runner/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/pyomo/opt/base/solvers.py", line 113, in __getattr__
    self._solver_error(attr)
  File "/home/runner/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/pyomo/opt/base/solvers.py", line 116, in _solver_error
    raise RuntimeError(
RuntimeError: Attempting to use an unavailable solver.

The SolverFactory was unable to create the solver "_cbc_shell"
and returned an UnknownSolver object.  This error is raised at the point
where the UnknownSolver object was used as if it were valid (by calling
method "set_problem_format").

The original solver was created with the following parameters:
	executable: D:\EiriniK\Downloads\amplbundle.mswin64\ampl.mswin64\cbc.exe
	type: _cbc_shell
	_args: ()
	options: {}

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[22], line 1
----> 1 values, solutions = jg_opt.Solve(
      2     household=population_count,
      3     current=current,
      4     potential=potential,
      5     accessibility=DISTANCE_TYPE,
      6     budgets=BUDGET,
      7     optimize=cbc_optimize,
      8     type="ID",
      9 )

File ~/work/Public-Infrastructure-Service-Access/Public-Infrastructure-Service-Access/optimization/jg_opt.py:362, in Solve(household, current, potential, accessibility, budgets, optimize, type)
    360 I = np.unique(list(IJ.keys()))
    361 J = np.unique(np.concatenate(list(IJ.values())))
--> 362 optimization = optimize(
    363     household,
    364     I,
    365     J,
    366     IJ,
    367     [max(budgets)],
    368     parsimonious=True,
    369     maxTimeInSeconds=5,
    370     mipGap=1e-15,
    371 )
    372 optimization["nof"] = [len(s) for s in optimization.solution]
    373 coverage = (optimization.value / household.sum() + percent_covered).to_frame()

File ~/work/Public-Infrastructure-Service-Access/Public-Infrastructure-Service-Access/optimization/jg_opt.py:126, in OpenOptimize(w, I, J, IJ, budget_list, parsimonious, maxTimeInSeconds, mipGap, trace, solver, solver_path)
    124 result.at[p, "modeling"] = pc() - start
    125 start = pc()
--> 126 solver_result = solver.solve(M, tee=trace)
    127 result.at[p, "solving"] = pc() - start
    128 result.at[p, "value"] = int(
    129     np.ceil(M.weighted_coverage() - np.finfo(np.float16).eps)
    130 )

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/pyomo/opt/base/solvers.py:99, in UnknownSolver.solve(self, *args, **kwds)
     97 def solve(self, *args, **kwds):
     98     """Perform optimization and return an SolverResults object."""
---> 99     self._solver_error('solve')

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/pyomo/opt/base/solvers.py:116, in UnknownSolver._solver_error(self, method_name)
    115     def _solver_error(self, method_name):
--> 116         raise RuntimeError(
    117             """Attempting to use an unavailable solver.
    118 
    119 The SolverFactory was unable to create the solver "%s"
    120 and returned an UnknownSolver object.  This error is raised at the point
    121 where the UnknownSolver object was used as if it were valid (by calling
    122 method "%s").
    123 
    124 The original solver was created with the following parameters:
    125 \t"""
    126             % (self.type, method_name)
    127             + "\n\t".join("%s: %s" % i for i in sorted(self._kwds.items()))
    128             + "\n\t_args: %s" % (self._args,)
    129             + "\n\toptions: %s" % (self.options,)
    130         )

RuntimeError: Attempting to use an unavailable solver.

The SolverFactory was unable to create the solver "cbc"
and returned an UnknownSolver object.  This error is raised at the point
where the UnknownSolver object was used as if it were valid (by calling
method "solve").

The original solver was created with the following parameters:
	executable: D:\EiriniK\Downloads\amplbundle.mswin64\ampl.mswin64\cbc.exe
	type: cbc
	_args: ()
	options: {'threads': 8, 'sec': 5.0, 'allowableGap': 1e-15}

values

solutions