Creating Model Scenarios#
The initial conditions of the simulation are dictated by demographics (e.g., population, age distribution, etc.). The laser-measles package provides a number of tools to help you generate demographics for your simulation. These can be used for the abm, compartmental, and biweekly models.
In this tutorial, we’ll download and process a shapefile of Ethiopia at administrative level 1 boundaries to estimate intitial populations per patch. We will also show how we can sub-divide each boundary shape into roughly equal-area patches.
Setup and plot the shapefile#
laser-measles provides some functionality for downloading and plotting GADM shapefiles. Below we will download the data, print it as a dataframe, and then plot it. Note that we have constructed a DOTNAME attribute has the format COUNTRY:REGION. The data is located in the local directory.
[1]:
from pathlib import Path
from IPython.display import display
from laser.measles.demographics import GADMShapefile
from laser.measles.demographics import get_shapefile_dataframe
from laser.measles.demographics import plot_shapefile_dataframe
# Name of the shapefile you want to use
shapefile = Path("ETH/gadm41_ETH_1.shp")
# We will check whether it exists and download it
if not shapefile.exists():
shp = GADMShapefile.download("ETH", admin_level=1)
print("Shapefile is now at", shp.shapefile)
else:
print("Shapefile already exists")
shp = GADMShapefile(shapefile=shapefile, admin_level=1)
# Access the shapfile and metadata as a polars dataframe
# This looks like geopandas but is more limited.
df = get_shapefile_dataframe(shp.shapefile)
print(df.head(n=2))
# Plot the shapefile
plot_shapefile_dataframe(df, plot_kwargs={"facecolor": "xkcd:sky blue"})
Downloading GADM shapefile for Ethiopia |████████████████████████████████████████| 0 in 0.7s (0.00/s)
Shapefile is now at /home/docs/checkouts/readthedocs.org/user_builds/laser-measles/checkouts/stable/docs/tutorials/ETH/gadm41_ETH_1.shp
shape: (2, 13)
┌─────────┬───────┬──────────┬─────────────┬───┬────────┬───────┬────────────────┬────────────┐
│ GID_1 ┆ GID_0 ┆ COUNTRY ┆ NAME_1 ┆ … ┆ HASC_1 ┆ ISO_1 ┆ DOTNAME ┆ shape │
│ --- ┆ --- ┆ --- ┆ --- ┆ ┆ --- ┆ --- ┆ --- ┆ --- │
│ str ┆ str ┆ str ┆ str ┆ ┆ str ┆ str ┆ str ┆ object │
╞═════════╪═══════╪══════════╪═════════════╪═══╪════════╪═══════╪════════════════╪════════════╡
│ ETH.1_1 ┆ ETH ┆ Ethiopia ┆ Addis Abeba ┆ … ┆ ET.AA ┆ NA ┆ ethiopia:addis ┆ Polygon #0 │
│ ┆ ┆ ┆ ┆ ┆ ┆ ┆ abeba ┆ │
│ ETH.2_1 ┆ ETH ┆ Ethiopia ┆ Afar ┆ … ┆ ET.AF ┆ ET-AF ┆ ethiopia:afar ┆ Polygon #1 │
└─────────┴───────┴──────────┴─────────────┴───┴────────┴───────┴────────────────┴────────────┘
[1]:
Population calculation#
For the simulation we will want to know the initial number of people in each region. First we’ll download our population file (~5.6MB) from worldpop using standard libraries:
[2]:
import requests
url = "https://data.worldpop.org/GIS/Population/Global_2000_2020_1km_UNadj/2010/ETH/eth_ppp_2010_1km_Aggregated_UNadj.tif"
output_path = Path("ETH/eth_ppp_2010_1km_Aggregated_UNadj.tif")
if not output_path.exists():
response = requests.get(url, stream=True)
if response.status_code == 200:
with open(output_path, "wb") as f:
for chunk in response.iter_content(chunk_size=8192):
f.write(chunk)
print("Download complete.")
else:
print(f"Failed to download. Status code: {response.status_code}")
Download complete.
We use the RasterPatchGenerator to sum the population in each of the shapes. This is saved into a dataframe that we can use to initialize a simulation.
[3]:
import sciris as sc
from laser.measles.demographics import RasterPatchGenerator
from laser.measles.demographics import RasterPatchParams
# Setup demographics generator
config = RasterPatchParams(
id="ETH_ADM1",
region="ETH",
shapefile=shp.shapefile,
population_raster=output_path,
)
# Create the generator
generator = RasterPatchGenerator(config)
# Time the population calculation
with sc.Timer() as t:
# Generate the demographics (in this case the population per patch)
generator.generate_demographics()
print(f"Total population: {generator.population['pop'].sum() / 1e6:.2f} million") # Should be ~90.5M
# the result is stored in a polars dataframe and can be accessed via `population`
generator.population.head(n=2)
on 0: Loading data...
on 0: [Field(name="DeletionFlag", field_type=FieldType.C, size=1, decimal=0), Field(name="GID_1", field_type=FieldType.C, size=10, decimal=0), Field(name="GID_0", field_type=FieldType.C, size=10, decimal=0), Field(name="COUNTRY", field_type=FieldType.C, size=10, decimal=0), Field(name="NAME_1", field_type=FieldType.C, size=31, decimal=0), Field(name="VARNAME_1", field_type=FieldType.C, size=35, decimal=0), Field(name="NL_NAME_1", field_type=FieldType.C, size=10, decimal=0), Field(name="TYPE_1", field_type=FieldType.C, size=10, decimal=0), Field(name="ENGTYPE_1", field_type=FieldType.C, size=10, decimal=0), Field(name="CC_1", field_type=FieldType.C, size=10, decimal=0), Field(name="HASC_1", field_type=FieldType.C, size=10, decimal=0), Field(name="ISO_1", field_type=FieldType.C, size=10, decimal=0), Field(name="DOTNAME", field_type=FieldType.C, size=50, decimal=0)]
on 0: Clipping:
on 0: 1 of 11 (9%) ethiopia:addis abeba (np.float64(38.785538550534746), np.float64(8.98048294246869)) {'lat': np.float64(8.98048294246869), 'lon': np.float64(38.785538550534746), 'pop': 3189938}
on 0: 2 of 11 (18%) ethiopia:afar (np.float64(40.767449344797555), np.float64(12.037909353819343)) {'lat': np.float64(12.037909353819343), 'lon': np.float64(40.767449344797555), 'pop': 1619009}
on 0: 3 of 11 (27%) ethiopia:amhara (np.float64(38.04732030224195), np.float64(11.562494255356375)) {'lat': np.float64(11.562494255356375), 'lon': np.float64(38.04732030224195), 'pop': 20078508}
on 0: 4 of 11 (36%) ethiopia:benshangul-gumaz (np.float64(35.44237778797806), np.float64(10.50800671019417)) {'lat': np.float64(10.50800671019417), 'lon': np.float64(35.44237778797806), 'pop': 942015}
on 0: 5 of 11 (45%) ethiopia:dire dawa (np.float64(42.00302663044891), np.float64(9.60626902481112)) {'lat': np.float64(9.60626902481112), 'lon': np.float64(42.00302663044891), 'pop': 403830}
on 0: 6 of 11 (55%) ethiopia:gambela peoples (np.float64(34.328623182923636), np.float64(7.7073679551810645)) {'lat': np.float64(7.7073679551810645), 'lon': np.float64(34.328623182923636), 'pop': 391565}
on 0: 7 of 11 (64%) ethiopia:harari people (np.float64(42.172525982646384), np.float64(9.289660181001201)) {'lat': np.float64(9.289660181001201), 'lon': np.float64(42.172525982646384), 'pop': 219787}
on 0: 8 of 11 (73%) ethiopia:oromia (np.float64(38.76544158459778), np.float64(7.507237452966563)) {'lat': np.float64(7.507237452966563), 'lon': np.float64(38.76544158459778), 'pop': 32432831}
on 0: 9 of 11 (82%) ethiopia:somali (np.float64(43.33707510619852), np.float64(6.913614872025026)) {'lat': np.float64(6.913614872025026), 'lon': np.float64(43.33707510619852), 'pop': 5284320}
on 0: 10 of 11 (91%) ethiopia:southern nations, nationalities (np.float64(36.811052337330274), np.float64(6.466449317786186)) {'lat': np.float64(6.466449317786186), 'lon': np.float64(36.811052337330274), 'pop': 17839278}
on 0: 11 of 11 (100%) ethiopia:tigray (np.float64(38.441497467078264), np.float64(13.776739407422726)) {'lat': np.float64(13.776739407422726), 'lon': np.float64(38.441497467078264), 'pop': 5087591}
Clipping population raster to shapefile |████████████████████████████████████████| 0 in 6.2s (0.00/s)
Total population: 87.49 million
Elapsed time: 6.25 s
[3]:
shape: (2, 4)
| dotname | lat | lon | pop |
|---|---|---|---|
| str | f64 | f64 | i64 |
| "ethiopia:addis abeba" | 8.980483 | 38.785539 | 3189938 |
| "ethiopia:afar" | 12.037909 | 40.767449 | 1619009 |
laser-measles demographics uses caching to save results. Now we will run the calculation again with a new instance of the RasterPatchGenerator.
[4]:
new_generator = RasterPatchGenerator(config)
with sc.Timer() as t:
# # Generate the demographics (in this case the population)
new_generator.generate_demographics()
print(f"Total population: {new_generator.population['pop'].sum() / 1e6:.2f} million") # Should be ~90.5M
# Note how the time to run the `generate_demographics` method a second time is greatly improved.
Total population: 87.49 million
Elapsed time: 17.6 ms
You can access the cache directory using the associated module
[5]:
from laser.measles.demographics import cache
print(f"Cache directory: {cache.get_cache_dir()}")
Cache directory: /home/docs/.cache/laser.measles
Sub-divide the regions#
Now we will generate roughtly equal area patches of 700 km using the original shp shapefile. Now each shape has a unique identifier with the form COUNTRY:REGION:ID. We will also time how long this takes.
[6]:
# Set the patch size
patch_size = 700 # sq km
# Create the GADMShapefile using the original shapefile
new_shp = GADMShapefile(shapefile=shp.shapefile, admin_level=1)
# Subdivide the original shapefile (this is costly)
new_shp.shape_subdivide(patch_size_km=patch_size)
print("Shapefile is now at", new_shp.shapefile)
# Get the results as a polars dataframe
new_df = get_shapefile_dataframe(new_shp.shapefile)
display(new_df.head(n=2))
# Plot the resulting shapes
import matplotlib.pyplot as plt
plt.figure()
ax = plt.gca()
plot_shapefile_dataframe(new_df, plot_kwargs={"facecolor": "xkcd:sky blue", "edgecolor": "gray"}, ax=ax)
plot_shapefile_dataframe(df, plot_kwargs={"fill": False}, ax=ax)
on 0: [Field(name="DeletionFlag", field_type=FieldType.C, size=1, decimal=0), Field(name="GID_1", field_type=FieldType.C, size=10, decimal=0), Field(name="GID_0", field_type=FieldType.C, size=10, decimal=0), Field(name="COUNTRY", field_type=FieldType.C, size=10, decimal=0), Field(name="NAME_1", field_type=FieldType.C, size=31, decimal=0), Field(name="VARNAME_1", field_type=FieldType.C, size=35, decimal=0), Field(name="NL_NAME_1", field_type=FieldType.C, size=10, decimal=0), Field(name="TYPE_1", field_type=FieldType.C, size=10, decimal=0), Field(name="ENGTYPE_1", field_type=FieldType.C, size=10, decimal=0), Field(name="CC_1", field_type=FieldType.C, size=10, decimal=0), Field(name="HASC_1", field_type=FieldType.C, size=10, decimal=0), Field(name="ISO_1", field_type=FieldType.C, size=10, decimal=0), Field(name="DOTNAME", field_type=FieldType.C, size=50, decimal=0)]
Subdividing shapefile gadm41_ETH_1 |████████████████████████████████████████| 0 in 0.1s (0.00/s)
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[6], line 8
5 new_shp = GADMShapefile(shapefile=shp.shapefile, admin_level=1)
7 # Subdivide the original shapefile (this is costly)
----> 8 new_shp.shape_subdivide(patch_size_km=patch_size)
9 print("Shapefile is now at", new_shp.shapefile)
11 # Get the results as a polars dataframe
File ~/checkouts/readthedocs.org/user_builds/laser-measles/envs/stable/lib/python3.12/site-packages/laser/measles/demographics/admin_shapefile.py:52, in AdminShapefile.shape_subdivide(self, patch_size_km)
48 self.add_dotname()
49 with alive_progress.alive_bar(
50 title=f"Subdividing shapefile {self.shapefile.stem}",
51 ) as _:
---> 52 shape_subdivide(
53 shape_stem=self.shapefile,
54 out_dir=self.get_shapefile_parent(),
55 out_suffix=f"{patch_size_km}km",
56 box_target_area_km2=patch_size_km,
57 )
58 self.shapefile = out_file
59 return out_file
File ~/checkouts/readthedocs.org/user_builds/laser-measles/envs/stable/lib/python3.12/site-packages/rastertoolkit/shape.py:526, in shape_subdivide(shape_stem, out_dir, out_suffix, output_centers, top_n, shape_attr, box_target_area_km2, points_per_box, random_seed, verbose)
524 # Add the new shape to the shapefile; splat the record
525 sf1new.poly(poly_as_list)
--> 526 sf1new.record(**new_recs)
528 if output_centers and new_recs is not None:
529 for i, p in enumerate([Point(xy) for xy in sub_node]):
File ~/checkouts/readthedocs.org/user_builds/laser-measles/envs/stable/lib/python3.12/site-packages/shapefile.py:3720, in Writer.record(self, *recordList, **recordDict)
3717 else:
3718 # Blank fields for empty record
3719 record = ["" for _ in range(fieldCount)]
-> 3720 self.__dbfRecord(record)
File ~/checkouts/readthedocs.org/user_builds/laser-measles/envs/stable/lib/python3.12/site-packages/shapefile.py:3729, in Writer.__dbfRecord(self, record)
3724 f = self.__getFileObj(self.dbf)
3725 if self.recNum == 0:
3726 # first records, so all fields should be set
3727 # allowing us to write the dbf header
3728 # cannot change the fields after this point
-> 3729 self.__dbfHeader()
3730 # first byte of the record is deletion flag, always disabled
3731 f.write(b" ")
File ~/checkouts/readthedocs.org/user_builds/laser-measles/envs/stable/lib/python3.12/site-packages/shapefile.py:3555, in Writer.__dbfHeader(self)
3551 if headerLength >= 65535:
3552 raise ShapefileException(
3553 "Shapefile dbf header length exceeds maximum length."
3554 )
-> 3555 recordLength = sum(field.size for field in fields) + 1
3556 header = pack(
3557 "<BBBBLHH20x",
3558 version,
(...) 3564 recordLength,
3565 )
3566 f.write(header)
File ~/checkouts/readthedocs.org/user_builds/laser-measles/envs/stable/lib/python3.12/site-packages/shapefile.py:3555, in <genexpr>(.0)
3551 if headerLength >= 65535:
3552 raise ShapefileException(
3553 "Shapefile dbf header length exceeds maximum length."
3554 )
-> 3555 recordLength = sum(field.size for field in fields) + 1
3556 header = pack(
3557 "<BBBBLHH20x",
3558 version,
(...) 3564 recordLength,
3565 )
3566 f.write(header)
AttributeError: 'tuple' object has no attribute 'size'