python空间统计分析——汇总

空间点模式描述——空间分布特征识别功能实现，包括

平均中心
	
中位数中心
	
中心要素
	
标准差圆/标准距离
	
标准差椭圆/分布方向
python代码

# -*- coding: utf-8 -*-
# @Author : Zhao HL
# @contact: huilin16@qq.com
# @Time   : 2020/10/28 20:06
# @File   : spatial_point_analysis.py
from __future__ import print_function
from __future__ import division

import numpy as np
import geopandas as gpd
from shapely.geometry import Point
from shapely.affinity import scale, rotate
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist

shp_point_path = r'E:\Course\spatial_statistic\exper23\source_data\eqpoint.shp'
shp_polygon_path = r'E:\Course\spatial_statistic\exper23\source_data\world_prj.shp'


def get_shpDF(shp_path):
    '''
    read shp file, and return shp_df
    读取shp文件，返回shp的dataframe
    :param shp_path: shp file
    :return:
    '''
    shp_df = gpd.read_file(shp_path, encoding='gb18030')
    return shp_df


def shpDF_2_pos_forPoint(shp_df):
    '''
    give a Point shp_df, return x list, y list
    读取点文件的shp的dataframe，返回所有x的列表、所有y的列表
    :param shp_df:
    :return:
    '''
    if 'geometry' not in shp_df.columns:
        return None
    else:
        geometry_list = shp_df['geometry']
        xs = []
        ys = []
        for geo_pos in geometry_list:
            xs.append(geo_pos.x)
            ys.append(geo_pos.y)
        return [xs, ys]


def get_posMeanCenter(pos):
    '''
    平均中心，x、y分别为所有点的平均值
    :param pos:
    :return:
    '''
    xs, ys = pos
    pos_x = np.mean(xs)
    pos_y = np.mean(ys)
    pos_point = Point(pos_x, pos_y)
    return pos_point


def get_posMedianCenter(pos):
    '''
    中位数中心，x、y为到所有要素距离累加和最小的位置，可以由公式逐步迭代得到
    https://glenbambrick.com/tag/median-center/
    https://blog.csdn.net/allenlu2008/article/details/47752943

    The is no single formula or equation for calculating an exact Median Center, according to Andy Mitchell
    it is an iterative process involving calculating the Mean Center, summing the distances from it to each feature,
     offsetting the center slightly and summing the distances again until
     it eventually zones in on the optimum location that has the lowest sum.

    :param pos:
    :return:
    '''

    def geometric_median(X, eps=1e-5):
        y = np.mean(X, 0)

        while True:
            D = cdist(X, [y])
            nonzeros = (D != 0)[:, 0]
            Dinv = 1 / D[nonzeros]
            Dinvs = np.sum(Dinv)
            W = Dinv / Dinvs
            T = np.sum(W * X[nonzeros], 0)
            num_zeros = len(X) - np.sum(nonzeros)
            if num_zeros == 0:
                y1 = T
            elif num_zeros == len(X):
                return y
            else:
                R = (T - y) * Dinvs
                r = np.linalg.norm(R)
                rinv = 0 if r == 0 else num_zeros / r
                y1 = max(0, 1 - rinv) * T + min(1, rinv) * y
            if np.linalg.norm(y - y1) < eps:
                return y1
            y = y1

    xy_arr = np.array(pos).transpose(1, 0)
    test_median = geometric_median(xy_arr)
    pos_point = Point(test_median[0], test_median[1])
    return pos_point


def get_posCentralFeature(pos):
    '''
    中心要素，x、y为到其他要素距离累加和最小的要素
    :param pos:
    :return:
    '''
    xs, ys = pos
    xy_arr = np.array(pos).transpose(1, 0)
    dsts = np.zeros(len(xs))
    for index, (pos_x, pos_y) in enumerate(zip(xs, ys)):
        dst_all = cdist(xy_arr, [[pos_x, pos_y]])
        dst_sum = np.sum(dst_all)
        dsts[index] += dst_sum
    dst_index = np.argmin(dsts)
    pos_point = Point(pos[0][dst_index], pos[1][dst_index])
    return pos_point


def get_circleDistandard(pos):
    '''
    标准差圆/标准距离
    https://desktop.arcgis.com/zh-cn/arcmap/10.3/tools/spatial-statistics-toolbox/h-how-standard-distance-spatial-statistic-works.htm
    :param pos:
    :return:
    '''
    meanCenter = get_posMeanCenter(pos)
    xy_arr = np.array(pos).transpose(1, 0)
    dst_all = cdist(xy_arr, [[meanCenter.x, meanCenter.y]])
    radius = np.sqrt(np.sum(dst_all ** 2) / (len(xy_arr) - 2))
    circle = meanCenter.buffer(radius)
    return circle


