How to Place Vertiport Locations in Any City Using Geospatial Machine Learning

airfield built for aircraft that take off and land straight up and down, such as electric air taxis. If a city has population, road, and airspace data but no reliable information about where people travel, how can we identify the best locations for placing vertiports?

That is the situation in Lagos, Nigeria, and in many fast-growing cities. The first attempt is obvious: use the population map, find the busiest areas, and place proposed vertiports there. To test that idea, I used a population raster, a map divided into small squares, each holding an estimated population count. I turned each square into one point, giving more importance to squares with more people, and ran K-means, a method that groups nearby points and places one center in each group.

The result spread across the parts of the city where the most people live. At a glance, it looked like a reasonable answer.

Here is what that first attempt produced. The left panel shows the population-only result; the right panel shows the later version after adding the screening rules.

The model was not failing randomly. It did exactly what I asked it to do: optimize population concentration. The problem was that population density alone was an incomplete objective. It did not tell me whether a location was accessible, buildable, or safe for aircraft.

The first result is not a usable answer. When I checked the proposed locations against two simple rules, some fell inside a safety distance drawn around Murtala Muhammed International Airport, while others landed on water or wetland. The model found groups of people, not usable sites. That gap, between where people are and where a vertiport can actually go, is the central problem this article solves.

This article uses a repeatable method that starts with population as a stand-in for travel demand, since real trip data is not available. It then adds closeness to roads and transit hubs, removes locations that fail a set of screening rules, groups the remaining points using real-world distances instead of raw map degrees, and checks every final location again after spacing adjustments. The output is not a construction plan, but a ranked shortlist another reader can reproduce and investigate further.

Lagos is a useful test case because its data is incomplete. There is no clean public table showing where people travel between neighborhoods, so this article shows how to build a useful first-pass model without pretending that population counts are the same as travel demand. The airport and military circles used later in the article are illustrative screening buffers, not official airspace maps. A real deployment study would replace them with geometry supplied by the aviation authorities.

The short version

The approach follows four decisions:

1. Where are the people? Use a population map as a starting estimate of likely demand.
2. Which places are easier to reach? Give a small boost to locations near major roads and transit hubs. This is the accessibility score: a simple measure of how close a place is to useful transport connections.
3. Which places should be ruled out? Remove water, wetland, mangrove, and locations inside the configured airport or military buffers.
4. How do the remaining places get grouped and ranked? Use K-means, the grouping method defined above, then check and rank every resulting site.

The order matters because a clustering algorithm is good at finding concentration but not at understanding airports, land ownership, flood risk, or aviation law, so those rules have to be made explicit instead of left for the algorithm to infer.

One vocabulary note before going further: throughout this article, demand means population count, used as a stand-in for real travel numbers, not measured trips. An exclusion is a place the model removes before choosing sites. A proposed site is a location worth checking next, not a confirmed one. Weighted means some map squares count more than others because they contain more people or have better access.

For the first full comparison, I used 100 proposed sites so the maps and validation table would show a substantial network. The site count analysis later in the article tests whether 100 is justified.

Now the map has context. The red crosses are not final recommendations. They are the result of combining the four decisions above, and the rest of the article explains how each decision changes the result.

What you will build

By the end, you will have a working first-pass siting pipeline you can point at another city. You will be able to:

• Assemble a boundary, population surface, land cover layer, road network, and transit hubs from public sources.
• Build a population-weighted demand surface from a raster, with no synthetic jittering: every demand point is a real pixel from the source data.
• Add transport and road accessibility scoring, so sites near existing transit rank higher than sites that only have raw population behind them.
• Add hard exclusion layers for airports, military installations, and unbuildable land such as open water, wetland, and mangrove.
• Generate and rank proposed sites by how much weighted demand each one actually captures, not just by raw population.
• Swap in a different city’s boundary, population raster, roads, transit stops, and constraints, and rerun the identical script.

The practical question behind these four decisions is how to turn incomplete city data into a shortlist that is useful without overstating what the model knows. The full script is in the Complete pipeline section, and the companion repository has the settings, data files, and sample outputs used here.

Before running the model, let us look at the two inputs it will use first: Lagos’s study area and the population map beneath it. Nothing has been grouped or ranked yet.

First try: use population alone

The first version asks one simple question: where do the most people live? It places proposed locations near the busiest parts of the population map.

To do that, the method treats each populated square as a small area of demand. A square estimated to contain 800 people matters eight times more than a square estimated to contain 100 people. The method then looks for nearby areas with many people and places one proposed location near the center of each area. I told it to return 100 locations because I wanted to test what a 100-site network would look like. That number was my test setting, not a conclusion discovered from the data.

This first version knows only where people live. It does not know that a location may be near an airport, on a lagoon, inside a military area, or far from a usable road.

The result shows the problem clearly. Eleven of the 100 proposed locations fail the screening checks: some fall inside the airport buffer near Ikeja, while others land on water or wetland. Adding roads and transport hubs without adding the screening checks does not solve the problem. It changes which busy areas receive more attention, but it still does not know which locations must be ruled out.

The lesson is simple: a busy area is not automatically a suitable site. The rest of the article adds the missing information one decision at a time.

What information do we need?

To choose good vertiport locations, we need several kinds of information. Each one tells us something different about the city:

One important map is missing: flood risk. A PIAHS study on flood vulnerability in Lagos shows why flood exposure matters, but I did not find a downloadable flood map suitable for this run. Flood risk is therefore a limitation of this first pass and should be added before using the method for real planning.

Building the pipeline: eight steps, any city

The method behind the Lagos map above has eight steps. The city specific inputs can change, but the steps stay the same:

Every city-specific detail, including the boundary file, population map, roads, transport stops, screening rules, number of sites, and spacing, lives in a JSON (JavaScript Object Notation) configuration file. JSON is simply a plain text settings file. None of these city details is hidden inside the Python functions, which is what makes the Nairobi run possible without changing the main script.

Steps 1 and 2: map where people live

When detailed population data is unavailable, one possible shortcut is to place a point at the center of each neighborhood and randomly spread copies around it. That creates the appearance of precise locations, but the extra points are invented. We avoid that shortcut. Instead, every populated square in the population map becomes one point, and its importance matches the population estimated in that square:

def load_population_points(paths: Paths, config: dict, boundary: gpd.GeoDataFrame) -> pd.DataFrame:
    """Turn every valid raster pixel inside the city boundary into a
    (lon, lat, population) demand point. No jittering: every point is a real
    pixel from the source raster."""
    raster_path = paths.root / config["population_raster"]
    with rasterio.open(raster_path) as src:
        boundary_reproj = boundary.to_crs(src.crs)
        geoms = [g.__geo_interface__ for g in boundary_reproj.geometry]
        clipped, transform = mask(src, geoms, crop=True, nodata=src.nodata)
        band = clipped[0]
        rows, cols = np.where(band != src.nodata)
        pop = band[rows, cols]
        xs, ys = rasterio.transform.xy(transform, rows, cols)

    df = pd.DataFrame({"lon": xs, "lat": ys, "population": pop})
    return df[df["population"] > 0].reset_index(drop=True)

For Lagos, this produces 133,490 real demand points totaling roughly 9.75 million people, a plausible bottom-up estimate for Lagos State. For Nairobi, the same function against Kenya’s population raster produces 45,942 demand points totaling roughly 5.38 million people across Nairobi County’s 17 sub-counties.

Step 3: favor places people can reach

A site with 50,000 people around it but no nearby transit is a worse choice than a site with 40,000 people right next to a rail station, because a vertiport only helps people who can actually reach it. The pipeline scores every demand point’s distance to the nearest transit hub and the nearest major road, using real OpenStreetMap data: for Lagos, 50 rail and bus rapid transit stations plus 124 ferry terminals, and 1,295 motorway, trunk, and primary road segments. For Nairobi, 20 commuter rail stations and 1,017 major road segments.

