Example Network Gravity Calculations

Example Network Gravity Calculations#

This notebook demonstrates how to calculate network gravity using the gravity model. The gravity model is a simple model that estimates the flow of people, goods, or information between two locations. The model is based on the idea that the flow between two locations is proportional to the product of their masses and inversely proportional to the distance between them.

Import Libraries and Load Data#

Import Libraries#

We import the following libraries:

  • os to define paths

  • rasterio to read raster files

  • geopandas to read vector files

  • pandas to manipulate data

  • pickle to load the network data

  • GOSTnets to access custom GOST functions

import os
import rasterio

import geopandas as gpd
import pandas as pd
import pickle as pkl

# import the GOSTnets library
import GOSTnets as gn

GDAL is not installed - OGR functionality not available

Load Data#

Define the path information and load the data.

# Define input data
pth = "./"
# Read in cleaned pickle from earlier analysis and convert to time
G = pkl.load(
    open(
        os.path.join(pth, "tutorial_outputs", r"iceland_network_clean.pickle"),
        "rb",
    )
)
G_time = gn.convert_network_to_time(
    G, distance_tag="length", road_col="infra_type", factor=1000
)
# Define origins and destinations
rek_grid_file = os.path.join(pth, "tutorial_data", "rek_grid.shp")
rek_pop_grid_file = rek_grid_file.replace(".shp", "_pop.shp")
churches_file = os.path.join(pth, "tutorial_data", "churches.shp")

Gravity calculations#

When we discuss gravity calculations, we are combining two measures to identify the attractiveness between pairs of locations

  1. Travel time between origins and destinations

  2. Weight of the value of interest at both the origin and destination

‘Market access’ is a measure of accessibility from one origin to all destinations, based purely on time. The outcome of this analysis is often visualized as an ‘isochrone’, centered on a given origin point. Isochrones depict how far away from the origin you can move, in all directions, within a certain time frame. (E.g. 1 hour).

‘Market potential’ can be assessed by weighting all potential destinations by a factor designed to represent their attractiveness to the origin point. In the context of performing commercial site selection, such factors might include GDP per capita or population (positively related), or the number of other firms (negatively related).

# Sample rek_grid with population - the data should already be in the tutorial inputs, so this attribution is not necessary
if not os.path.exists(rek_pop_grid_file):
    population_data = "R:\GLOBAL\POP&DEMO\LandScan_2012\lspop2012.tif"

    in_grid = gpd.read_file(rek_grid_file)
    in_pop = rasterio.open(population_data)
    in_grid = in_grid.to_crs(in_pop.crs)

    geoms = ([x.x, x.y] for x in in_grid["geometry"])
    vals = in_pop.sample(geoms)
    in_grid["Pop"] = [x[0] + 1 for x in vals]

    in_grid.to_file(rek_pop_grid_file)

Read the data and calculate origin-destination table#

We will use the GOSTnets library to calculate the origin-destination table. This table will be used to calculate the gravity model.

# Read in the input
rek_grid = gpd.read_file(rek_pop_grid_file)
in_churches = gpd.read_file(churches_file)
in_churches = in_churches.to_crs(rek_grid.crs)

# calculate the origins and destinations by snapping to the road network
origins_df = gn.pandana_snap_c(
    G_time,
    rek_grid,
    source_crs="epsg:4326",
    target_crs="epsg:4326",
    add_dist_to_node_col=True,
)
origins = list(set(origins_df["NN"]))
destinations_df = gn.pandana_snap_c(
    G_time,
    in_churches,
    source_crs="epsg:4326",
    target_crs="epsg:4326",
    add_dist_to_node_col=True,
)
destinations = list(set(destinations_df["NN"]))

OD = gn.calculate_OD(G_time, origins, destinations, fail_value=9999999)

Calculate the origin-destination table and gravity#

We will calculate the origin-destination table and gravity using the GOSTnets library.

# The calculate_od_raw contains functions for performing OD and gravity calculations
import GOSTnets.calculate_od_raw as calcOD

# For this calculation the origins are all weighted equally with a value of 1
gravity = calcOD.calculate_gravity(
    OD, oWeight=[1] * OD.shape[0], dWeight=destinations_df["Pop"]
)
gravity["NN"] = origins

# merge the gravity results with the origins
origins_gravity = pd.merge(origins_df, gravity, on="NN")
origins_gravity.to_file(rek_pop_grid_file.replace(".shp", "_gravity.shp"))

/tmp/ipykernel_2025/4100008754.py:3: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.
  origins_gravity.to_file(rek_pop_grid_file.replace(".shp", "_gravity.shp"))

# Now let's introduce the population numbers into the origins
origins = origins_df["NN"]
destinations = destinations_df["NN"]
OD = gn.calculate_OD(G_time, origins, destinations, fail_value=9999999)
gravity = calcOD.calculate_gravity(
    OD, oWeight=origins_df["Pop"], dWeight=destinations_df["Pop"]
)
gravity["NN"] = origins

# merge the gravity results with the origins
origins_gravity = pd.merge(origins_df, gravity, on="NN")
origins_gravity.to_file(rek_pop_grid_file.replace(".shp", "_gravity2.shp"))

/tmp/ipykernel_2025/1925333871.py:3: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.
  origins_gravity.to_file(rek_pop_grid_file.replace(".shp", "_gravity2.shp"))

Examine the data#

You can explore and examine the data by opening it up and loading it in QGIS.

Contents