3 — Score locations end to end¶
Build a scoring model out of the catalog data, rate two candidate locations against it, then score a whole area as a grid.
The model is deliberately small — three factors over two layers — because the point is the mechanics, not the model. Petala's flagship case is a model trained on your outcomes; this is the hand-weighted version you start from.
1. Point the client at an instance¶
These are the only two lines you change. PETALA_URL is the instance root;
PETALA_API_KEY is an organization key (ptk_...), minted by an org owner or
in the app (account menu → API keys), or with the admin token at
POST /api/orgs/{org_id}/api-keys. Both values can also come from the environment, which is how this
notebook was executed for the committed outputs.
An empty key only works against a dev-mode instance (one started with no
PETALA_ADMIN_TOKEN). The key is never printed here, and must never be pasted
into a notebook you commit.
import os
PETALA_URL = os.environ.get("PETALA_URL", "http://localhost:8020")
PETALA_API_KEY = os.environ.get("PETALA_API_KEY", "")
from petala import Petala
p = Petala(PETALA_URL, api_key=PETALA_API_KEY or None)
meta = p.meta()
print("connected to", p.base_url)
print("server ", meta["name"], meta["version"], "| auth_required:", meta["auth_required"])
print("api key ", "set (ptk_..., not printed)" if PETALA_API_KEY else "none - dev-mode instance only")
2. Locate the shipped open data¶
Petala ships a small open-data catalog under data/catalog/. This notebook uses
poi/id_yo_health.geojson — health facilities in Daerah Istimewa Yogyakarta,
extracted from OpenStreetMap (ODbL; see LICENSES-DATA.md). Nothing here is
customer data.
The default path assumes the notebook runs from docs/notebooks/. Set
PETALA_CATALOG if you moved it. The assert is deliberate: a missing catalog
should stop the notebook, not silently produce an empty map.
from pathlib import Path
CATALOG = Path(os.environ.get("PETALA_CATALOG", "../../data/catalog")).resolve()
assert CATALOG.is_dir(), f"catalog not found at {CATALOG} - set PETALA_CATALOG"
POI_FILE = CATALOG / "poi" / "id_yo_health.geojson"
# Tail of the path only: an absolute path is noise in a committed output.
print("catalog:", Path(*CATALOG.parts[-2:]), "|", POI_FILE.name,
f"| {POI_FILE.stat().st_size / 1024:.0f} KiB")
3. Get or create the project¶
p.project(name) is get-or-create by exact name: it lists projects, returns the
match, and only creates when there is none. Re-running this notebook therefore
reuses the same workspace instead of accumulating copies of it.
PROJECT = "Notebook demo"
proj = p.project(PROJECT)
print(f"project {proj['name']!r} -> {proj['id']}")
print("layers already here:", [layer["name"] for layer in p.layers(proj)])
4. Upload the layer — idempotently¶
external_ref is your own stable identifier for the layer, and mode="replace"
tells the server to swap the features of the layer carrying that ref rather than
create a second one. The layer keeps its id, name and style, so anything already
pointing at it — a map, a scoring factor, an embed — keeps working.
replaced in the response tells you which happened: False on the first run,
True on every run after. That is the whole idempotency story, and it is why
this notebook is safe to re-run.
poi = p.upload(
proj,
POI_FILE,
name="Health facilities (OSM, Yogyakarta)",
external_ref="notebooks/id_yo_health",
mode="replace",
)
print("layer ", poi["name"], "->", poi["id"])
print("geometry", poi["geom_type"], "| features:", poi["feature_count"])
print("replaced", poi["replaced"], "| skipped null geometry:", poi["skipped_null_geometry"])
5. A second layer: kabupaten population¶
The catalog's demography/idn_adm2_population.geojson is national — 522
kabupaten, ~5 MB. We only need the five in Daerah Istimewa Yogyakarta, so filter
locally and push the subset with import_geojson, the same external_ref +
replace contract the file upload uses.
Dropping to four columns is not only tidiness: the source carries date columns that are not JSON-serialisable once geopandas has parsed them.
import json
import geopandas as gpd
adm2 = gpd.read_file(CATALOG / "demography" / "idn_adm2_population.geojson")
diy = adm2.loc[
adm2["adm1_name"] == "Daerah Istimewa Yogyakarta",
["adm2_name", "population", "area_sqkm", "geometry"],
]
kab = p.import_geojson(
proj,
"Kabupaten population (DIY)",
json.loads(diy.to_json()),
external_ref="notebooks/idn_adm2_population_diy",
mode="replace",
)
print("layer", kab["name"], "->", kab["id"], "| features:", kab["feature_count"],
"| replaced:", kab["replaced"])
diy.drop(columns="geometry").sort_values("population", ascending=False)
6. Define the model¶
A factor is one reading at a point, and there are three kinds:
| kind | reads | needs |
|---|---|---|
density |
how many features of a layer lie within radius_m |
radius_m |
proximity |
metres to the nearest feature of a layer | max_distance_m (optional clamp) |
attribute |
a numeric property of the feature containing the point | property |
direction says which way is good — lower_better flips the normalized value,
which is how "closer is better" is expressed for a distance. Weights are plain
non-negative numbers and do not have to sum to 1: the score is the weighted
sum divided by the total weight, so it lands in 0–1 regardless.
There is no update call on the client yet, so this is get-or-create by name. The
layer ids stay valid across re-runs precisely because the uploads above used
mode="replace".
MODEL = "Clinic siting demo"
config = {
"normalization": "minmax",
"factors": [
{
"name": "health_density_2km",
"kind": "density",
"layer_id": poi["id"],
"params": {"radius_m": 2000},
"direction": "higher_better",
"weight": 2.0,
},
{
"name": "nearest_facility_m",
"kind": "proximity",
"layer_id": poi["id"],
"params": {"max_distance_m": 5000},
"direction": "lower_better",
"weight": 1.0,
},
{
"name": "kabupaten_population",
"kind": "attribute",
"layer_id": kab["id"],
"params": {"property": "population"},
"direction": "higher_better",
"weight": 1.0,
},
],
}
existing = [m for m in p.scoring_models(proj) if m["name"] == MODEL]
model = existing[0] if existing else p.create_scoring_model(proj, MODEL, config)
print(f"model {model['name']!r} -> {model['id']}")
for factor in model["config"]["factors"]:
print(f" {factor['name']:<22} {factor['kind']:<10} w={factor['weight']:<4} {factor['direction']}")
7. Score two points¶
Malioboro is the centre of Yogyakarta city. The second point sits on the Gunungkidul coast, about 40 km south-east — same province, very different surroundings.
refresh_domains=True on the first call recomputes each factor's min/max from
the layers as they stand right now. Domains are cached against the layer as it
was, and we just replaced both layers, so the first read should not trust the
cache. Subsequent calls can.
import pandas as pd
PROBES = {
"Malioboro (city centre)": (110.3656, -7.7925),
"Gunungkidul coast (remote)": (110.6300, -8.0700),
}
results = {}
for i, (label, (lng, lat)) in enumerate(PROBES.items()):
results[label] = p.score(model, lng=lng, lat=lat, refresh_domains=(i == 0))
for label, r in results.items():
print(f"{label:<28} score = {r['score']:.3f} skipped: {r['skipped'] or 'none'}")
The breakdown is where a score becomes explainable. raw is the reading in its
own unit, normalized is that reading mapped into 0–1 against the factor's
domain (already direction-flipped), and weighted is what it contributed.
rows = []
for label, r in results.items():
for f in r["factors"]:
rows.append({
"location": label,
"factor": f["name"],
"raw": f["raw"],
"normalized": round(f["normalized"], 4) if f["normalized"] is not None else None,
"weighted": round(f["weighted"], 4) if f["weighted"] is not None else None,
})
pd.DataFrame(rows).pivot(index="factor", columns="location",
values=["raw", "normalized", "weighted"])
Read the middle factor first: the city point sits 18 m from a facility and has ~100 within 2 km, the coastal point has none within 2 km and nothing inside the 5 km clamp. That is the whole gap.
kabupaten_population pulls the other way — Kota Yogyakarta is the smallest
DIY kabupaten by headcount because it is tiny in area, so raw population
punishes the city point. It is an honest reminder that a raw catalog attribute
is not automatically the factor you meant: population density, or a
catchment sum, would say what was intended here.
8. Score a grid¶
score_grid runs the same model over a hexagon lattice and writes the result as
a new layer. Each cell carries every factor's normalized value plus score (the
weighted sum) and score_norm (that sum over the total weight).
One caveat, stated in the product itself: grid cells are ranked against each other; the pin is rated against its layers as a whole. The two scores are on different scales — do not compare them. A grid min-max-normalizes over the cells it just computed; a pin uses the model's cached domains.
The bbox is [west, south, east, north] and is chosen to contain both probe
points.
SCORE_GRID_LAYER = "Clinic siting demo (1 km score grid)"
BBOX = [110.30, -8.10, 110.66, -7.74]
def layer_named(project, name):
for layer in p.layers(project):
if layer["name"] == name:
return layer
return None
scored = layer_named(proj, SCORE_GRID_LAYER)
if scored is None:
scored = p.score_grid(model, cell_m=1000, bbox=BBOX, name=SCORE_GRID_LAYER)
print("cells:", scored["cell_count"], "| projected in", scored["crs"],
"| bbox from", scored["bbox_source"])
else:
print("reusing existing score grid", scored["id"])
grid_gdf = p.read(scored)
print(grid_gdf.shape, "|", ", ".join(c for c in grid_gdf.columns if c != "geometry"))
grid_gdf[["score", "score_norm"]].describe().T
9. Map the score¶
One hue, light to dark: score_norm is a magnitude, and the only question the
map has to answer is where it is high. Both probe points are marked so the pin
numbers above can be located on the surface — remembering that the cell colours
rank against each other, not against the pins.
%matplotlib inline
import matplotlib.pyplot as plt
# Recessive axes so the data carries the figure, not the chrome.
plt.rcParams.update({
"figure.dpi": 110,
"savefig.dpi": 110,
"font.size": 9,
"axes.edgecolor": "#c8ccd4",
"axes.labelcolor": "#3c4350",
"text.color": "#3c4350",
"xtick.color": "#6b7280",
"ytick.color": "#6b7280",
"axes.spines.top": False,
"axes.spines.right": False,
})
import matplotlib.patheffects as pe
# Which side each label sits on, so neither runs into the colour bar.
SIDES = {"Malioboro (city centre)": (11, "left"),
"Gunungkidul coast (remote)": (-11, "right")}
fig, ax = plt.subplots(figsize=(7.4, 6.2))
grid_gdf.plot(ax=ax, column="score_norm", cmap="Greens", vmin=0, vmax=1,
linewidth=0.1, edgecolor="#ffffff",
legend=True,
legend_kwds={"label": "score_norm (rank within this grid)", "shrink": 0.72})
for label, (lng, lat) in PROBES.items():
dx, ha = SIDES[label]
ax.plot(lng, lat, marker="o", markersize=9, linestyle="none",
markerfacecolor="#ffffff", markeredgecolor="#1f2933",
markeredgewidth=1.6, zorder=5)
ax.annotate(label.split(" (")[0], (lng, lat), textcoords="offset points",
xytext=(dx, 7), ha=ha, fontsize=8, color="#1f2933", zorder=6,
path_effects=[pe.withStroke(linewidth=2.5, foreground="#ffffff")])
ax.set_title(f"Clinic siting score, 1 km hexagons ({len(grid_gdf)} cells)")
ax.set_xlabel("longitude")
ax.set_ylabel("latitude")
fig.tight_layout()
Where this goes next¶
The weights here were typed by hand, which is exactly what every generic siting score does. The version worth having fits them to outcomes you already own — revenue per outlet, claims per clinic — so the model states what actually predicted success in your network rather than what seemed reasonable.
Everything above is the same HTTP surface a scheduled job would use: keep the
external_refs, keep the model name, and re-running is a refresh rather than a
pile of duplicates.