def add_accessibility_features(df: pd.DataFrame, paths: Paths, config: dict) -> pd.DataFrame:
    """Add transport_hub_km and major_road_km columns: distance from each
    demand point to the nearest transport hub / major road, in kilometers.
    Skips a column entirely if the config does not supply that layer."""
    df = df.copy()
    if config.get("transport_hubs"):
        hubs = gpd.read_file(paths.root / config["transport_hubs"])
        hub_lons = hubs.geometry.x.to_numpy()
        hub_lats = hubs.geometry.y.to_numpy()
        pt_lons = df["lon"].to_numpy()
        pt_lats = df["lat"].to_numpy()
        chunk = 5000
        nearest = np.empty(len(df))
        for start in range(0, len(df), chunk):
            end = start + chunk
            d = haversine_km(pt_lons[start:end, None], pt_lats[start:end, None], hub_lons[None, :], hub_lats[None, :])
            nearest[start:end] = d.min(axis=1)
        df["transport_hub_km"] = nearest
    if config.get("road_network"):
        roads = gpd.read_file(paths.root / config["road_network"])
        metric_crs = roads.estimate_utm_crs()
        pts_metric = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["lon"], df["lat"]), crs="EPSG:4326").to_crs(metric_crs)
        roads_metric = roads.to_crs(metric_crs)
        nearest = gpd.sjoin_nearest(pts_metric, roads_metric[["geometry"]], distance_col="dist_m")
        nearest = nearest[~nearest.index.duplicated(keep="first")].sort_index()
        df["major_road_km"] = (nearest["dist_m"] / 1000.0).to_numpy()
    return df

This function does the actual measuring: for every one of the 133,490 demand points, it finds the straight line distance to the nearest transit hub and the nearest major road. The next function turns those two distances into a single weight:

def compute_weighted_demand(df: pd.DataFrame, weights: dict) -> pd.Series:
    """Combine population with accessibility bonuses into one weight used
    both to fit K means and to rank the resulting sites. Each accessibility
    signal contributes a bonus that decays smoothly with distance, so a
    point right next to a hub gets close to the full weight bump and a
    point far away gets close to none."""
    decay_km = weights.get("accessibility_decay_km", 3.0)
    multiplier = np.ones(len(df))
    if "transport_hub_km" in df.columns and weights.get("transport_access", 0):
        multiplier += weights["transport_access"] * np.exp(-df["transport_hub_km"].to_numpy() / decay_km)
    if "major_road_km" in df.columns and weights.get("road_access", 0):
        multiplier += weights["road_access"] * np.exp(-df["major_road_km"].to_numpy() / decay_km)
    return df["population"] * multiplier

The distance signal here is straight line distance to the nearest hub or road, not routed travel time along the actual street network, and the difference can be large: two points 500 meters apart in straight line distance could be a 15 minute walk apart if a canal or a highway sits between them with no crossing nearby. Building genuine routed travel time would mean running a routing engine such as OSRM or Valhalla over the full road network, which this version of the pipeline does not do. Straight line distance to transit and roads is a reasonable first approximation, and it is explicitly labeled as one, not as measured travel time.

Step 4: remove places that fail a basic screen

This is the step that turns those 11 unsafe sites into zero. Two kinds of exclusion run here. Point buffer exclusions remove any demand point within a set radius of a named hazard, such as an airport. Land cover exclusions remove any point whose ESA WorldCover class is water, wetland, or mangrove.

def apply_point_buffer_exclusion(df: pd.DataFrame, boundary: gpd.GeoDataFrame, layer: dict) -> tuple[pd.DataFrame, gpd.GeoDataFrame]:
    """Drop demand points within `buffer_km` of any point the layer lists
    (airports, military installations, or any other named hazard point)."""
    metric_crs = boundary.estimate_utm_crs()
    points_gdf = gpd.GeoDataFrame(
        layer["points"],
        geometry=[Point(p["lon"], p["lat"]) for p in layer["points"]],
        crs="EPSG:4326",
    )
    buffers_metric = points_gdf.to_crs(metric_crs)
    buffers_metric["geometry"] = buffers_metric.geometry.buffer(layer["buffer_km"] * 1000)
    buffers_wgs84 = buffers_metric.to_crs("EPSG:4326")

    demand_points = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["lon"], df["lat"]), crs="EPSG:4326")
    union_buffer = buffers_wgs84.geometry.union_all()
    inside = demand_points.geometry.within(union_buffer)
    return df.loc[~inside.values].reset_index(drop=True), buffers_wgs84

EPSG:4326 Above is the standard GPS format for latitude and longitude. A coordinate reference system is simply the rule used to turn a place on Earth into numbers. A distance such as five kilometers cannot be drawn reliably in GPS degrees, so the code temporarily converts the map to a local meter-based system before drawing each buffer, then converts back for display.

For Lagos, the point locations come from named airport and military features in the project data. The circles around them are configurable screening buffers, not official no-fly zones. Together with the ESA WorldCover water, wetland, and mangrove exclusion, these rules removed 13,610 of Lagos’s 133,490 demand points, about 10 percent, before clustering happened. This is deliberately cautious, because a first pass tool should flag a site for review rather than imply that a dense population pocket is automatically buildable.

Steps 5 and 6: combine the evidence and group places

Once accessibility scoring and exclusions are both in place, population and accessibility combine into a single number per point, weighted_demand. Population supplies the demand signal; proximity to transit and major roads supplies modest bonuses that decay with distance. These weights are assumptions, not learned truths, which is why the article later sweeps the settings instead of presenting one map as definitive.

There is one detail that changes the results here: K means measures straight line distance, but latitude and longitude are angular coordinates, not flat map distances, so fitting on raw degrees can distort cluster shapes. The fix converts the points into a local meter based map before fitting, then converts the resulting centers back to latitude and longitude for reporting:

def to_metric_xy(df: pd.DataFrame, metric_crs) -> np.ndarray:
    points = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["lon"], df["lat"]), crs="EPSG:4326").to_crs(metric_crs)
    return np.column_stack([points.geometry.x.to_numpy(), points.geometry.y.to_numpy()])


def site_vertiports(df: pd.DataFrame, k: int, seed: int, metric_crs) -> tuple[np.ndarray, np.ndarray]:
    X = to_metric_xy(df, metric_crs)
    model = KMeans(n_clusters=k, random_state=seed, n_init=10)
    labels = model.fit_predict(X, sample_weight=df["weighted_demand"].to_numpy())
    centers_metric = gpd.GeoDataFrame(
        geometry=gpd.points_from_xy(model.cluster_centers_[:, 0], model.cluster_centers_[:, 1]), crs=metric_crs
    ).to_crs("EPSG:4326")
    centers = np.column_stack([centers_metric.geometry.x.to_numpy(), centers_metric.geometry.y.to_numpy()])
    return labels, centers

K means here is doing one specific job: finding where the remaining, already filtered demand is spatially concentrated. It is not deciding whether a location is buildable, connected, or legal. Those decisions happen in the exclusion stage and in the final validation stage.

A K means problem worth checking

Even after every demand point going into K means is screened, a cluster’s centroid can still be unsafe. A centroid is a weighted average position, not one of the members itself. If safe points sit on opposite sides of a narrow excluded strip, their average can land inside it.

The fix is a verification step that checks every centroid against the city boundary and the same exclusion layers one more time. Both this check and the one after merging in the next step share the same underlying test, a function that checks a batch of sites and returns which ones are actually safe:

def screened_site_mask(sites: pd.DataFrame, paths: Paths, boundary: gpd.GeoDataFrame,
                       exclusion_layers: list[dict]) -> np.ndarray:
    """Return one True/False value per site, preserving the original IDs."""
    probe = sites[["lon", "lat"]].copy().reset_index(drop=True)
    probe["population"] = 1.0
    probe["__site_id"] = probe.index
    safe_probe, _, _ = apply_exclusion_layers(probe, paths, boundary, exclusion_layers)
    safe_ids = set(safe_probe["__site_id"])
    study_area = boundary.geometry.union_all()
    return np.array([
        site_id in safe_ids and study_area.covers(Point(row["lon"], row["lat"]))
        for site_id, row in probe.iterrows()
    ])

The study_area.covers(...) check is doing real work here: it catches the case where a centroid drifts entirely outside the city boundary, something the exclusion layers alone would never flag, since they only know about named hazards, not the edge of the study area itself. Any site that fails this check gets replaced with the nearest real, already screened demand point from its own cluster:

def snap_unsafe_sites(sites: pd.DataFrame, demand: pd.DataFrame, labels: np.ndarray, paths: Paths,
                       boundary: gpd.GeoDataFrame, exclusion_layers: list[dict]) -> pd.DataFrame:
    """Move any site whose own coordinate fails the exclusion check to the
    nearest real, already-verified-safe member point in its cluster."""
    sites = sites.reset_index(drop=True)
    safe_mask = screened_site_mask(sites, paths, boundary, exclusion_layers)
    unsafe_ids = set(np.where(~safe_mask)[0])

    sites = sites.copy()
    for site_id in unsafe_ids:
        members = demand.loc[labels == site_id]
        if members.empty:
            continue
        site_lon, site_lat = sites.loc[site_id, "lon"], sites.loc[site_id, "lat"]
        dist = haversine_km(site_lon, site_lat, members["lon"].to_numpy(), members["lat"].to_numpy())
        nearest = members.iloc[int(np.argmin(dist))]
        sites.loc[site_id, "lon"] = nearest["lon"]
        sites.loc[site_id, "lat"] = nearest["lat"]
    return sites

