使用Python函数式编程进行空间回归分析

2 0 0 9. 0 6 No. 3 WWW. GWN. ORG. CN地理信息世界 GEOMATICS WORLD 学术研究 Abstract: As the traditional programming languages do not suite for spatial data analysis, functional programming is more advanced and more practical. This paper put forward to intro－ duce Python’s functional programming into spatial regression analysis. It analyzed and dis－ cussed application prospect of script language combined with functional programming method in GIS data analysis. The experiments proved that it can effectively improve the programming efficiency, reduce the codes, and focus the researcher’s attention on the problem using the Python’s functional programming combined with the existed open source software packages. Keywords: Python; functional programming; Spatial Regression Analysis; Open Source Soft－ ware; data mining 0 引 言 从目前涉及空间数据分析的 主要教材、 书 关于书的成语关于读书的排比句社区图书漂流公约怎么写关于读书的小报汉书pdf 籍和 论文 政研论文下载论文大学下载论文大学下载关于长拳的论文浙大论文封面下载 来看，对于 空间分析的实现和解算，一部分研 究人员采用封装了统计功能的 GIS，RS软件模块进行操作 [1]；另一 部分人员采用诸如 SPSS，Matlab等 统计、数学计算软件实现[2]。既有的 软件只包含有限的功能，可定制性 不强；数学软件不易操作 GIS特有 的数据结构，数据准备的过程较为 烦琐。而在学校教学中，学生多数 使用 C语言等传统的编程语言完 成空间分析程序的实现，过程烦 琐。对于较为灵活的探索性空间数 据分析来说，需要不断尝试各种模 型进行比较分析，因此要求编程工 具具备灵活性和强大数值计算能 力。Python是一种面向对象的动态 语言，采用 Python进行函数式编程 恰好可以解决此类问题，在编程效 率和可定制性之间找到一个平衡， 也便于空间分析模块整合在更复 杂的程序中。同时，Python平台支持 从数据库连接到 Web 服务器的全 面支持，以 Python开发的程序模块 可以和系统整合参与实际应用。 1 编程工具简介 Python 是一种底层基于 C语 言实现的脚本语言，诞生于 1991 年。作为一种脚本语言，Python采用 解释运行的方式，编写后无需进行 编译即可直接通过解释器执行，具 有典型的动态语言特点，编程效率 极高。Python[3]目前版本为 2.6，其已 从 Linux平台上的一个脚本语言发 展成为跨操作系统、多种底层实现 的通用语言，并且可以实现对既有 的 C/C++库调用，大大丰富了 Pyt- hon的功能。 函数式编程（Functional Pro- gramming）是一种编程范式[4]。以函 数为程序基础，通过向函数传入参 数，返回输出而获得计算结果，函 数本身由多个函数组成，直至函数 的最底层仅仅采用程序的基本运 算功能实现[5]。采用函数为基础，抛 弃了传统过程式语言的控制流和 状态变量，对于复杂的序列计算， 使用 Python函数式编程进行空间回归分析 文章编号：1672- 1586（2009）02- 0077- 04 中图分类号：TP311.132.4 文献标识码：A 摘要： 针对使用传统的编程方式进行空间回归分析过于繁杂的问题，本文提出将 Python 函数式编程引入空间回归分析方法，分析和探讨脚本语言+函数式编程方式在 GIS数据分 析中的推广应用前景。编程实践证明，采用 Python函数式编程，结合既有的开源软件包，可 以有效提高编程效率，减少代码量，将实践人员的注意力集中在问题本身。 关键词：Python；函数式编程；空间回归分析；开源软件；数据挖掘 孙 宁，蒲英霞 （南京大学 地理与海洋科学学院，江苏 南京 210093） SUN Ning，PU Ying-xia （School of Geographic and Oceanographic Sciences，Nanjing University，Nanjing 210093，China） Spatial Regression Analysis with Python’s Functional Programming 孙 宁（1987-），男， 江苏南京人，本科学生， 南京大学地理与海洋科 学学院地理信息科学系， 研究方向：GIS 开放 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 及开源软件应用，脚本语 言发展及应用。 基金项目：国家自然科学 基金项目资助（40601074）。 E- mail：class icning @163.com 收稿日期：2008- 06- 28 66 2 0 0 9. 0 6 No.3 WWW. GWN. ORG. CN 地理信息世界 GEOMATICS WORLD 学术研究 有效提高了编程效率。实现了函数 式编程的语言包括 Python，Haskell， Scheme，Erlang和各种数学软件专 用语言（如 S，ML等）。 Python函数式编程支持[6]： 1）匿名函数 lambda 在 Python中，可以使用 lambda 关键字定义匿名函数，并将此函数 作为参数传入其他函数中。 lambda argument1，argument2，. .. argumentN : expression using argu- ments 2）map（） 映射函数，可以按照参数中指 定的函数对序列进行映射，生成新 的序列。例如： >>> map（lambda x: x*2，range （5）） [0，2，4，6，8] 3）reduce（） reduce操作通过指定的函数对 序列进行迭代运算。例如： >>> reduce （lambda x，y: x· （y+1），range（9），1） 362880 4）filter（） 过滤器用户从一个序列中取 出满足条件的子序列。例如： >>> filter（lambda x: x%2==0， range（10）） [0，2，4，6，8] 2 空间回归分析 空间数据分析（Spatial Data Analysis，SDA）是 GIS进行分析、决 策的重要工具。利用空间数据分析 技术，探索和挖掘现象分布的空间 关系是 SDA的核心所在。SDA的实 现分为以下四个层次：①认知，对 空间数据进行有效获取和科学的 组织描述，利用空间数据来再现事 物本身；②解释，理解并解释空间 数据的背景过程，认识事件的本质 规律；③预报，了解、掌握事件发生 的规律后，运用预测模型对未来的 状况做出合理推测；④宏观决策和 调控，根据 SDA结果做出合理决 策，调控地理空间上发生的事件[7]。 空间数据分析通过对已存在的现象 匹配最为接近的地理空间模型，建 立可进行判别和预测的现象模型。 空间回归分析，是将传统线性 回归分析结合空间权重矩阵，从而 引入空间依赖的一种分析方法。由 于一部分地理现象在空间上的分 布是连续的，传统的线性回归模型 针对离散的数据进行方程的拟合， 没有对数据之间的关联影响进行 处理，因而计算出的模型具有一定 的局限性。选择适当的空间关联模 型，有助于地理现象关联模式的探 索，结合空间依赖关系的回归模 型，本身保持了一定的连续性，因 此在描述地理现象上具有更大的 优势。 针对空间滞后性的特点，主要 有三种空间滞后模型，划分的标准 即在于空间滞后性的出现主要发 生在因变量还是自变量[8]。本文以 自变量的空间回归模型为例，假设 空间滞后出现在自变量 X上： Y=Xβ1+WXβ2+ε 式中，X为 k×n的自变量矩阵， W 为进行了行标准化的 n×n 空间 权重矩阵，β为 k×1的系数向量，ρ 是高阶模型的空间自回归系数，ε 是 n×1回归误差向量。 3 空间回归分析模型解算 下面以全国各省为单元，解算 各省一年内主要工业资源常量（煤 炭、电力、钢铁）与工业产值关系的 空间回归模型。 3.1 统一量纲 在进行回归分析计算之前，由 于各项指标的量纲不同，需要统一 进行数据变换，如中心化变换、极 差标准化变换、规格化变换等。规 格化变换的 公式 小学单位换算公式大全免费下载公式下载行测公式大全下载excel公式下载逻辑回归公式下载 如下： X ' ij= Xij - min{Xij} max{Xij}- min{Xij} Xij表示第 i 行、第 j 列数据， min{Xij}，max{Xij}分别表示第 j列最 小值和最大值。 def f（m）: return [[（float（m[i][j]）- min（zip （*m）[j]））/ （max （zip （*m）[j]）- min（zip （*m）[j]）） for j in range（len（m[i]））] for i in range（len（m））] 采用如下标准化变换公式，可 转化为均值为 0，方差为 1的数据： X ' ij =（Xij -X j）/Sj Xij表示第 i行、第 j列数据，X j， Sj分别表示第 j列数据的均值和方 差。 def f2（m）: return [[ （m[i][j]- sum（zip（*m）[j]）/float （len（zip（*m）[j]）））/ （sum（[（k- sum（zip（*m）[j]） /float（len（zip（*m）[j]）））**2 for k in zip （*m）[j]]）/len（zip （*m）[j]）） for j in range（len（m[i]））] for i in range（len（m））] 3.2 建立空间权重矩阵 空间权重矩阵表达了不同空 间单元之间的空间关联关系。本文 所用数据为国家基础地理信息数 67 2 0 0 9. 0 6 No. 3 WWW. GWN. ORG. CN地理信息世界 GEOMATICS WORLD 学术研究 据（可从国家基础地理信息网站下 载），通过 python 的 ogr 库结合 shapely库可以计算出各省之间的 邻接关系。 用 ogr读取 shp文件，必要时 需将同省的多边形进行合并，把结 果存储在一个以各省名为 key的 dictionary（名值对）数据结构中备 用。以下是读取 shp文件放入 dic- tionary中的代码。 from osgeo import ogr #从 bou2_4p.shp中读入 shp = ogr.Open（'bou2_4p.shp'） layer = shp.GetLayer（0） features = [layer.GetFeature（i） for i in range（layer.GetFeatureCount （））] provinces = {} for f in features: if f.GetFieldAsString（'name'） not in provinces.keys（）: provinces[f.GetFieldAsString （'name'）] = f else: geo = provinces[f.GetFieldAsS- tring（'name'）].GetGeometryRef（） geo2 = f.GetGeometryRef（） provinces[f.GetFieldAsString （'name'）].SetGeometry（geo.Union （geo2）） Python的 ogr库是 C++ GDAL /OGR库的一个包装，用于实现矢 量数据的存取和简单的几何操作。 然后即可对存储在 provinces 字典中的各省几何数据进行两两 的空间关系运算。这里采用一种名 为 DE- 9IM的空间位置关系判别方 法[9]。根据 DE- 9IM的命名规则，多 边形间的邻接关系以字符串 “FF2F11212”表示。通过 shapely中 几何对象的 relate方法可以比较两 图形的位置关系，将以上获得的各 省几何数据进行此操作即可获得 邻接关系矩阵[8]。 由 ogr 的 Geometry 对象到 shapely的 Geometry对象的数据交 换过程可以采用WKT（Well- Known Text）作为中间介质，如： >>> from shapely import wkt >>> heilongjiang = wkt.loads（str （f.GetGeometryRef（））） >>> print heilongjiang.relate （neimenggu） FF2F11212 这样，通过比较两两几何对象 的 relate方法返回值即可确定其位 置关系是否满足邻接，由此建立全 国各省的几何关系邻接矩阵。 在这个工业产值和资源的模 型中，仅仅依靠邻接关系建立的空 间权重矩阵并不一定可靠。由于现 代交通运输发达，即使不是邻接的 空间对象，依然可能产生相互影 响。因此可以尝试建立更加切合实 际情况的距离权重矩阵。考虑到计 算的简便，距离权重矩阵也采用了 二值化的方式，阈值为最大距离的 1/8。进而采用 ogr库读取基础地理 信息，进行距离的计算，建立矩阵。 3.3 全局空间自相关 Moran’s I 检验 应用全局空间自相关 Moran’s I检验函数如下： I= nS0 · n i Σ n j = 1 Σωij（Xi-X）（Xj-X） n i Σ（Xi-X） 2 若采用传统的编程方式利用 上述公式，需要进行循环，多次存 储状态变量，编程效率较低。若采 用 Python的函数式编程，则可直接 以类似公式的形式书写代码，不必 纠缠于状态变量的维护，部分代码 如下： I = [] for c in range（3）: avg = sum（[i[c] for i in source]） /n total = sum（[sum（map（ lambda x，y，z:x×y×z， weight[i]， [source[i][c]- avg]·n， [s[c]- avg for s in source]））for i in range（n）]） m = sum（map（lambda x: （x- avg）**2，[i[c] for i in source]）） I.append（n/s0*total/m） 在以上的函数式编程代码中， 诸如这样复杂的计算可以在一行 中计算完成。代码的表达和数学公 式的表达恰好是对应的， 下面给出焦炭产量、发电量和 钢铁产量三个因子在不同量纲统 一方式和不同空间权重矩阵组合 下的Moran’s I计算结果。 规格化量纲、邻接矩阵模型： [0.399 070，0.287 646，0.410 155] 规格化量纲、距离权重矩阵模型： [0.487 644，0.271 192，0.412 432] 标准化量纲、邻接矩阵模型： [0.399 309，0.287 250，0.409 710] 标准化量纲、距离权重矩阵模型： [0.488 278，0.270 333，0.412 223] 检测结果显示第四组模型的 空间自相关程度最高。因此采用第 四组模型进行后续的空间回归分 析。 3.4 空间滞后模型 本文选择标准化量纲距离权 重矩阵模型进行空间滞后模型拟 合，同时采用 Python函数式编程方 68 2 0 0 9. 0 6 No.3 WWW. GWN. ORG. CN 地理信息世界 GEOMATICS WORLD 学术研究 法，结合既有的矩阵运算类库，实 现系数向量的求解。 对于本模型三列的情况，可以 使用以下方式实现每一行的空间 权重赋值。 for i in range（n）: wx1.append（sum（map（lambda x， y: x·y，x1，weight[i]））） wx2.append（sum（map（lambda x， y: x·y，x2，weight[i]））） wx3.append（sum（map（lambdax，y :x·y，x3，weight[i]））） wy.append（sum（map（lambda x，y: x·y，y，weight[i]））） 对于最小二乘法的计算模型： β=（XTX）- 1XTY 那么对于自变量的空间滞后 模型，X，Y矩阵的形式为： X，Y=zip（*[x1，x2，x3，wx1，wx2， wx3，const]），zip（*[y]） 这样，结合一个矩阵运算的函 数库即可解算出 β矩阵，得自变量 的空间回归模型： y=X· - 0.202 6 0.830 3 0.164 ! " " " " " " " " " # $ % % % % % % % % % &7 +WX· - 0.522 3 0.080 1 0.449 ! " " " " " " " " " # $ % % % % % % % % % &2 +0.009 4 4 运算脚本性能分析 对于函数式编程的性能，可以 通过 Python 自带的 Profile 进行运 行时间的统计。以计算Moran’s I的 脚本为例，我们可以比较函数式编 程和采用传统方法编程的效率。 在执行 Python 程序时，通过 python m cProfile ***.py 的方式运 行，调用 cProfile 模块进行 Profile 统计。 对于函数式编程，Profile 结果 如下： 3440 function calls（3409 prim- itive calls）in 0.019 CPU seconds 对于传统的编程方式，采用循 环遍历、状态变量存储，Profile结果 如下： 461 function calls（430 primitive calls）in 0.009 CPU seconds 函数式编程执行了更多的函 数调用，在效率上有一定损失。 5 结 语 通过以上示例，Python 函数式 编程的优缺点已经非常清晰。在较 为复杂的空间分析运算中，面对大 量的公式和数据，采用传统的语言 编程过程过于复杂，而纯粹通过数 学软件来实现需要完成一定工作 量的数据准备，数学功能相对独 立，不易与既有系统衔接。 Python作为 ArcGIS 9.0的使用 脚本语言，在 GIS上已经得到了一 定程度的推广，而其开源 GIS类库 直接继承自比较成熟 C++ 系， GDAL/OGR，Shapely，RTree 等都处 在一个良好的发展态势。此外，除 GIS相关功能以外，Python 完全可 以胜任从 Web 开发到跨平台桌面 开发、游戏开发等各个方面。随着 脚本语言 Ruby在Web开发上暴发 式的巨大成功和 Linux操作系统的 进一步推广，包括 Python在内的其 他一些语言也得到了更大的关注， 社区发展迅速，未来必然将会出现 更多的 Python项目。对于长期遭受 诟病的解释型脚本语言的效率，尽 管 Python的运行效率低于 C语言， 但是随着新版本的开发，Python 的 性能不断提升。另一方面，由于可 以和编译后的 C模块共同使用，可 将性能要求较高的部分以 C语言 实现。 Python的最大优势在于极高的 编程效率，如 SDA这样需要不断进 行探索计算，进行模型的变换、参 数的修改而获得最优解的研究方 法，对于编程工具的敏捷性要求则 更高，而这正是脚本语言的优势所 在。作为一门结合了函数式编程的 脚本语言，Python的优势更加明显。 因此，Python 语言在 GIS研究中的 地位也必将随着人们对它的了解 加深而显得越发重要。 参考文献 [1] 周兰霞. ArcGIS9.0在苯丙酮尿症空间 分析中的应用 [J]. 中国卫生统计， 2007，24（6）: 604. [2] 刘承良.武汉都市圈经济社会要素流 的空间分析 [J]. 人文地理，2007，22 （6）: 30. [3] Python Official Documentations[EB/OL]. 2008- 04- 16. http:// www.python.org/a- bout/ . [4] Wikipedia: Functional Programming. [EB/OL].2008- 04- 16.http://en.wikiped- ia.org/wiki/Functional_programming [5] John Hughes. Why Functional Program- ming Matters [EB/OL].2008- 04- 18.http: // www.cschalmers.se/~rjmb/papers/vhy- fp.html. [6] Mark Lutz . Learning Python 3rd Edition [M]. O’Reilly，2007: 350- 353. [7] 朱欣娟.基于 GIS的空间分析及其发 展研究 [J]. 计算机 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 与应用 2002，18（1）: 1. [8] L.Anselin. SPATIAL ECONOMETRICS [M]. Springer，1988: 5- 11. [9] OpenGIS Implementation Specification for Geographic information - Simple fea- ture access - Part 1: Common architec- ture [EB/OL]. 2008- 03- 20.http://www. epengeospatial.org/standasds/sfa. 69