运筹系列56：python空间分析库pysal.spaghtti

1. 介绍

Pysal是一个面向地理空间数据科学的开源跨平台库，重点是用python编写的地理空间矢量数据。它支持空间分析高级应用程序的开发，例如

• explore -用于对空间和时空数据进行探索性分析的模块，包括对点、网络和多边形格的统计测试。还包括空间不等式和分布动力学的方法。
• viz -可视化空间数据中的模式，以检测集群、异常值和热点。
• model -用各种线性、广义线性、广义加性和非线性模型对数据中的空间关系进行建模。
• lib -解决各种各样的计算几何问题：从多边形格、线和点构建图形；空间权重矩阵与图形的构建与交互编辑；α形状、空间指数和空间拓扑关系的计算；读写稀疏图形数据，以及纯python空间矢量数据阅读器。

SPAtial GrapHs: nETworks, Topology, & Inference是pysal下的一个库，主要用于绘图。

如果在mac下安装，要执行如下操作：

brew install spatialindex
sudo pip3 install osmnx

1. 基本绘图案例

1.1 基本路网

import geopandas
import libpysal
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
import spaghetti
%matplotlib inline
ntw = spaghetti.Network(in_data=libpysal.examples.get_path("streets.shp"))
vertices_df, arcs_df = spaghetti.element_as_gdf(ntw, vertices=True, arcs=True)
base = arcs_df.plot(color="k", alpha=0.25, figsize=(12, 12), zorder=0)
# create vertices keyword arguments for matplotlib
vertices_df.plot(color="k", markersize=5, alpha=1, ax = base)

出现下图：

1.2 点吸附

使用snapobservations方法

ntw.snapobservations( libpysal.examples.get_path("schools.shp"), "schools", attribute=False)
print("observation 1\ntrue coords:\t%s\nsnapped coords:\t%s" % (
    ntw.pointpatterns["schools"].points[0]["coordinates"],
    ntw.pointpatterns["schools"].snapped_coordinates[0]
))

返回

observation 1
true coords:    (727082.0462136, 879863.260705768)
snapped coords: (727287.6644417326, 879867.3863186113)

1.3 绘图

vertices_df, arcs_df = spaghetti.element_as_gdf(ntw, vertices=True, arcs=True)
#args = [], []
kwargs = {"c":"k","lw":0}

# set legend entry
arcs = mlines.Line2D(*args, label="Network Arcs", alpha=0.5)
vertices = mlines.Line2D( *args, **kwargs, ms=2.5, marker="o", label="Network Vertices")
tschools = mlines.Line2D(*args, **kwargs, ms=25, marker="X", label="School Locations")
sschools = mlines.Line2D(*args, **kwargs, ms=12, marker="o", label="Snapped Schools")

# draw map
base = arcs_df.plot(color="k", alpha=0.25, figsize=(12, 12), zorder=0)
vertices_df.plot(color="k", markersize=5, alpha=1, ax = base)
kwargs.update({"cmap":"tab20", "column":"id", "zorder":2})
true_schools_df.plot(marker="X", markersize=500,ax = base, **kwargs)
snapped_schools_df.plot(markersize=200, ax = base, **kwargs)
plt.legend(handles=(arcs,vertices,tschools,sschools),
    fancybox=True,
    framealpha=0.8,
    scatterpoints=1,
    fontsize="xx-large",
    bbox_to_anchor=(1.04, 0.75),
    borderpad=2.,
    labelspacing=2.)

2. 和求解器结合：TSP问题

• TSP问题的案例是用pulp求解的
• 选址问题更复杂一点，使用ortools+cbc求解，可以参考：https://pysal.org/spaghetti/notebooks/facility-location.html
• 还有一个运输问题的例子，使用python-mip进行求解，可以参考：https://pysal.org/spaghetti/notebooks/transportation-problem.html

2.1 定义TSP