This step is not the whole safety story, since the minimum spacing step can create the same problem again.

Step 7: stop proposed sites from sitting on top of one another

Two vertiports 200 meters apart is not two proposed sites, it is one site reported twice. After clustering, any pair of sites closer than a configured minimum distance gets merged into a single demand weighted centroid, and the process repeats until every remaining pair clears the minimum. At Lagos’s 1.5 kilometer minimum spacing, all 100 requested sites survived without needing a merge. A separate test run with 300 requested sites at a 2 kilometer minimum merged down to 201, which is the expected direction: tighter spacing at higher density triggers more merges. The full merge function is in the complete script below.

Merging two sites has the same blind spot, because the merged position is another weighted average, so it can land in an excluded zone even when both parent sites were safe. The pipeline therefore checks every site again after merging and moves any failing site to the nearest screened demand point in the filtered surface.

def snap_to_nearest_safe_point(sites: pd.DataFrame, demand: pd.DataFrame, paths: Paths,
                                boundary: gpd.GeoDataFrame, exclusion_layers: list[dict]) -> pd.DataFrame:
    sites = sites.reset_index(drop=True)
    safe_mask = screened_site_mask(sites, paths, boundary, exclusion_layers)
    unsafe_ids = set(np.where(~safe_mask)[0])
    study_area = boundary.geometry.union_all()

    sites = sites.copy()
    for site_id in unsafe_ids:
        site_lon, site_lat = sites.loc[site_id, "lon"], sites.loc[site_id, "lat"]
        inside_boundary_mask = [study_area.covers(Point(lon, lat)) for lon, lat in zip(demand.lon, demand.lat)]
        safe_demand = demand.loc[inside_boundary_mask].reset_index(drop=True)
        if safe_demand.empty:
            continue
        dist = haversine_km(site_lon, site_lat, safe_demand["lon"].to_numpy(), safe_demand["lat"].to_numpy())
        nearest = safe_demand.iloc[int(np.argmin(dist))]
        sites.loc[site_id, "lon"] = nearest["lon"]
        sites.loc[site_id, "lat"] = nearest["lat"]
    return sites

Running this check against both this article’s final results confirms zero sites fail the exclusion test after both safety nets, the centroid check right after clustering and this one right after merging, which is what the metrics files in the companion repository actually record, not an assumption about how the method should behave.

Step 8: rank the shortlist

The last step sorts sites by how much weighted demand each one actually captures, so the output is a priority order, not an unordered dump of coordinates. The table below shows the top sites from the 100 site comparison scenario used in this section:

Rank 6 shows why ranking by weighted demand, not raw population, matters: it captures more raw population than rank 5, 166,958 versus 156,208, but sits 4.10 kilometers from the nearest transit hub instead of 1.24, so its accessibility adjusted score comes out lower and it ranks behind a smaller crowd that is easier to actually reach. The ranking is doing exactly what it is supposed to do, trading some raw population coverage for meaningfully better transit access.

Here is the 100 site scenario in full, with the top 6 from the table above marked directly on the map. The next section tests whether Lagos needs this many sites.

Validating the method: three versions compared

Showing one failed population only attempt and one improved map is useful, but it is not proof that the improvement works across the full set of proposed sites. To check that, I run three complete versions of the method on the same city and compare every site, not just a few.

Read each column in the direction that matters. “Sites checked unsafe” counts proposed sites outside the study boundary or inside one of this pipeline’s screening rules; lower is better, and variant C’s zero is the point of adding those rules. “Population captured” is the raw number of people living near the top 20 proposed sites; higher is generally useful, but it is not automatically better if those people are near a place the screening rules remove. “Mean distance to transit” and “mean distance to a road” are straight line distances in kilometers from the population points in the top 20 groups to the nearest hub or road; lower means better connected.

The accessibility column does not move in a single clean direction, and that is worth reporting honestly rather than smoothing over. Variant B improves mean transit distance over variant A, 3.17 kilometers versus 3.48, but its mean road distance actually gets slightly worse, 1.31 kilometers versus 1.21, because the areas with the strongest transit access are not always the areas closest to a major road. Variant C is the only version that improves both at once, 2.88 kilometers to transit and 1.16 to a road, which happens partly because removing unsafe, densely populated areas near the airport redistributes the top ranked sites toward other well connected parts of the city.

The population tradeoff turned out smaller than a first look at that comparison would suggest. Variant C captures about 1.7 percent less population in its top 20 sites than variant A does, 2.92 million versus 2.97 million. The real cost shows up against variant B instead: accessibility weighting without exclusions captures about 9.5 percent more population than the fully screened version, 3.23 million versus 2.92 million. That is the cost of requiring proposed sites to pass the screening rules, and it is exactly the tradeoff a planner needs to see.

How sensitive is this to the settings you pick

The pipeline was rerun at several values for two settings, to check how much each one moves the final answer: the number of requested sites, and the airport exclusion buffer radius.

Neither result is surprising once you see it, and that is exactly why it is worth showing rather than assuming: a planner picking K or a buffer radius is not making a free choice. Every unit of extra safety margin or extra site count has a measurable population cost, and this pipeline reports that cost instead of hiding it behind a single confident looking map.

How many proposed sites does Lagos need?

The number 100 was useful for comparing maps, but it was not an answer. K means requires a number of groups before it starts, so choosing 100 at the beginning would simply make the method repeat our assumption. To choose the number more carefully, I ran the constrained method with 25, 50, 75, 100, 125, 150, and 200 proposed sites.

I used two checks. The first asks what happens when we add more sites. For each possible site count, I measured how far each populated map square was from the proposed site serving it. Adding more sites should bring people closer to a site, but each extra site should eventually produce a smaller improvement. The point where the gains begin to level off is often called the elbow. It helps show when adding more sites may no longer be worth the extra cost, but it is not precise enough to make the decision by itself.

The second check asks a question that is easier to explain: how close is the nearest proposed site to the people it is meant to serve? I measured the average distance and the share of the screened population within 5 kilometers. The 5 kilometer threshold is an assumption for this demonstration, not a universal planning standard. I selected the smallest tested network that covered at least 90 percent of the screened population within 5 kilometers while keeping the average nearest site distance below 4 kilometers.

The result is 25 proposed sites under these assumptions. At 25 sites, the average nearest site distance is 3.37 kilometers and 92.2 percent of the screened population is within 5 kilometers. More sites improve the numbers, but they are not needed to meet the target: 50 sites raises 5 kilometer coverage to 98.3 percent, while the average distance falls to 2.32 kilometers. A city with a smaller service radius, a different budget, or a requirement for better access within 3 kilometers would reasonably choose a larger network.

This is why the article calls 25 the recommended number for this scenario, not the universally optimal number for Lagos. The data provides evidence about coverage and spacing; planners still have to provide the service standard, available land, construction budget, and aviation requirements. The method makes those assumptions visible and lets another city replace them in the configuration file.

The recommended 25 site result is shown below. Unlike the earlier 100 site map, this is the network produced after applying the coverage rule.

The top 6 sites in the recommended network are listed below, ranked the same way as the earlier 100 site example: by weighted demand, not raw population. The area name is the Local Government Area containing the model’s coordinate. It is a useful map reference, not a confirmed address or a promise that land is available there.

The rank labels on the map correspond to this table. A reader can reproduce each point by copying its latitude and longitude into a map, then checking the surrounding roads, land, flood exposure, ownership, and aviation restrictions before treating it as a real planning option.

The same pattern from Step 8 shows up here too. Rank 1 captures less raw population than rank 2, 622,095 versus 716,908, but it sits far closer to transit and a major road, 1.35 kilometers versus 4.66, so its accessibility adjusted score puts it first anyway. With only 25 sites instead of 100, each one now covers a much larger, more populated area, which is why these captured population numbers run several times higher than the 100 site table above.

Changing the city and rerunning, for real

The strongest test of a “reusable” pipeline claim is pointing the same, unmodified script at a second city’s data and reporting what happens, including the parts that do not go as expected, rather than a paragraph that only promises it works elsewhere.

Nairobi, Kenya was that second run. Its config file swaps in Kenyan boundary and population data, a completely different set of verified airports and a military airbase, and real Nairobi commuter rail stations and major roads, none of which required touching the pipeline’s Python code:

{
  "city_name": "Nairobi, Kenya",
  "city_boundary": "data/nairobi_subcounties.geojson",
  "population_raster": "data/nairobi_worldpop_gridded_v2_0.tif",
  "road_network": "data/nairobi_major_roads.geojson",
  "transport_hubs": "data/nairobi_transport_hubs.geojson",
  "exclusion_layers": [
    {
      "type": "point_buffer",
      "name": "airports_and_airbase",
      "buffer_km": 5.0,
      "points": [
        {"name": "Jomo Kenyatta International Airport (HKJK/NBO)", "lat": -1.3169486, "lon": 36.9288569},
        {"name": "Wilson Airport (HKNW/WIL)", "lat": -1.3241997, "lon": 36.8134204},
        {"name": "Moi Air Base (HKRE)", "lat": -1.2731598, "lon": 36.8652199}
      ]
    },
    {
      "type": "landcover",
      "name": "water_wetland_mangrove",
      "source": "esa-worldcover-2021",
      "cache_path": "data/nairobi_worldcover_2021.tif",
      "excluded_classes": [80, 90, 95]
    }
  ],
  "number_of_sites": 60,
  "minimum_site_distance_km": 1.5,
  "objective_weights": {
    "transport_access": 0.4,
    "road_access": 0.2,
    "accessibility_decay_km": 3.0
  }
}

Running python code/vertiport_siting.py --config code/configs/nairobi.json against this file produced 57 final sites from 60 requested, with a silhouette score of 0.33. A silhouette score is a rough measure of whether the groups are reasonably separate; higher is better, but it does not prove that the locations are useful. Zero Nairobi sites failed the configured screening check.

The buffer radius that was fine for Lagos was not fine for Nairobi

The Nairobi config reused Lagos’s 5 kilometer airport buffer as a starting point. That single unchanged number removed 14,195 of Nairobi’s 45,942 demand points, close to 31 percent. The same buffer radius in Lagos removed about 10 percent.

This happens because of geography, not because anything in the pipeline is broken. Jomo Kenyatta International Airport and Wilson Airport sit roughly 13 kilometers apart, both close to the center of a compact 17 sub county area. Their two 5 kilometer buffers overlap and together cover a large share of Nairobi’s total land area. Lagos’s airports are spread across a much larger, more elongated metropolitan footprint, so the same radius removes a smaller fraction of the whole.

That is a real, unresolved tension between two reasonable positions: keeping one fixed safety distance that travels with the pipeline versus recalibrating the radius to each city’s specific airspace geometry before the model ever sees the data. The honest answer probably depends on whether you are running a first pass planning exercise or preparing something an aviation regulator will actually review, and it is worth returning to once the rest of this pipeline’s limitations are on the table. The limitations section below has more on where that recalibration would need to happen, and the conclusion closes on this exact tradeoff.

What this does not model yet

Three gaps are worth stating plainly, because a tutorial that hides its own limitations is more dangerous than one that lists them.

The exclusion buffers are illustrative, not verified airspace. A 5 kilometer circle around an airport is a reasonable stand in for “somewhere near an active runway is unsafe,” but it is not the same as Nigeria’s or Kenya’s actual controlled airspace, which has an irregular shape tied to real flight paths, not a simple circle. A deployment ready version of this pipeline would need to source that shape data directly from each country’s civil aviation authority.

Distance to transit and roads is straight line, not routed travel time. As flagged earlier, a rider’s actual travel time depends on the street network, traffic, and whether a direct path even exists, none of which straight line distance captures.

No flood risk layer exists for Lagos in this pipeline, for the simple reason that no open, downloadable flood hazard dataset could be found during this project, only academic figures. Given that roughly a third of the area studied in the Lagos flood vulnerability paper cited earlier falls into a high risk category, this is arguably a larger gap than the airspace buffer’s imprecision, and it should be the first addition anyone extending this work for real deployment planning makes.

The complete pipeline: one script, any city

This is the full script behind every result in this article. It takes a config file as its only required argument and contains no city specific coordinates or filenames.

#!/usr/bin/env python3
"""Find proposed vertiport locations in any city using geospatial machine learning.

This script does not hard code Lagos. Every city specific input comes from a
JSON config file (see configs/lagos.json for the worked example):

  city_boundary        path to a polygon file (GeoJSON/shapefile) for the city
  population_raster    path to a gridded population raster clipped to the city
  road_network         optional path to a line layer of major roads
  transport_hubs       optional path to a point layer of rail/BRT/ferry stops
  exclusion_layers      list of hard exclusion rules, each either:
                          {"type": "point_buffer", "points": [...], "buffer_km": n}
                          {"type": "landcover", "source": "esa-worldcover-2021",
                           "excluded_classes": [...], "cache_path": "..."}
  number_of_sites       how many vertiports to place (K for K-Means)
  site_count_values     site counts to compare before choosing a network size
  site_count_policy     coverage and distance targets for the recommendation
  minimum_site_distance_km  minimum spacing enforced between final sites
  objective_weights     {"transport_access": w, "road_access": w,
                          "accessibility_decay_km": d}

The method, independent of city:

  1. Load the city boundary.
  2. Turn a population raster into weighted demand points (one point per
     pixel with population greater than zero -- no synthetic jittering).
  3. Score each demand point's transport and road accessibility, if those
     layers are supplied.
  4. Remove demand points inside any hard exclusion layer (airports,
     military installations, water, wetland, mangrove, or anything else the
     config lists).
  5. Combine population and accessibility into one weighted demand score per
     point, controlled by objective_weights.
  6. Run weighted K-Means to find `number_of_sites` demand clusters.
  7. Merge clusters closer than `minimum_site_distance_km` so no two final
     sites sit unrealistically close together.
  8. Rank the surviving sites by how much weighted demand each one captures.

To choose `number_of_sites` instead of assuming it, run
`python code/compare_variants.py --config code/configs/lagos.json`. That command
compares the values in `site_count_values`, marks the smallest one that meets
`site_count_policy`, and saves the supporting chart and CSV.

The central lesson this pipeline is built to demonstrate: it does not produce
one universally optimal set of vertiport locations. It produces the best
proposed locations under a specific objective, dataset, and set of
constraints -- change the config and the ranked list changes with it.

The worked example is Lagos, Nigeria (see configs/lagos.json and the article
for the caveats specific to that run: illustrative aviation buffers rather
than verified NCAA controlled airspace, and no open flood-risk dataset was
found for Lagos, so flood exposure is discussed qualitatively, not modeled).

Run from the article folder:

    python code/vertiport_siting.py --config code/configs/lagos.json
"""

from __future__ import annotations

import argparse
import json
from dataclasses import dataclass
from pathlib import Path
from typing import Any, Optional

import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio
from rasterio.mask import mask
from rasterio.merge import merge as rio_merge
from shapely.geometry import Point
from sklearn.cluster import KMeans
from sklearn.metrics import davies_bouldin_score, silhouette_score

SEED = 42
EARTH_RADIUS_KM = 6371.0088


@dataclass(frozen=True)
class Paths:
    root: Path
    data: Path
    outputs: Path
    media: Path

    @classmethod
    def from_root(cls, root: Path, city_slug: str) -> "Paths":
        outputs = root / "outputs" / city_slug
        media = root / "media" / city_slug
        outputs.mkdir(parents=True, exist_ok=True)
        media.mkdir(parents=True, exist_ok=True)
        return cls(root=root, data=root / "data", outputs=outputs, media=media)


def slugify_city(city_name: str) -> str:
    first_word = city_name.split(",")[0].strip().lower()
    return "".join(c if c.isalnum() else "_" for c in first_word).strip("_")


def haversine_km(lon1: np.ndarray, lat1: np.ndarray, lon2: np.ndarray, lat2: np.ndarray) -> np.ndarray:
    p1, p2 = np.radians(lat1), np.radians(lat2)
    dphi = np.radians(lat2 - lat1)
    dlambda = np.radians(lon2 - lon1)
    a = np.sin(dphi / 2) ** 2 + np.cos(p1) * np.cos(p2) * np.sin(dlambda / 2) ** 2
    return 2 * EARTH_RADIUS_KM * np.arcsin(np.sqrt(a))


def load_config(config_path: Path) -> dict:
    return json.loads(config_path.read_text())


def load_city_boundary(paths: Paths, config: dict) -> gpd.GeoDataFrame:
    gdf = gpd.read_file(paths.root / config["city_boundary"])
    if gdf.crs is None:
        gdf = gdf.set_crs(epsg=4326)
    return gdf.to_crs(epsg=4326)


