Changing the zone system for a model

In this example we show how to change the zone system for a Polaris model.

from pathlib import Path
import shutil
from tempfile import gettempdir

import geopandas as gpd
import numpy as np
import pandas as pd
from shapely.geometry import Polygon, MultiPolygon

from polaris import Polaris
from polaris.network.network import Network
from polaris.project.project_restorer import restore_project_from_csv
from polaris.utils.database.db_utils import commit_and_close

original_model = Path("/tmp/Bloomington")
target = original_model.parent / "bloomington_new_zones_from_census_block"

shutil.copytree(original_model, target, dirs_exist_ok=True)

tract_path = Path(gettempdir()) / "illinois_census_block_groups_portion.gpkg"

model = Polaris.from_dir(original_model)
model.upgrade()
from_net = model.network
db_name = model.run_config.db_name

zones = from_net.data_tables.get_geo_layer("zone")

Read the tracts and put them in the same CRS as the zones

tracts = gpd.read_file(tract_path)
tracts = tracts.to_crs(zones.crs)

# We will use GEOID because we know it is an unique identifier, but we could use any unique field in this table.
id_field = "GEOID"

Compute the overlay between the zone and the block group geometries

ovrl = zones.overlay(tracts, how="intersection")[["zone", id_field, "geometry"]]

We compute the area for each of the overlaps We will use the largest overlap for each new zone to copy the data for a few fields from the original zones Most fields will be aggregated with a sum, however

ovrl = ovrl.assign(overlay_area=ovrl.area)
idx = ovrl.groupby("zone")["overlay_area"].idxmax()
df_max = ovrl.loc[idx].reset_index(drop=True)

zones_with_tract = zones.merge(df_max[["zone", id_field]], on="zone", how="left")

copy_fields = ["area_type", "electric_grid_transmission", "electricity_provider"]
data = zones_with_tract[copy_fields + [id_field]].set_index(id_field)

data = data[~data.index.duplicated(keep="first")]

We will dissolve the zones by the new id_field and sum the data for all fields

new_zones = zones_with_tract.dissolve(by=id_field, aggfunc="sum")

# We drop the fields that we should not have summed
new_zones.drop(columns=copy_fields, inplace=True)

# And copy the data from the original zones, based on the largest overlap, as defined above
new_zones = new_zones.join(data)
new_zones.loc[:, "zone"] = np.arange(1, new_zones.shape[0] + 1)
new_zones = new_zones[zones.columns].reset_index(drop=True)

We will save the new zones to a CSV file For that, we save the geometry as WKT in the appropriate format, drop the geometry column And make sure that the data is in the appropriate schema

new_zones.geo = new_zones.geo.apply(lambda geom: MultiPolygon([geom]) if isinstance(geom, Polygon) else geom)
zone_data = new_zones.assign(geo_wkt=new_zones.geo.to_wkt())
zone_data.drop(columns=["geo"], inplace=True)

with commit_and_close(from_net.path_to_file, spatial=True) as conn:
    schema = pd.read_sql_query("pragma table_info(zone)", conn)
    fields = ["geo_wkt" if fld == "geo" else fld for fld in schema["name"].to_list()]

zone_data[fields].to_csv(target / "supply" / "zone.csv", index=False)

We rebuild the new model

restore_project_from_csv(target, target, db_name, True)

new_net = Network.from_file(target / from_net.path_to_file.name)
new_net.data_tables.get_geo_layer("zone").plot()
# It is HIGHLY advisable to inspect the result of this procedure, as
# undesirable artifacts from the spatial operations may have been introduced

<Axes: >

and update geo_consistency –> THIS IS A CRITICAL STEP

new_net.geo_consistency.update_all()

/builds/polaris/code/polarislib/polaris/network/data/data_table_cache.py:42: UserWarning: THIS TABLE WAS CACHED, USE WITH CAUTION
  warnings.warn("THIS TABLE WAS CACHED, USE WITH CAUTION")