def get_ellipseDistandard(pos):
    '''
    标准差椭圆/方向分布
    https://desktop.arcgis.com/zh-cn/arcmap/10.3/tools/spatial-statistics-toolbox/h-how-directional-distribution-standard-deviationa.htm
    https://blog.csdn.net/allenlu2008/article/details/47780405
    https://gis.stackexchange.com/questions/243459/drawing-ellipse-with-shapely
    :param pos:
    :return:
    '''
    meanCenter = get_posMeanCenter(pos)
    xs, ys = pos
    xs = np.expand_dims(np.array(xs), -1)
    ys = np.expand_dims(np.array(ys), -1)

    delta_x = xs - meanCenter.x
    delta_y = ys - meanCenter.y
    delta_x_square = np.sum(delta_x ** 2)
    delta_y_square = np.sum(delta_y ** 2)
    delta_xy = np.sum(delta_x * delta_y)

    tan_theta_a = delta_x_square - delta_y_square
    tan_theta_b = np.sqrt(np.square(delta_x_square - delta_y_square) + 4 * np.square(delta_xy))
    tan_theta_c = 2 * delta_xy
    tan_theta = (tan_theta_a + tan_theta_b) / tan_theta_c
    angle_theta = np.arctan(tan_theta)
    long_axis = np.sqrt(2 / len(xs) * np.sum((delta_x * np.cos(angle_theta) - delta_y * np.sin(angle_theta)) ** 2))
    short_axis = np.sqrt(2 / len(ys) * np.sum((delta_x * np.sin(angle_theta) + delta_y * np.cos(angle_theta)) ** 2))

    ellipse = Point((meanCenter.x, meanCenter.y)).buffer(1)
    ellipse = scale(ellipse, xfact=long_axis, yfact=short_axis)
    ellipse = rotate(ellipse, angle=-angle_theta, use_radians=True)
    # print(meanCenter.x, meanCenter.y, long_axis, short_axis, angle_theta, angle_theta / np.pi * 180)
    return ellipse


def shp_show(features, colors=None, cmaps=None, markers=None):
    fig, ax = plt.subplots()
    for index, feature in enumerate(features):
        if not isinstance(feature, gpd.GeoDataFrame) and not isinstance(feature, gpd.GeoSeries):
            feature = gpd.GeoSeries(feature)
        if colors != None and markers != None:
            feature.plot(ax=ax, color=colors[index], marker=markers[index], cmap=cmaps[index])
        else:
            feature.plot(ax=ax)
    plt.show()


if __name__ == '__main__':
    pass
    shp_df_polygon = get_shpDF(shp_polygon_path)
    shp_df_point = get_shpDF(shp_point_path)
    pos = shpDF_2_pos_forPoint(shp_df_point)

    meanCenter = get_posMeanCenter(pos)
    medianCenter = get_posMedianCenter(pos)
    centralFeature = get_posCentralFeature(pos)
    circleDistandard = get_circleDistandard(pos)
    ellipseDistandard = get_ellipseDistandard(pos)
    # shp_show([shp_df_polygon,shp_df_point,meanCenter],colors=['lightgreen','b','r'],cmaps=['YlGn',None,None],
    #          markers=[None,'.','*'])
    # shp_show([shp_df_polygon, shp_df_point, medianCenter], colors=['lightgreen', 'b', 'r'], cmaps=['YlGn', None, None],
    #          markers=[None, '.', '*'])
    # shp_show([shp_df_polygon, shp_df_point, centralFeature], colors=['lightgreen', 'b', 'r'], cmaps=['YlGn', None, None],
    #          markers=[None, '.', '*'])
    # shp_show([shp_df_polygon, circleDistandard, shp_df_point], colors=['lightgreen', 'r', 'b'], cmaps=['YlGn', None, None],
    #          markers=[None, None, '.'])
    # shp_show([shp_df_polygon, ellipseDistandard, shp_df_point], colors=['lightgreen', 'r', 'b'], cmaps=['YlGn', None, None],
    #          markers=[None, None, '.'])
    shp_show(
        [shp_df_polygon, circleDistandard, ellipseDistandard, shp_df_point, meanCenter, medianCenter, centralFeature],
        colors=['lightgreen', 'deepskyblue', 'orange', 'pink', 'r', 'g', 'b'], cmaps=['YlGn', None, None, None, None, None, None],
        markers=[None, None, None, '.', '*', '*', '*'])


 