def load_population_points(paths: Paths, config: dict, boundary: gpd.GeoDataFrame) -> pd.DataFrame:
    """Turn every valid raster pixel inside the city boundary into a
    (lon, lat, population) demand point. No jittering: every point is a real
    pixel from the source raster."""
    raster_path = paths.root / config["population_raster"]
    with rasterio.open(raster_path) as src:
        boundary_reproj = boundary.to_crs(src.crs)
        geoms = [g.__geo_interface__ for g in boundary_reproj.geometry]
        clipped, transform = mask(src, geoms, crop=True, nodata=src.nodata)
        band = clipped[0]
        rows, cols = np.where(band != src.nodata)
        pop = band[rows, cols]
        xs, ys = rasterio.transform.xy(transform, rows, cols)

    df = pd.DataFrame({"lon": xs, "lat": ys, "population": pop})
    return df[df["population"] > 0].reset_index(drop=True)


def ensure_worldcover_raster(paths: Paths, boundary: gpd.GeoDataFrame, layer: dict) -> Path:
    """Return a path to a city-clipped ESA WorldCover 2021 raster, fetching
    and clipping it from Microsoft Planetary Computer on first run if it is
    not already on disk. ESA WorldCover is public with no restriction of
    use; Planetary Computer just gates the storage account behind a short
    lived, free signed URL, so this call needs no API key."""
    out_path = paths.root / layer["cache_path"]
    if out_path.exists():
        return out_path

    import requests

    minx, miny, maxx, maxy = boundary.total_bounds
    search = requests.post(
        "https://planetarycomputer.microsoft.com/api/stac/v1/search",
        json={
            "collections": ["esa-worldcover"],
            "bbox": [minx, miny, maxx, maxy],
            "query": {"esa_worldcover:product_version": {"eq": "2.0.0"}},
        },
        timeout=30,
    ).json()

    srcs = []
    for feature in search["features"]:
        href = feature["assets"]["map"]["href"]
        signed = requests.get(
            "https://planetarycomputer.microsoft.com/api/sas/v1/sign",
            params={"href": href}, timeout=30,
        ).json()["href"]
        srcs.append(rasterio.open("/vsicurl/" + signed))

    mosaic, out_transform = rio_merge(srcs, bounds=(minx, miny, maxx, maxy))
    out_meta = srcs[0].meta.copy()
    out_meta.update({
        "height": mosaic.shape[1], "width": mosaic.shape[2],
        "transform": out_transform, "compress": "deflate", "predictor": 2,
    })
    out_path.parent.mkdir(parents=True, exist_ok=True)
    with rasterio.open(out_path, "w", **out_meta) as dst:
        dst.write(mosaic)
    for s in srcs:
        s.close()
    return out_path


def apply_point_buffer_exclusion(df: pd.DataFrame, boundary: gpd.GeoDataFrame, layer: dict) -> tuple[pd.DataFrame, gpd.GeoDataFrame]:
    """Drop demand points within `buffer_km` of any point the layer lists
    (airports, military installations, or any other named hazard point)."""
    metric_crs = boundary.estimate_utm_crs()
    points_gdf = gpd.GeoDataFrame(
        layer["points"],
        geometry=[Point(p["lon"], p["lat"]) for p in layer["points"]],
        crs="EPSG:4326",
    )
    buffers_metric = points_gdf.to_crs(metric_crs)
    buffers_metric["geometry"] = buffers_metric.geometry.buffer(layer["buffer_km"] * 1000)
    buffers_wgs84 = buffers_metric.to_crs("EPSG:4326")

    demand_points = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["lon"], df["lat"]), crs="EPSG:4326")
    union_buffer = buffers_wgs84.geometry.union_all()
    inside = demand_points.geometry.within(union_buffer)
    return df.loc[~inside.values].reset_index(drop=True), buffers_wgs84


def apply_landcover_exclusion(df: pd.DataFrame, raster_path: Path, excluded_classes: list[int]) -> tuple[pd.DataFrame, dict]:
    """Drop demand points that sample to an excluded land cover class, such
    as open water, wetland, or mangrove."""
    with rasterio.open(raster_path) as src:
        coords = list(zip(df["lon"], df["lat"]))
        classes = np.array([v[0] for v in src.sample(coords)])
    excluded_mask = np.isin(classes, excluded_classes)
    dropped_counts = {int(c): int((classes == c).sum()) for c in sorted(set(classes[excluded_mask].tolist()))}
    return df.loc[~excluded_mask].reset_index(drop=True), dropped_counts


def apply_exclusion_layers(df: pd.DataFrame, paths: Paths, boundary: gpd.GeoDataFrame,
                            exclusion_layers: list[dict]) -> tuple[pd.DataFrame, dict, dict]:
    """Apply every exclusion layer in the config in order. Returns the
    filtered demand points, the point-buffer geometries (for plotting), and
    a report of how many points each layer removed."""
    buffers_by_layer: dict[str, gpd.GeoDataFrame] = {}
    report: dict[str, Any] = {}
    for layer in exclusion_layers:
        n_before = len(df)
        if layer["type"] == "point_buffer":
            df, buffers = apply_point_buffer_exclusion(df, boundary, layer)
            buffers_by_layer[layer["name"]] = buffers
            report[layer["name"]] = n_before - len(df)
        elif layer["type"] == "landcover":
            raster_path = ensure_worldcover_raster(paths, boundary, layer)
            df, dropped_counts = apply_landcover_exclusion(df, raster_path, layer["excluded_classes"])
            report[layer["name"]] = dropped_counts
        else:
            raise ValueError(f"Unknown exclusion layer type: {layer['type']}")
    return df, buffers_by_layer, report


def screened_site_mask(sites: pd.DataFrame, paths: Paths, boundary: gpd.GeoDataFrame,
                       exclusion_layers: list[dict]) -> np.ndarray:
    """Return one True/False value per site, preserving the original IDs."""
    probe = sites[["lon", "lat"]].copy().reset_index(drop=True)
    probe["population"] = 1.0
    probe["__site_id"] = probe.index
    safe_probe, _, _ = apply_exclusion_layers(probe, paths, boundary, exclusion_layers)
    safe_ids = set(safe_probe["__site_id"])
    study_area = boundary.geometry.union_all()
    return np.array([
        site_id in safe_ids and study_area.covers(Point(row["lon"], row["lat"]))
        for site_id, row in probe.iterrows()
    ])


def add_accessibility_features(df: pd.DataFrame, paths: Paths, config: dict) -> pd.DataFrame:
    """Add transport_hub_km and major_road_km columns: distance from each
    demand point to the nearest transport hub / major road, in kilometers.
    Skips a column entirely if the config does not supply that layer."""
    df = df.copy()
    if config.get("transport_hubs"):
        hubs = gpd.read_file(paths.root / config["transport_hubs"])
        hub_lons = hubs.geometry.x.to_numpy()
        hub_lats = hubs.geometry.y.to_numpy()
        pt_lons = df["lon"].to_numpy()
        pt_lats = df["lat"].to_numpy()
        # Chunk the demand points so the (points x hubs) distance matrix
        # stays a manageable size in memory for large demand surfaces.
        chunk = 5000
        nearest = np.empty(len(df))
        for start in range(0, len(df), chunk):
            end = start + chunk
            d = haversine_km(pt_lons[start:end, None], pt_lats[start:end, None], hub_lons[None, :], hub_lats[None, :])
            nearest[start:end] = d.min(axis=1)
        df["transport_hub_km"] = nearest
    if config.get("road_network"):
        roads = gpd.read_file(paths.root / config["road_network"])
        metric_crs = roads.estimate_utm_crs()
        pts_metric = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["lon"], df["lat"]), crs="EPSG:4326").to_crs(metric_crs)
        roads_metric = roads.to_crs(metric_crs)
        nearest = gpd.sjoin_nearest(pts_metric, roads_metric[["geometry"]], distance_col="dist_m")
        nearest = nearest[~nearest.index.duplicated(keep="first")].sort_index()
        df["major_road_km"] = (nearest["dist_m"] / 1000.0).to_numpy()
    return df


def compute_weighted_demand(df: pd.DataFrame, weights: dict) -> pd.Series:
    """Combine population with accessibility bonuses into one weight used
    both to fit K-Means and to rank the resulting sites. Each accessibility
    signal contributes a bonus that decays smoothly with distance, so a
    point right next to a hub gets close to the full weight bump and a
    point far away gets close to none."""
    decay_km = weights.get("accessibility_decay_km", 3.0)
    multiplier = np.ones(len(df))
    if "transport_hub_km" in df.columns and weights.get("transport_access", 0):
        multiplier += weights["transport_access"] * np.exp(-df["transport_hub_km"].to_numpy() / decay_km)
    if "major_road_km" in df.columns and weights.get("road_access", 0):
        multiplier += weights["road_access"] * np.exp(-df["major_road_km"].to_numpy() / decay_km)
    return df["population"] * multiplier


