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, cijself.p = self.nodes.shape[0]self.rp_0n, self.rp_1n = range(self.p), range(1, self.p)self.id = self.nodes.nameself.xij_tag, self.ui_tag = xij_tag, ui_tag# instantiate a modelself.tsp = pulp.LpProblem("MTZ_TSP", pulp.LpMinimize)# create and set the tour decision variablesself.tour_dvs()# create and set the arbitraty real number decision variablesself.arn_dvs()# set the objective functionself.objective_func()# node entry constraintsself.entry_exit_constrs(entry=True)# node exit constraintsself.entry_exit_constrs(entry=False)# subtour prevention constraintsself.prevent_subtours()# solveself.tsp.solve()# origin-destination lookupself.get_decisions(display=display)# extract the sequence of nodes to construct the optimal tourself.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 = xijdef 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 = uidef 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 constraintselse: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]):continueif 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] = dif display:print("Status: %s" % pulp.LpStatus[self.tsp.status])self.cycle_ods = cycle_odsdef 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_pairsdef 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 tourtour = 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)

绘制结果如下：

运筹系列56：python空间分析库pysal.spaghtti相关推荐

python空间分析库_空间分析：5-1.空间分析库PySAL的使用
Pysal与geoda非常相似,一个通过写脚本来实现空间分析,一个通过软件操作来实现空间分析. Pysal的官网对于自己的介绍是,开源.跨平台的地理空间数据分析库. Pysal能干什么? 空间分析+可 ...
开源的前端GIS空间分析库介绍（三）turf与ol结合
前言 turf是mapbox出品的前端空间分析库,官网:http://turfjs.org/ turf库中包含的空间分析计算功能比较多,也非常简单易用.相比于jsts,turf的官方文档维护的非常好, ...
JavaScript 空间分析库——JSTS和Turf
前言项目中有管线的空间拓扑关系查询需求,在npm中检索到JSTS和Turf两个JavaScript 空间分析库. JSTS JSTS是一个符合OGC规范的简单要素空间位置判定函数JavaScript ...
python空间分析_读书笔记——《python地理空间分析指南》
本文为<Python地理空间分析指南(第2版)>的读书摘录,顺便挖个坑,进一步对python的几个包做学习整理. 本笔记的用途:了解python地理空间处理的技术框架和实现途径. 第三章 ...
开源的前端GIS空间分析库介绍（一）jsts与turf
文章目录 1 前言 2 JSTS 3 turf 4 安装使用 4.1 jsts 4.1.1 直接引入 4.1.2 NPM 4.2 turf 4.1.1 直接引入 4.1.2 NPM 5 空间分析 5. ...
Turf.js 地理空间分析库简介
Turf.js是一个轻量级的JavaScript库,用于地理空间分析和操作.它提供了许多强大的函数和算法,用于处理地理空间数据,如点.线.多边形和网格等.Turf.js的API简单易用,可以轻松地与其 ...
python网络爬虫系列教程——python中requests库应用全解
全栈工程师开发手册 (作者:栾鹏) python教程全解 python中requests库的基础应用,网页数据挖掘的常用库之一.也就是说最主要的功能是从网页抓取数据. 使用前需要先联网安装reques ...
python深度学习库系列教程——python调用opencv库教程
分享一个朋友的人工智能教程.零基础!通俗易懂!风趣幽默!还带黄段子!大家可以看看是否对自己有帮助:点击打开全栈工程师开发手册 (作者:栾鹏) python教程全解 OpenCV安装 pip inst ...
python 股票分析库_GitHub - reference-project/stock-1: stock，股票系统。使用python进行开发。...
pythonstock V1 项目简介特别说明:股市有风险投资需谨慎,本项目只能用于Python代码学习,股票分析,投资失败亏钱不负责,不算BUG. PythonStock V1 是基于Python ...

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

1. 介绍

1. 基本绘图案例

1.1 基本路网

1.2 点吸附

1.3 绘图

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

2.1 定义TSP

2.2 定义地图数据

2.3 求解，并提取最短路

2.4 绘图

运筹系列56：python空间分析库pysal.spaghtti相关推荐

最新文章

热门文章