Python空间数据处理环境搭建

作者：阿振

邮箱：tanzhenyugis@163.com

博客：https://blog.csdn.net/theonegis/article/details/80089375

修改时间：2018-04-26

优酷视频地址：http://v.youku.com/v_show/id_XMzU2NjAzMTk0NA==.html



Conda的下载和安装

什么是Conda? 官方定义：Package, dependency and environment management for any language—Python, R, Ruby, Lua, Scala, Java, JavaScript, C/ C++, FORTRAN

Conda就是一个虚拟环境和包（库）依赖管理工具

下载地址：Downloading conda

对于Windows版本的，确定Python版本和系统类型直接下载安装包进行安装

对于Linux和macOS系统，在Terminal中运行bash脚本进行安装即可。



Conda的使用



新建虚拟环境（Managing environments）

conda create -n osgeo



切换到新建的虚拟环境

source activate osgeo  （Linux和macOS）

activate osgeo  （Windows）



退出虚拟环境

source deactivate  （Linux和macOS）

deactivate  （Windows）

实用命令：


查看虚拟环境列表  conda env list 或者 conda info --envs
删除虚拟环境  conda remove --name <environment> --all
查看虚拟环境中的包列表  conda list
更新conda或者某个包  conda update conda 或者 conda update <package>
更新虚拟环境下的所有包  conda update --all
查看过时的包  conda search --outdated
搜索指定包  conda search <package>
删除某个包 conda remove <package>
添加channel到conda配置文件  conda config --add channels <channel> 或者 conda config --append channels <channel>




空间数据处理Python库的安装



常用的空间数据处理Python库


GDAL 全能型的基础空间数据处理库
fiona 基于GDAL的空间矢量数据处理库
rasterio 基于GDAL的空间栅格处理库
basemap 基于matplotlib的空间制图库
GeoPandas 基于pandas的空间数据分析库
RSGISlib 针对遥感数据及GIS分析的高级库




使用conda进行库的安装

打开命令行工具（Terminal），输入命令，进入虚拟环境


安装GDAL库

conda install -c conda-forge gdal
安装fiona库

conda install -c conda-forge fiona
安装rasterio库

conda install -c conda-forge rasterio

​




使用pip进行库的安装

什么是pip呢？pip是Python默认和推荐实用的包管理工具，可以利用pip从PyPI网络仓库自动下载Python包进行安装和管理。

对于Windows下的二进制库的预编译包，提供给大家一个网站：Unofficial Windows Binaries for Python Extension Packages

使用pip安装的时候，如果该Python包底层依赖一些C++库，则需要手动进行编译，或者安装指定平台下预编译好的库。


安装GDAL库

pip install GDAL‑2.2.4‑cp37‑cp37m‑win_amd64.whl
安装fiona库

pip install Fiona‑1.7.11.post1‑cp37‑cp37m‑win_amd64.whl
安装rasterio库

pip install rasterio‑1.0a12‑cp37‑cp37m‑win_amd64.whl
安装Jupyter

pip install jupyter

使用Jupyter Notebook进行编程

jupyter notebook

学python空间数据处理不一定要选课，哈哈哈
最近学校开了门python空间数据处理，感觉有不止1mol（1mol=6.02*10^23）人去上这门课，说句实在话，如果你知道椭球体投影和坐标那些东西，并系统地学过一门面向对象语言（不论它是C++,C#或者是java），并掌握了一些基础的数据结构知识，这门课你不用去上了。我个人猜测大家和我一样是更想要学分，并不是想学东西，因为你想学的东西的话，资源就多了而且并不一定要上课。不过我听过直播，宋老师真是个贴心的好人，连命令行窗口都手把手教。

闲话少说，为了给空间数据处理做准备，我准备了一些用的比较多的库，如果后续有补充，我将进行博客连载，当然方法最重要了，一通百通才是王道。