def to_metric_xy(df: pd.DataFrame, metric_crs) -> np.ndarray:
    """Project (lon, lat) degrees into a local metric CRS. K-Means and the
    cluster quality metrics both use plain Euclidean distance, and a degree
    of longitude is not the same real-world distance as a degree of
    latitude except at the equator, so fitting on raw degrees silently
    distorts cluster shapes more and more as a city's latitude increases."""
    points = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["lon"], df["lat"]), crs="EPSG:4326").to_crs(metric_crs)
    return np.column_stack([points.geometry.x.to_numpy(), points.geometry.y.to_numpy()])


def site_vertiports(df: pd.DataFrame, k: int, seed: int, metric_crs) -> tuple[np.ndarray, np.ndarray]:
    X = to_metric_xy(df, metric_crs)
    model = KMeans(n_clusters=k, random_state=seed, n_init=10)
    labels = model.fit_predict(X, sample_weight=df["weighted_demand"].to_numpy())
    centers_metric = gpd.GeoDataFrame(
        geometry=gpd.points_from_xy(model.cluster_centers_[:, 0], model.cluster_centers_[:, 1]), crs=metric_crs
    ).to_crs("EPSG:4326")
    centers = np.column_stack([centers_metric.geometry.x.to_numpy(), centers_metric.geometry.y.to_numpy()])
    return labels, centers


def summarize_clusters(df: pd.DataFrame, labels: np.ndarray, centers: np.ndarray) -> pd.DataFrame:
    """Build one row per K-Means cluster: its centroid and how much
    population / weighted demand / accessibility it actually captures."""
    rows = []
    for cluster_id, (lon, lat) in enumerate(centers):
        members = df.loc[labels == cluster_id]
        row = {
            "lon": lon, "lat": lat,
            "captured_population": float(members["population"].sum()),
            "captured_weighted_demand": float(members["weighted_demand"].sum()),
        }
        if "transport_hub_km" in df.columns:
            row["mean_transport_hub_km"] = float(np.average(members["transport_hub_km"], weights=members["population"]))
        if "major_road_km" in df.columns:
            row["mean_major_road_km"] = float(np.average(members["major_road_km"], weights=members["population"]))
        rows.append(row)
    return pd.DataFrame(rows)


def snap_unsafe_sites(sites: pd.DataFrame, demand: pd.DataFrame, labels: np.ndarray, paths: Paths,
                       boundary: gpd.GeoDataFrame, exclusion_layers: list[dict]) -> pd.DataFrame:
    """A K-Means centroid is the weighted mean position of a cluster's
    members, so it can still land inside an excluded zone even when every
    member point feeding it is already outside every exclusion layer, for
    example when a cluster's safe points sit on both sides of a creek. Any
    site whose own coordinate fails the same exclusion check gets moved to
    the nearest real, already-verified-safe member point in its cluster
    instead of being reported at a location that was never actually
    checked."""
    sites = sites.reset_index(drop=True)
    safe_mask = screened_site_mask(sites, paths, boundary, exclusion_layers)
    unsafe_ids = set(np.where(~safe_mask)[0])

    sites = sites.copy()
    for site_id in unsafe_ids:
        members = demand.loc[labels == site_id]
        if members.empty:
            continue
        site_lon, site_lat = sites.loc[site_id, "lon"], sites.loc[site_id, "lat"]
        dist = haversine_km(site_lon, site_lat, members["lon"].to_numpy(), members["lat"].to_numpy())
        nearest = members.iloc[int(np.argmin(dist))]
        sites.loc[site_id, "lon"] = nearest["lon"]
        sites.loc[site_id, "lat"] = nearest["lat"]
    return sites


def snap_to_nearest_safe_point(sites: pd.DataFrame, demand: pd.DataFrame, paths: Paths,
                                boundary: gpd.GeoDataFrame, exclusion_layers: list[dict]) -> pd.DataFrame:
    """A second safety net, run after minimum spacing merges. Merging two
    sites averages their positions the same way a K-Means centroid does, so
    a merged site can land in an excluded zone even when both sites that
    produced it were safe. Merged sites no longer track which original
    demand points they came from, so any site that still fails the
    exclusion check here gets moved to the nearest real demand point
    anywhere in the already filtered demand surface, not just its own
    former cluster."""
    sites = sites.reset_index(drop=True)
    safe_mask = screened_site_mask(sites, paths, boundary, exclusion_layers)
    unsafe_ids = set(np.where(~safe_mask)[0])
    study_area = boundary.geometry.union_all()

    sites = sites.copy()
    for site_id in unsafe_ids:
        site_lon, site_lat = sites.loc[site_id, "lon"], sites.loc[site_id, "lat"]
        inside_boundary_mask = [study_area.covers(Point(lon, lat)) for lon, lat in zip(demand.lon, demand.lat)]
        safe_demand = demand.loc[inside_boundary_mask].reset_index(drop=True)
        if safe_demand.empty:
            continue
        dist = haversine_km(site_lon, site_lat, safe_demand["lon"].to_numpy(), safe_demand["lat"].to_numpy())
        nearest = safe_demand.iloc[int(np.argmin(dist))]
        sites.loc[site_id, "lon"] = nearest["lon"]
        sites.loc[site_id, "lat"] = nearest["lat"]
    return sites


def enforce_minimum_spacing(sites: pd.DataFrame, min_distance_km: float) -> pd.DataFrame:
    """Greedily merge the closest pair of sites, combining their captured
    demand into one demand-weighted centroid, until every remaining pair is
    at least `min_distance_km` apart. A real vertiport cannot share a
    footprint with the one next door, so K-Means centroids that land too
    close together get folded into a single, stronger site instead of
    being reported as two separate sites."""
    df = sites.reset_index(drop=True).copy()
    weight_cols = ["captured_population", "captured_weighted_demand"]
    mean_cols = [c for c in ["mean_transport_hub_km", "mean_major_road_km"] if c in df.columns]

    while len(df) > 1:
        lons, lats = df["lon"].to_numpy(), df["lat"].to_numpy()
        dist = haversine_km(lons[:, None], lats[:, None], lons[None, :], lats[None, :])
        np.fill_diagonal(dist, np.inf)
        i, j = np.unravel_index(np.argmin(dist), dist.shape)
        if dist[i, j] >= min_distance_km:
            break

        a, b = df.iloc[i], df.iloc[j]
        w = a["captured_weighted_demand"] + b["captured_weighted_demand"]
        merged = {
            "lon": (a["lon"] * a["captured_weighted_demand"] + b["lon"] * b["captured_weighted_demand"]) / w,
            "lat": (a["lat"] * a["captured_weighted_demand"] + b["lat"] * b["captured_weighted_demand"]) / w,
        }
        for col in weight_cols:
            merged[col] = a[col] + b[col]
        pop_total = a["captured_population"] + b["captured_population"]
        for col in mean_cols:
            merged[col] = (a[col] * a["captured_population"] + b[col] * b["captured_population"]) / pop_total

        df = df.drop(df.index[[i, j]])
        df = pd.concat([df, pd.DataFrame([merged])], ignore_index=True)

    return df.reset_index(drop=True)


def rank_sites(sites: pd.DataFrame) -> pd.DataFrame:
    ranked = sites.sort_values("captured_weighted_demand", ascending=False).reset_index(drop=True)
    ranked.insert(0, "rank", range(1, len(ranked) + 1))
    return ranked


def add_admin_area_names(sites: pd.DataFrame, boundary: gpd.GeoDataFrame) -> pd.DataFrame:
    """Attach the administrative area containing each proposed coordinate."""
    name_col = next((c for c in ("adm2_name", "name", "LGA") if c in boundary.columns), None)
    if name_col is None:
        return sites
    points = gpd.GeoDataFrame(
        sites.copy(),
        geometry=[Point(lon, lat) for lon, lat in zip(sites["lon"], sites["lat"])],
        crs="EPSG:4326",
    )
    areas = boundary[[name_col, "geometry"]]
    joined = gpd.sjoin(points, areas, how="left", predicate="within")
    result = sites.copy()
    result.insert(3, "local_government_area", joined[name_col].fillna("Outside named area").to_numpy())
    return result


def evaluate_clusters(df: pd.DataFrame, labels: np.ndarray, metric_crs, seed: int, sample_size: int = 5000) -> dict:
    X = to_metric_xy(df, metric_crs)
    sample = min(sample_size, len(X))
    return {
        "silhouette": float(silhouette_score(X, labels, sample_size=sample, random_state=seed)),
        "davies_bouldin": float(davies_bouldin_score(X, labels)),
        "n_demand_points": int(len(X)),
    }


