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

Create an AdministrativeArea object with the country name and administrative level of interest.
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%|          | 32.8k/6.25M [00:00<00:24, 256kB/s]

  2%|▏         | 144k/6.25M [00:00<00:10, 600kB/s]

  6%|▌         | 377k/6.25M [00:00<00:04, 1.23MB/s]

  9%|▉         | 560k/6.25M [00:00<00:04, 1.36MB/s]

 17%|█▋        | 1.07M/6.25M [00:00<00:02, 2.50MB/s]

 26%|██▌       | 1.64M/6.25M [00:00<00:01, 3.32MB/s]

 45%|████▍     | 2.79M/6.25M [00:00<00:00, 5.60MB/s]

 54%|█████▍    | 3.37M/6.25M [00:00<00:00, 5.36MB/s]

 92%|█████████▏| 5.73M/6.25M [00:01<00:00, 10.3MB/s]

100%|██████████| 6.25M/6.25M [00:01<00:00, 5.72MB/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.

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[5], line 9
# 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,
)
----> 9 hospitals_df = Facilities(admin_area_boundaries=baucau).get_existing_facilities()

File ~/work/Public-Infrastructure-Service-Access/Public-Infrastructure-Service-Access/pisa/facilities.py:135, in Facilities.get_existing_facilities(self)
"""Retrieve existing facilities from the specified data source.

This method acts as a dispatcher that calls the appropriate data source-specific method based on the `data_src`
   (...)
implement new methods and update this dispatcher method accordingly, or refactor to use an abstract base class.
"""
if self.data_src == "osm":
--> 135     return self._get_existing_facilities_osm(
       admin_area_boundaries=self.admin_area_boundaries, osm_tags=self.osm_tags
   )
raise NotImplementedError(f"Data source '{self.data_src}' not implemented")

File ~/work/Public-Infrastructure-Service-Access/Public-Infrastructure-Service-Access/pisa/facilities.py:169, in Facilities._get_existing_facilities_osm(osm_tags, admin_area_boundaries)
# retrieves facilities GeodataFrame from osm
try:
--> 169     facilities_gdf = ox.features_from_polygon(polygon=admin_area_boundaries, tags=osm_tags)
except InsufficientResponseError:
   facilities_gdf = gpd.GeoDataFrame(
       pd.DataFrame(columns=["id", "element", "amenity", "geometry"]),
       geometry=[],
       crs="EPSG:4326",
   )

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/osmnx/features.py:332, in features_from_polygon(polygon, tags)
# retrieve the data from Overpass then turn it into a GeoDataFrame
response_jsons = _overpass._download_overpass_features(polygon, tags)
--> 332 return _create_gdf(response_jsons, polygon, tags)

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/osmnx/features.py:406, in _create_gdf(response_jsons, polygon, tags)
elements = []
response_count = 0
--> 406 for response_json in response_jsons:
   response_count += 1
   if not settings.cache_only_mode:

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/osmnx/_overpass.py:427, in _download_overpass_features(polygon, tags)
for polygon_coord_str in polygon_coord_strs:
   query_str = _create_overpass_features_query(polygon_coord_str, tags)
--> 427     yield _overpass_request(OrderedDict(data=query_str))

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/osmnx/_overpass.py:466, in _overpass_request(data, pause, error_pause)
# pause then request this URL
if pause is None:
--> 466     this_pause = _get_overpass_pause(settings.overpass_url)
hostname = _http._hostname_from_url(url)
msg = f"Pausing {this_pause} second(s) before making HTTP POST request to {hostname!r}"

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/osmnx/_overpass.py:177, in _get_overpass_pause(base_endpoint, recursive_delay, default_duration)
# try to retrieve the URL
try:
--> 177     response = requests.get(
       url,
       headers=_http._get_http_headers(),
       timeout=settings.requests_timeout,
       **settings.requests_kwargs,
   )
   response_text = response.text
except RequestsConnectionError as e:  # pragma: no cover
   # cannot reach status endpoint: log error and return default duration

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/requests/api.py:73, in get(url, params, **kwargs)
def get(url, params=None, **kwargs):
   r"""Sends a GET request.

   :param url: URL for the new :class:`Request` object.
   (...)
   :rtype: requests.Response
   """
