基于python的空间距离权重计算——arcgis中的arcpy

2019独角兽企业重金招聘Python工程师标准>>>

arcgis出了10后，python代替vb成为官方脚本语言，python中的gis库也是有很多，但因为本人最先接触的是arcgis所以就直接学习他自带的库arcpy了，主要以空间权重矩阵作为一个契机顺便学习arcpy。

代码及数据资料：https://github.com/shikanon/WeightDistance

因为还是不太习惯arcpy中的表操作，所以借助了pandas进行表操作，构建一个表转换：

#ListFields包含field类的数组
fields=arcpy.ListFields(path)

def GetTable():
    '''将arcpy表单变为pandas表单，还是喜欢pandas些~'''
    table=[]
    fieldname=[field.name for field in fields]
    #游标集合，用for 循环一次后没办法循环第二次!一个游标实例只能循环一次
    data=arcpy.SearchCursor(path)
    for row in data:
        #Shape字段中的要数是一个几何类
        r=[]
        for field in fields:
            r.append(row.getValue(field.name))
        table.append(r)
    return pd.DataFrame(table,columns=fieldname)

因为求欧几里距离在arcpy 中polyline直接有length属性所以只需要构建一个线类，然后直接求长度就行了

def getdistance(point1,point2,sparef):
    '''求两点的距离'''
    #两点构造一条线
    l=arcpy.Polyline(arcpy.Array([point1,point2]),sparef)
    #求线长
    return l.getLength()

在求长度时候需要考虑到shp的投影坐标，因此需要先得到一个spatialReference

#空间坐标投影参数
#随便抽取一个shape出来求得该图层的空间投影
sparef=data['Shape'][0].spatialReference
print sparef.name
#计算质心centre
#Shape字段下是一个geometry object对象,其下有中心点和面积等工具
#貌似由centroid抽出来的点没有空间投影属性,
#所以data['centre'][0].spatialReference将会报错
data['centre']=[d.centroid for d in data['Shape']]

#计算欧几里权重矩阵
weight=[]
for i in data['centre']:
    row=[]
    for j in data['centre']:
        row.append(getdistance(i,j,sparef))
    weight.append(row)
weight=pd.DataFrame(weight,index=data['FID'],columns=data['FID'])

最后导出csv:

#输出csv
weight.to_csv(path_or_buf='weight.csv',sep=',')

小小练习完成，查看一下结果，图为广东省市级单元的空间权重矩阵：

贴上完整代码：

# -*- coding: cp936 -*-

import arcpy
import numpy as np
import pandas as pd

path='data/gd.shp'
#ListFields包含field类的数组
fields=arcpy.ListFields(path)

def GetTable():
    '''将arcpy表单变为pandas表单，还是喜欢pandas些~'''
    table=[]
    fieldname=[field.name for field in fields]
    #游标集合，用for 循环一次后没办法循环第二次!一个游标实例只能循环一次
    data=arcpy.SearchCursor(path)
    for row in data:
        #Shape字段中的要数是一个几何类
        r=[]
        for field in fields:
            r.append(row.getValue(field.name))
        table.append(r)
    return pd.DataFrame(table,columns=fieldname)

def getdistance(point1,point2,sparef):
    '''求两点的距离'''
    #两点构造一条线
    l=arcpy.Polyline(arcpy.Array([point1,point2]),sparef)
    #求线长
    return l.getLength()

data=GetTable()
#空间坐标投影参数
#随便抽取一个shape出来求得该图层的空间投影
sparef=data['Shape'][0].spatialReference
print sparef.name
#计算质心centre
#Shape字段下是一个geometry object对象,其下有中心点和面积等工具
#貌似由centroid抽出来的点没有空间投影属性,
#所以data['centre'][0].spatialReference将会报错
data['centre']=[d.centroid for d in data['Shape']]

#计算欧几里权重矩阵
weight=[]
for i in data['centre']:
    row=[]
    for j in data['centre']:
        row.append(getdistance(i,j,sparef))
    weight.append(row)
weight=pd.DataFrame(weight,index=data['FID'],columns=data['FID'])
#输出csv
weight.to_csv(path_or_buf='weight.csv',sep=',')

http://my.oschina.net/Kanonpy/

------------窈----------------窈---------------窈------------窈------------窈-------------窈------------窈-------

转载于:https://my.oschina.net/Kanonpy/blog/361078