def plot_study_area_and_demand(paths: Paths, config: dict, boundary: gpd.GeoDataFrame, demand: pd.DataFrame) -> Path:
    """Two grounding figures shown before any clustering happens: the plain
    study area boundary, and the raw population weighted demand surface with
    no cluster colors or site markers yet."""
    import matplotlib.pyplot as plt

    fig, axes = plt.subplots(1, 2, figsize=(14, 6))

    boundary.plot(ax=axes[0], color="lightgrey", edgecolor="black")
    axes[0].set_title(f"Study area: {config['city_name']}")
    axes[0].set_xlabel("Longitude")
    axes[0].set_ylabel("Latitude")

    boundary.boundary.plot(ax=axes[1], color="black", linewidth=0.6)
    scatter = axes[1].scatter(demand["lon"], demand["lat"], c=demand["population"], cmap="inferno_r", s=2, alpha=0.6)
    fig.colorbar(scatter, ax=axes[1], label="Estimated population per pixel")
    axes[1].set_title("Raw population weighted demand surface")
    axes[1].set_xlabel("Longitude")

    fig.tight_layout()
    out_path = paths.media / "study_area_and_demand.png"
    fig.savefig(out_path, dpi=150, bbox_inches="tight")
    plt.close(fig)
    return out_path


def plot_static_map(paths: Paths, config: dict, boundary: gpd.GeoDataFrame, df: pd.DataFrame,
                     labels: np.ndarray, sites: pd.DataFrame, buffers_by_layer: dict) -> Path:
    import matplotlib.pyplot as plt

    fig, ax = plt.subplots(figsize=(11, 11))
    boundary.boundary.plot(ax=ax, color="black", linewidth=0.6)
    ax.scatter(df["lon"], df["lat"], c=labels, cmap="tab20", s=1, alpha=0.35)

    colors = ["crimson", "darkorange", "purple", "teal"]
    for color, (name, buffers) in zip(colors, buffers_by_layer.items()):
        buffers.boundary.plot(ax=ax, color=color, linestyle="--", linewidth=1.1, label=f"{name} buffer")

    ax.scatter(sites["lon"], sites["lat"], c="red", marker="x", s=45, label="Ranked proposed site", zorder=5)
    for _, site in sites.head(5).iterrows():
        ax.annotate(
            f"#{int(site['rank'])}",
            (site["lon"], site["lat"]),
            xytext=(5, 5),
            textcoords="offset points",
            fontsize=9,
            fontweight="bold",
            color="black",
            bbox={"boxstyle": "round,pad=0.15", "facecolor": "white", "alpha": 0.8, "edgecolor": "none"},
        )
    ax.set_title(f"Proposed vertiport sites for {config['city_name']}nranked by weighted demand, with exclusion buffers")
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    ax.legend(loc="lower left", fontsize=8)
    out_path = paths.media / "vertiport_sites_static.png"
    fig.savefig(out_path, dpi=150, bbox_inches="tight")
    plt.close(fig)
    return out_path


def plot_interactive_map(paths: Paths, config: dict, boundary: gpd.GeoDataFrame,
                          sites: pd.DataFrame, buffers_by_layer: dict) -> Path:
    import folium

    minx, miny, maxx, maxy = boundary.total_bounds
    fmap = folium.Map(location=[(miny + maxy) / 2, (minx + maxx) / 2], zoom_start=10, tiles="cartodbpositron")

    id_cols = [c for c in boundary.columns if "name" in c.lower()]
    boundary_outline = boundary[[id_cols[0], "geometry"]] if id_cols else boundary[["geometry"]]
    folium.GeoJson(boundary_outline.to_json(), name="City boundary",
                    style_function=lambda _: {"color": "black", "weight": 1, "fillOpacity": 0}).add_to(fmap)

    for _, site in sites.iterrows():
        rank = int(site["rank"])
        folium.CircleMarker(
            location=[site["lat"], site["lon"]], radius=6, color="crimson", fill=True, fill_opacity=0.9,
            popup=f"Rank #{rank} | population captured: {site['captured_population']:,.0f}",
        ).add_to(fmap)
        folium.map.Marker(
            [site["lat"], site["lon"]],
            icon=folium.DivIcon(
                icon_size=(24, 24), icon_anchor=(12, 8),
                html=f'
{rank}
',
            ),
        ).add_to(fmap)

    colors = ["crimson", "darkorange", "purple", "teal"]
    for color, (name, buffers) in zip(colors, buffers_by_layer.items()):
        folium.GeoJson(buffers.to_json(), name=f"{name} buffer",
                        style_function=lambda _, c=color: {"color": c, "dashArray": "4", "fillOpacity": 0.05}).add_to(fmap)

    folium.LayerControl().add_to(fmap)
    out_path = paths.outputs / "vertiport_sites_map.html"
    fmap.save(str(out_path))
    return out_path


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument("--config", type=Path, default=Path(__file__).parent / "configs" / "lagos.json")
    parser.add_argument("--k", type=int, default=None, help="Override number_of_sites from the config")
    parser.add_argument("--min-distance-km", type=float, default=None, help="Override minimum_site_distance_km from the config")
    parser.add_argument("--seed", type=int, default=SEED)
    args = parser.parse_args()

    config = load_config(args.config)
    k = args.k or config["number_of_sites"]
    min_distance_km = args.min_distance_km if args.min_distance_km is not None else config["minimum_site_distance_km"]

    article_root = Path(__file__).resolve().parents[1]
    paths = Paths.from_root(article_root, slugify_city(config["city_name"]))

    print(f"City: {config['city_name']}")
    print("Loading city boundary...")
    boundary = load_city_boundary(paths, config)
    print(f"  bounds={boundary.total_bounds}")

    print("Building population-weighted demand surface...")
    demand = load_population_points(paths, config, boundary)
    n_start = len(demand)
    print(f"  {n_start:,} demand points, total population ~{demand['population'].sum():,.0f}")

    study_area_path = plot_study_area_and_demand(paths, config, boundary, demand)
    print(f"  Saved study area and demand surface figure to {study_area_path}")

    print("Applying exclusion layers...")
    demand, buffers_by_layer, exclusion_report = apply_exclusion_layers(demand, paths, boundary, config["exclusion_layers"])
    print(f"  {n_start - len(demand):,} demand points removed: {exclusion_report}")

    print("Scoring transport and road accessibility...")
    demand = add_accessibility_features(demand, paths, config)

    demand["weighted_demand"] = compute_weighted_demand(demand, config["objective_weights"])

    metric_crs = boundary.estimate_utm_crs()

    print(f"Fitting weighted K-Means with K={k}...")
    labels, centers = site_vertiports(demand, k, args.seed, metric_crs)

    print("Evaluating cluster quality...")
    metrics = evaluate_clusters(demand, labels, metric_crs, args.seed)
    print(f"  Silhouette: {metrics['silhouette']:.4f}  |  Davies-Bouldin: {metrics['davies_bouldin']:.4f}")

    sites = summarize_clusters(demand, labels, centers)

    print("Checking that every site clears the exclusion layers on its own coordinate...")
    sites = snap_unsafe_sites(sites, demand, labels, paths, boundary, config["exclusion_layers"])

    print(f"Enforcing minimum site spacing of {min_distance_km} km...")
    sites = enforce_minimum_spacing(sites, min_distance_km)
    print(f"  {k} requested sites merged down to {len(sites)} after spacing enforcement")

    print("Rechecking sites after spacing merges, since a merged centroid can land in an excluded zone...")
    sites = snap_to_nearest_safe_point(sites, demand, paths, boundary, config["exclusion_layers"])

    sites = rank_sites(sites)
    sites = add_admin_area_names(sites, boundary)

    final_unsafe_count = int((~screened_site_mask(
        sites, paths, boundary, config["exclusion_layers"]
    )).sum())
    if final_unsafe_count:
        print(f"  Warning: {final_unsafe_count} final site(s) still fail the exclusion check after spacing merges")

    metrics_path = paths.outputs / "vertiport_metrics.json"
    metrics_path.write_text(json.dumps({
        "city_name": config["city_name"],
        "requested_sites": k,
        "final_sites_after_spacing": len(sites),
        "final_sites_unsafe_count": final_unsafe_count,
        "minimum_site_distance_km": min_distance_km,
        "objective_weights": config["objective_weights"],
        "seed": args.seed,
        "demand_points_before_exclusions": n_start,
        "demand_points_after_exclusions": len(demand),
        "exclusion_report": exclusion_report,
        **metrics,
    }, indent=2))
    print(f"  Saved metrics to {metrics_path}")

    sites_csv = paths.outputs / "vertiport_sites.csv"
    sites.to_csv(sites_csv, index=False)
    print(f"  Saved {len(sites)} ranked sites to {sites_csv}")

    static_path = plot_static_map(paths, config, boundary, demand, labels, sites, buffers_by_layer)
    print(f"  Saved static map to {static_path}")

    interactive_path = plot_interactive_map(paths, config, boundary, sites, buffers_by_layer)
    print(f"  Saved interactive map to {interactive_path}")