---> 73     return request("get", url, params=params, **kwargs)

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/requests/api.py:59, in request(method, url, **kwargs)
# By using the 'with' statement we are sure the session is closed, thus we
# avoid leaving sockets open which can trigger a ResourceWarning in some
# cases, and look like a memory leak in others.
with sessions.Session() as session:
---> 59     return session.request(method=method, url=url, **kwargs)

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/requests/sessions.py:589, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
send_kwargs = {
   "timeout": timeout,
   "allow_redirects": allow_redirects,
}
send_kwargs.update(settings)
--> 589 resp = self.send(prep, **send_kwargs)
return resp

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/requests/sessions.py:703, in Session.send(self, request, **kwargs)
start = preferred_clock()
# Send the request
--> 703 r = adapter.send(request, **kwargs)
# Total elapsed time of the request (approximately)
elapsed = preferred_clock() - start

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/requests/adapters.py:667, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
   timeout = TimeoutSauce(connect=timeout, read=timeout)
try:
--> 667     resp = conn.urlopen(
       method=request.method,
       url=url,
       body=request.body,
       headers=request.headers,
       redirect=False,
       assert_same_host=False,
       preload_content=False,
       decode_content=False,
       retries=self.max_retries,
       timeout=timeout,
       chunked=chunked,
   )
except (ProtocolError, OSError) as err:
   raise ConnectionError(err, request=request)

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/urllib3/connectionpool.py:787, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
response_conn = conn if not release_conn else None
# Make the request on the HTTPConnection object
--> 787 response = self._make_request(
   conn,
   method,
   url,
   timeout=timeout_obj,
   body=body,
   headers=headers,
   chunked=chunked,
   retries=retries,
   response_conn=response_conn,
   preload_content=preload_content,
   decode_content=decode_content,
   **response_kw,
)
# Everything went great!
clean_exit = True

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/urllib3/connectionpool.py:534, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
# Receive the response from the server
try:
--> 534     response = conn.getresponse()
except (BaseSSLError, OSError) as e:
   self._raise_timeout(err=e, url=url, timeout_value=read_timeout)

File ~/.cache/pypoetry/virtualenvs/pisa-abw-OIyr1Vno-py3.10/lib/python3.10/site-packages/urllib3/connection.py:516, in HTTPConnection.getresponse(self)
_shutdown = getattr(self.sock, "shutdown", None)
# Get the response from http.client.HTTPConnection
--> 516 httplib_response = super().getresponse()
try:
   assert_header_parsing(httplib_response.msg)

File /opt/hostedtoolcache/Python/3.10.22/x64/lib/python3.10/http/client.py:1423, in HTTPConnection.getresponse(self)
try:
   try:
-> 1423         response.begin()
   except ConnectionError:
       self.close()

File /opt/hostedtoolcache/Python/3.10.22/x64/lib/python3.10/http/client.py:330, in HTTPResponse.begin(self)
# read until we get a non-100 response
for _ in range(_MAXINTERIMRESPONSES):
--> 330     version, status, reason = self._read_status()
   if status != CONTINUE:
       break

File /opt/hostedtoolcache/Python/3.10.22/x64/lib/python3.10/http/client.py:291, in HTTPResponse._read_status(self)
def _read_status(self):
--> 291     line = str(self.fp.readline(_MAXLINE + 1), "iso-8859-1")
   if len(line) > _MAXLINE:
       raise LineTooLong("status line")

File /opt/hostedtoolcache/Python/3.10.22/x64/lib/python3.10/socket.py:717, in SocketIO.readinto(self, b)
while True:
   try:
--> 717         return self._sock.recv_into(b)
   except timeout:
       self._timeout_occurred = True

File /opt/hostedtoolcache/Python/3.10.22/x64/lib/python3.10/ssl.py:1320, in SSLSocket.recv_into(self, buffer, nbytes, flags)
   if flags != 0:
       raise ValueError(
         "non-zero flags not allowed in calls to recv_into() on %s" %
         self.__class__)
-> 1320     return self.read(nbytes, buffer)
else:
   return super().recv_into(buffer, nbytes, flags)

File /opt/hostedtoolcache/Python/3.10.22/x64/lib/python3.10/ssl.py:1176, in SSLSocket.read(self, len, buffer)
try:
   if buffer is not None:
-> 1176         return self._sslobj.read(len, buffer)
   else:
       return self._sslobj.read(len)

KeyboardInterrupt: 

plot_facilities(hospitals_df, baucau)

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

plot_population_heatmap(population_gdf, baucau)

plot_population(population_gdf, baucau, random_sample_n=1000)

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

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

plot_isochrones(isopolygons_existing_facilities, baucau)

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:

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.
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.
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",
)

values

solutions