1. 主流的python库安装方法不一定成功
这里主流的python库安装方法指的是：在安装pip的前提下，直接使用pip install 库名称`就好了，当然了，系统默认的镜像源慢得像蜗牛一样，但是清华的镜像源又下架了，我推荐使用豆瓣的镜像源，网址是：https://pypi.doubanio.com/simple，非常简单，你只需要pip install 库名称==版本名 -i 镜像源地址，大致就这么个简单的语法规则，以numpy库安装为例，在命令行窗口输入下面的语句即可完成numpy库的安装

pip install numpy -i https://pypi.doubanio.com/simple


这里我要说明一下，如scipy，pandas，matplotlib都可以用这种方法安装成功，其实我也是比较推荐使用这种方法的，实在不行再用下面的安装方法。

2.whl文件的安装方法
目前，就小编我的孤陋寡闻，whl文件安装只要你下载文件版本对应，成功率100%吧，哈哈哈。我在这里推荐一个很牛逼的whl文件下载网址：

https://www.lfd.uci.edu/~gohlke/pythonlibs/

这个网址基本全了你要的库，安装方法很简单，复杂的话我就不敢拿出来了，找到你要的库，下载到本地电脑文件夹中，这里以GDAL库安装为例：



在网页中找到你要的whl文件，这里cp37代表支持python3.7，cp36代表支持python3.6，对应下载即可。

这里以GDAL库安装为例（因为主流方法应该是不能成功了），因为我是python3.7版本，我下载了GDAL-2.4.1-cp37-cp37m-win_amd64.whl这个文件，然后呢，打开你的命令行窗口，找到该文件所在文件夹（一般我是通过右键在文件属性里复制粘贴来找的），然后输入

cd C:\Users\Sun Strong\Desktop

，C:\Users\Sun Strong\Desktop是该文件所在的文件夹，然后输入

pip install GDAL-2.4.1-cp37-cp37m-win_amd64.whl

，其中GDAL-2.4.1-cp37-cp37m-win_amd64.whl是文件全名（文件名+后缀）。



由于我的库里面已经有了GDAL，所示pip提示我已经安装了。

pyproj和scikit_learn都可以使用该方法安装成功。

3.ospybook是特殊的存在
值得一提的是，ospybook你可以用下列命令完成安装：

pip install https://github.com/cgarrard/osgeopy-code/raw/master/ospybook-latest.zip

到这儿，仅仅一小部分的python库算是安装好了，下次博哥我给大家带来我认为最好用的python集成开发环境。

Conda的下载和安装
什么是Conda? 官方定义：Package, dependency and environment management for any language—Python, R, Ruby, Lua, Scala, Java, JavaScript, C/ C++, FORTRAN
Conda就是一个虚拟环境和包（库）依赖管理工具
对于Windows版本的，确定Python版本和系统类型直接下载安装包进行安装
对于Linux和macOS系统，在Terminal中运行bash脚本进行安装即可。
Conda的使用
conda create -n osgeo
切换到新建的虚拟环境
source activate osgeo  （Linux和macOS）
activate osgeo  （Windows）
退出虚拟环境
source deactivate  （Linux和macOS）
deactivate  （Windows）
实用命令：查看虚拟环境列表  conda env list 或者 conda info --envs
删除虚拟环境  conda remove --name  --all
查看虚拟环境中的包列表  conda list
更新conda或者某个包  conda update conda 或者 conda update 
更新虚拟环境下的所有包  conda update --all
查看过时的包  conda search --outdated
搜索指定包  conda search 
删除某个包 conda remove 
添加channel到conda配置文件  conda config --add channels  或者 conda config --append channels 
空间数据处理Python库的安装
常用的空间数据处理Python库GDAL 全能型的基础空间数据处理库
fiona 基于GDAL的空间矢量数据处理库
rasterio 基于GDAL的空间栅格处理库
basemap 基于matplotlib的空间制图库
GeoPandas 基于pandas的空间数据分析库
RSGISlib 针对遥感数据及GIS分析的高级库
使用conda进行库的安装
打开命令行工具（Terminal），输入命令，进入虚拟环境安装GDAL库
conda install -c conda-forge gdal
安装fiona库
conda install -c conda-forge fiona
安装rasterio库
conda install -c conda-forge rasterio
​
使用pip进行库的安装
什么是pip呢？pip是Python默认和推荐实用的包管理工具，可以利用pip从PyPI网络仓库自动下载Python包进行安装和管理。
使用pip安装的时候，如果该Python包底层依赖一些C++库，则需要手动进行编译，或者安装指定平台下预编译好的库。安装GDAL库
pip install GDAL‑2.2.4‑cp37‑cp37m‑win_amd64.whl
安装fiona库
pip install Fiona‑1.7.11.post1‑cp37‑cp37m‑win_amd64.whl
安装rasterio库
pip install rasterio‑1.0a12‑cp37‑cp37m‑win_amd64.whl
安装Jupyter
pip install jupyter
使用Jupyter Notebook进行编程
jupyter notebook