if __name__ == "__main__":
    main()

Setup and how to run this

Python 3.11 or later is required. Clone the companion repository, then from its root folder:

pip install -r code/requirements.txt
python code/vertiport_siting.py --config code/configs/lagos.json

The requirements.txt file pins the exact package versions used to generate the results. These are the map, table, and machine learning libraries listed in the file. No API key is needed for this pipeline.

The boundary and road files use common map formats such as GeoJSON and shapefiles. You do not need to edit them by hand; replace them with equivalent files for your city and point to them from the settings file.

The boundary, population, road, and transit files for both cities are bundled directly in the data/ folder of the companion repository, listed with their exact source links, license, and a SHA-256 checksum in data/README.md, so the first run does not depend on any dataset provider staying online. The one exception is the ESA WorldCover land cover raster: that file is not bundled, and vertiport_siting.py fetches and clips it automatically from Microsoft Planetary Computer, a public, unauthenticated endpoint, the first time it runs for a given city, then caches it in data/ for every run after that.

Expected folder layout, all produced automatically on first run:

articles/lagos-vertiport-siting/
  code/
    vertiport_siting.py
    compare_variants.py
    requirements.txt
    configs/
      lagos.json
      nairobi.json
  data/
    lagos_lgas.geojson
    lagos_worldpop_gridded_v2_0.tif
    lagos_major_roads.geojson
    lagos_transport_hubs.geojson
    lagos_worldcover_2021.tif        (built automatically on first run)
    nairobi_subcounties.geojson
    nairobi_worldpop_gridded_v2_0.tif
    nairobi_major_roads.geojson
    nairobi_transport_hubs.geojson
  outputs/
    lagos/vertiport_sites.csv, vertiport_metrics.json, vertiport_sites_map.html
    nairobi/vertiport_sites.csv, vertiport_metrics.json, vertiport_sites_map.html
  media/
    lagos/study_area_and_demand.png, vertiport_sites_static.png, recommended_25_sites.png, naive_vs_constrained.png, sensitivity_sweep.png, site_count_analysis.png
    nairobi/study_area_and_demand.png, vertiport_sites_static.png

Every run in this article used a fixed random seed of 42, so the site coordinates, silhouette scores, and ranking reported here should reproduce exactly on a rerun with the same settings and package versions. The three version comparison, sensitivity sweep, elbow chart, and site count table all live in code/compare_variants.py, a second script that imports and reuses every function shown above rather than duplicating the method. Running python code/compare_variants.py --config code/configs/lagos.json regenerates the comparison table, both sensitivity charts, and the site count analysis used above.

What actually transfers to your own city

The reusable lesson here is that K means cannot choose infrastructure sites on its own, and the first population only version proved that concretely: 11 proposed sites failed the screening check when population was the only input.

The same distinction showed up in an earlier project fine tuning a robot control model on a limited Colab graphics processing unit (GPU): a short training run proves the training path is wired correctly, not that a robot is ready to deploy. This pipeline is the same kind of evidence. It proves the siting method is wired correctly, from real population data through exclusion checks to a ranked list. It does not prove Lagos or Nairobi should build vertiports at these exact coordinates; that would need the routed travel time, verified airspace, and flood risk data this article already flagged as missing.

What actually transfers is the combination: population used openly as a stand in for demand, clearly labeled as an estimate rather than real commuting data; accessibility scoring against transit and roads that a reader can inspect line by line; exclusion layers built from named, verifiable sources rather than assumed constraints; and a site count analysis that connects the number of locations to a stated coverage target. Together, those four components turn incomplete city data into a defensible first pass planning tool, not a final answer.

Three findings mattered most in this specific run. Adding exclusion layers changed the safety outcome completely, from 11 unsafe sites to zero. Adding accessibility weighting changed the ranking meaningfully but did not fix safety on its own, so exclusions have to run regardless of what else the model considers, no matter what order the other steps happen in. And the airport buffer radius that worked fine for Lagos removed three times the proportional demand in Nairobi, which means a constraint tuned for one city is not automatically calibrated for the next one.

If you are applying this to your own city, replace the boundary, population, roads, transit stops, and screening layers first. Then set your own coverage target in site_count_policy, run the site count analysis, and use the smallest tested number that meets that target. Run the sensitivity sweep as well, before trusting the first map the method gives you.

What would actually change your mind about where the buffer radius should come from: a fixed safety standard, or a number recalculated for every city’s specific airport geometry? That question does not have a clean answer yet, and it is exactly the kind of open question worth arguing about in the comments.

Glossary

Urban Air Mobility (UAM): the general term for using small electric aircraft, including air taxis, for short trips within or around a city, instead of ground transport.

Vertiport: a small airfield built for aircraft that take off and land straight up and down, like electric air taxis, the way a helipad works for helicopters, but sized and equipped for electric aircraft.

K means: an algorithm that groups a set of points into a fixed number of clusters, then reports the center point of each cluster. In this article, the points being clustered are population weighted locations, and the resulting centers become proposed vertiport sites.

Weighted K means: the same algorithm, but each point counts more or less toward where a cluster’s center lands, based on a weight value. Here, that weight combines population with transit and road accessibility.

Silhouette score: a number between negative 1 and 1 that measures how well separated a set of clusters are. Higher is better. A score in the 0.3 to 0.4 range, as in both this article’s runs, indicates moderate, usable separation rather than sharply distinct clusters.

Davies Bouldin index: a second cluster quality measure, where lower is better, based on how compact each cluster is relative to how far apart clusters are from each other.

Stand in for demand: a measurable substitute used in place of something that cannot be directly measured. Population density stands in for travel demand here, since no origin destination commuting dataset exists for either city in this article.

Coordinate reference system (CRS): a defined way of mapping locations on the earth to numerical coordinates. EPSG:4326 is the CRS used for standard GPS latitude and longitude; distance and area calculations in this pipeline temporarily switch to a local metric CRS so that buffer radii can be measured in real kilometers.

Cloud optimized GeoTIFF: a version of the GeoTIFF raster image format organized so a program can read a specific, bounded region of a very large file over the internet without downloading the whole thing. This is how the pipeline pulls only Lagos’s or Nairobi’s slice of the global ESA WorldCover dataset without downloading the entire multi gigabyte tile.

Exclusion layer: a rule that removes proposed locations from consideration entirely, regardless of how much population or accessibility they would otherwise score. This article uses two kinds: a radius around a named point, and a land cover classification.

Proposed site: a location the pipeline suggests for further review, not a final decision. Every ranked site in this article’s output is a proposed site, not a confirmed vertiport location.

Selected Sources

• Nigeria Civil Aviation Authority, Nigeria Civil Aviation Regulations Part 14: official regulatory source for the identification of prohibited, restricted, and danger airspace areas referenced in the opening and the exclusion layer design.
• Proceedings of the International Association of Hydrological Sciences, Urban flood vulnerability mapping of part of the Lagos metropolis (2020): source for the Lagos flood risk figures cited in the datasets and limitations sections.
• Humanitarian Data Exchange, Nigeria Subnational Administrative Boundaries and Kenya Subnational Administrative Boundaries: source for both cities’ Local Government Area and sub county boundary files, licensed CC BY IGO.
• GRID3, Bottom up gridded population estimates for Nigeria, version 2.0 and the equivalent Kenya gridded population estimates, version 2.0: source for the population raster used as the stand in for demand in both cities, licensed CC BY.
• European Space Agency, WorldCover 2021: source for the land cover classification used to exclude water, wetland, and mangrove, accessed through Microsoft Planetary Computer.
• OpenStreetMap contributors, via the Overpass API and Nominatim: source for verified airport, military installation, transit stop, and major road coordinates in both cities, licensed under the Open Database License.
• scikit-learn documentation, KMeans: reference for the weighted K means implementation used throughout this pipeline.
• I Tried Fine Tuning a Robot AI Model on Colab, Here Is What Worked: source for the smoke test versus deployment readiness distinction applied to this pipeline’s own conclusion.