import geopandas
from libpysal import examples
import matplotlib
import numpy
import pulp
import spaghetti
class MTZ_TSP:
    def __init__(self, nodes, cij, xij_tag="x_%s,%s", ui_tag="u_%s", display=True):
        self.nodes, self.cij = nodes, cij
        self.p = self.nodes.shape[0]
        self.rp_0n, self.rp_1n = range(self.p), range(1, self.p)
        self.id = self.nodes.name
        self.xij_tag, self.ui_tag = xij_tag, ui_tag

        # instantiate a model
        self.tsp = pulp.LpProblem("MTZ_TSP", pulp.LpMinimize)
        # create and set the tour decision variables
        self.tour_dvs()
        # create and set the arbitraty real number decision variables
        self.arn_dvs()
        # set the objective function
        self.objective_func()
        # node entry constraints
        self.entry_exit_constrs(entry=True)
        # node exit constraints
        self.entry_exit_constrs(entry=False)
        # subtour prevention constraints
        self.prevent_subtours()
        # solve
        self.tsp.solve()
        # origin-destination lookup
        self.get_decisions(display=display)
        # extract the sequence of nodes to construct the optimal tour
        self.construct_tour()

    def tour_dvs(self):
        def _name(_x):
            return self.nodes[_x].split("_")[-1]

        xij = numpy.array([[pulp.LpVariable(self.xij_tag % (_name(i), _name(j)), cat="Binary") for j in self.rp_0n]
                for i in self.rp_0n])
        self.xij = xij

    def arn_dvs(self):
        """Create arbitrary real number decision variables - eq (6)."""
        ui = numpy.array([pulp.LpVariable(self.ui_tag % (i), lowBound=0) for i in self.rp_0n])
        self.ui = ui

    def objective_func(self):
        """Add the objective function - eq (1)."""
        self.tsp += pulp.lpSum([self.cij[i, j] * self.xij[i, j]
                for i in self.rp_0n for j in self.rp_0n if i != j])

    def entry_exit_constrs(self, entry=True):
        """Add entry and exit constraints - eq (2) and (3)."""
        if entry:
            for i in self.rp_0n:
                self.tsp += (pulp.lpSum([self.xij[i, j] for j in self.rp_0n if i != j]) == 1)
        # exit constraints
        else:
            for j in self.rp_0n:
                self.tsp += (pulp.lpSum([self.xij[i, j] for i in self.rp_0n if i != j]) == 1)

    def prevent_subtours(self):
        """Add subtour prevention constraints - eq (4)."""
        for i in self.rp_1n:
            for j in self.rp_1n:
                if i != j:
                    self.tsp += (self.ui[i] - self.ui[j] + self.p * self.xij[i, j] <= self.p - 1)

    def get_decisions(self, display=True):
        """Fetch the selected decision variables."""
        cycle_ods = {}
        for var in self.tsp.variables():
            if var.name.startswith(self.ui_tag[0]):
                continue
            if var.varValue > 0:
                if display:
                    print("%s: %s" % (var.name, var.varValue))
                od = var.name.split("_")[-1]
                o, d = [int(tf) for tf in od.split(",")]
                cycle_ods[o] = d
        if display:
            print("Status: %s" % pulp.LpStatus[self.tsp.status])

        self.cycle_ods = cycle_ods

    def construct_tour(self):
        """Construct the tour."""
        tour_pairs = []
        for origin in self.rp_0n:
            tour_pairs.append([])
            try:
                tour_pairs[origin].append(next_origin)
                next_origin = self.cycle_ods[next_origin]
                tour_pairs[origin].append(next_origin)
            except NameError:
                next_origin = self.cycle_ods[origin]
                tour_pairs[origin].append(origin)
                tour_pairs[origin].append(next_origin)

        tour_pairs = {idx: sorted(tp) for idx, tp in enumerate(tour_pairs)}
        self.tour_pairs = tour_pairs

    def extract_tour(self, paths, id_col, leg_label="leg"):

        paths[leg_label] = int
        # set label of journey leg for each OD pair.
        for leg, tp in self.tour_pairs.items():
            paths.loc[paths[id_col] == tuple(tp), leg_label] = leg

        # extract only paths in the tour
        tour = paths[paths[leg_label] != int].copy()
        tour.sort_values(by=[leg_label], inplace=True)

        return tour

2.2 定义地图数据

streets = geopandas.read_file(examples.get_path("streets.shp"))
streets.crs = "esri:102649"
streets = streets.to_crs("epsg:2762")
all_crimes = geopandas.read_file(examples.get_path("crimes.shp"))
all_crimes.crs = "esri:102649"
all_crimes = all_crimes.to_crs("epsg:2762")
numpy.random.seed(1960)
koenigsberg_cases = 7 * 2
subset_idx = numpy.random.choice(all_crimes.index, koenigsberg_cases, replace=False)
crimes_scenes = all_crimes[all_crimes.index.isin(subset_idx)].copy()
ntw = spaghetti.Network(in_data=streets)
ntw.snapobservations(crimes_scenes, "crime_scenes")
pp_obs = spaghetti.element_as_gdf(ntw, pp_name="crime_scenes")
pp_obs_snapped = spaghetti.element_as_gdf(ntw, pp_name="crime_scenes", snapped=True)

2.3 求解，并提取最短路

d2d_dist, tree = ntw.allneighbordistances("crime_scenes", gen_tree=True)
pp_obs["dv"] = pp_obs["id"].apply(lambda _id: "x_%s" % _id)
mtz_tsp = MTZ_TSP(pp_obs["dv"], d2d_dist)
paths = ntw.shortest_paths(tree, "crime_scenes")
paths_gdf = spaghetti.element_as_gdf(ntw, routes=paths)
tour = mtz_tsp.extract_tour(paths_gdf.copy(), "id")

求解结果为：

x_0,5: 1.0
x_1,0: 1.0
x_10,12: 1.0
x_11,6: 1.0
x_12,11: 1.0
x_13,9: 1.0
x_2,1: 1.0
x_3,2: 1.0
x_4,3: 1.0
x_5,7: 1.0
x_6,4: 1.0
x_7,13: 1.0
x_8,10: 1.0
x_9,8: 1.0
Status: Optimal

2.4 绘图

def tour_labels(t, b):
    def _lab_loc(_x):return _x.geometry.interpolate(0.5, normalized=True).coords[0]
    kws = {"size": 20, "ha": "center", "va": "bottom", "weight": "bold"}
    t.apply(lambda x: b.annotate(text=x.leg, xy=_lab_loc(x), **kws), axis=1)


def obs_labels(o, b):
    def _lab_loc(_x):return _x.geometry.coords[0]
    kws = {"size": 14, "ha": "left", "va": "bottom", "style": "oblique", "color": "k"}
    o.apply(lambda x: b.annotate(text=x.id, xy=_lab_loc(x), **kws), axis=1)
base = arcs.plot(alpha=0.2, linewidth=1, color="k", figsize=(10, 10), zorder=0)
tour.plot(ax=base, column="leg", cmap="tab20", alpha=0.50, linewidth=10, zorder=1)
vertices.plot(ax=base, markersize=1, color="r", zorder=2)
pp_obs.plot(ax=base, markersize=20, color="k", zorder=3)
pp_obs_snapped.plot(ax=base, markersize=20, color="k", marker="x", zorder=2)
# tour leg labels
tour_labels(tour, base)
# crime scene labels
obs_labels(pp_obs, base)

绘制结果如下：