基于python的栅格数据标准差椭圆

基于python的栅格数据标准差椭圆

方法介绍

标准差椭圆空间统计方法（standard deviational elipse,SDE）能够精确地揭示地理要素空间分布整体特征。该方法基于研究对象的空间区位和空间结构，从全局性空间的角度定量解释地理要素空间分布的中心性、展布性、方向性、空间形态等特征，揭示地理要素的空间分布及时空演化过程。其计算过程如下：

1. 计算平均/加权平均中心：
   • 对于点集 $(x_i,y_i)$ 和对应权重 $w_i$，加权平均中心计算为：$$\begin{array}{r}\bar{x}=\frac{{\sum }_{i=1}^n{w}_i{x}_i}{{\sum }_{i=1}^n{w}_i}\\ \bar{y}=\frac{{\sum }_{i=1}^n{w}_i{y}_i}{{\sum }_{i=1}^n{w}_i}\end{array}$$
2. **计算加权协方差矩阵：$$\begin{array}{rl}\Sigma & =\left[\begin{array}{cc}Cov(X,X)& Cov(X,Y)\\ Cov(Y,X)& Cov(Y,Y)\end{array}\right]\\ Cov(X,Y)& =\sum _{i=1}^n{w}_i\cdot (x_i-\bar{x})\cdot (y_i-\bar{y})\end{array}$$
3. 特征值和特征向量：
   • 求解特征方程：$$|\Sigma -\lambda I|=0$$
     得到特征值 ${\lambda }_1,{\lambda }_2$ 和对应的特征向量 $\overrightarrow{v_1},\overrightarrow{v_2}$
4. 椭圆参数计算:
   • 长短半轴：$$\begin{array}{r}\alpha =\sqrt{{\lambda }_1{\chi }_{2,p}^2}\\ \beta =\sqrt{{\lambda }_2{\chi }_{2,p}^2}\end{array}$$ 其中 ${\chi }_{2,p}^2$ 是自由度为 2，置信水平为 p 的卡方分布临界值（两倍标准差即为 0.63）
   • 旋转角度：$$\theta =\arctan (\frac{v_{12}}{v_{11}})$$ 其中 $v_{11},v_{12}$ 是第一特征向量的分量
5. 计算椭圆参数方程：$$\begin{array}{rl}x(t)& =\bar{x}+\alpha \cos (t)\cos (\theta )-\beta \sin (t)\sin (\theta )\\ y(t)& =\bar{y}+\alpha \cos (t)\sin (\theta )+\beta \sin (t)\cos (\theta )\end{array}$$ 其中 $t\in [0,2\pi ]$

这些公式主要由以下函数来计算：

• np.average() 计算加权平均
• np.cov() 计算加权协方差矩阵
• np.linalg.eig() 计算特征值和特征向量
• chi2.ppf() 计算卡方分布临界值

代码展示

• 首先需要读取栅格数据为数组形式进行计算，这里使用 rasterio 库：

def raster_to_points(raster_path):# 读取栅格数据with rasterio.open(raster_path) as src:# 读取栅格数据的第一波段band1 = src.read(1)nodata = src.nodata# 获取栅格数据的行列号rows, cols = np.where(band1 != nodata)# 将行列号转换为地理坐标x_coords, y_coords = src.xy(rows, cols)# 将坐标组合成点数据points = np.column_stack((x_coords, y_coords))# 假设权重为栅格值weights = band1[rows, cols]return points, weights, src.crs

• 然后计算椭圆的参数，这里使用 numpy 和 scipy.stats 库：

def calculate_ellipse(data, weight, confidence_level):# 计算均值mean = np.average(data, axis=0, weights=weight)# 计算加权协方差矩阵cov_matrix = np.cov(data, rowvar=False, aweights=weight)# 协方差矩阵的特征值和特征向量eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)# 计算椭圆长短半轴chi2_val = chi2.ppf(confidence_level, df=2)alpha = np.sqrt(eigenvalues[0] * chi2_val)beta = np.sqrt(eigenvalues[1] * chi2_val)# 计算椭圆角度angle = np.degrees(np.arctan2(*eigenvectors[:, 0][::-1]))# 绘制椭圆ellipse = Ellipse(xy=mean, width=alpha * 2, height=beta * 2, angle=angle, fc='None', lw=2)ellipse_info = {'mean': mean,'alpha': alpha,'beta': beta,'angle': angle}return ellipse, ellipse_info

• 最后生成椭圆方程，并用 geopandas 库给椭圆添加地理信息，保存为矢量文件：

def create_ellipse_polygon(mean, alpha, beta, angle, num_points=100):# 生成椭圆上的点theta = np.linspace(0, 2 * np.pi, num_points)# 椭圆参数方程x = mean[0] + alpha * np.cos(theta) * np.cos(np.radians(angle)) - beta * np.sin(theta) * np.sin(np.radians(angle))y = mean[1] + alpha * np.cos(theta) * np.sin(np.radians(angle)) + beta * np.sin(theta) * np.cos(np.radians(angle))# 创建多边形return Polygon(zip(x, y))
def process_raster(raster_path, output_shp_path):...gdf = gpd.GeoDataFrame({'confidence': ['63%'],  # 置信度'center_x': [mean[0]],  # 中心点 X 坐标'center_y': [mean[1]],  # 中心点 Y 坐标'major_axis': [alpha],  # 长半轴'minor_axis': [beta],   # 短半轴'angle': [angle],       # 旋转角度'geometry': [ellipse_polygon]}, crs=crs)gdf.to_file(output_shp_path, driver='ESRI Shapefile')...

完整代码

import numpy as np
from matplotlib.patches import Ellipse
from scipy.stats import chi2
import rasterio
import geopandas as gpd
from shapely.geometry import Polygon
import os
def calculate_ellipse(data, weight, confidence_level):# 计算均值mean = np.average(data, axis=0, weights=weight)# 计算加权协方差矩阵cov_matrix = np.cov(data, rowvar=False, aweights=weight)# 协方差矩阵的特征值和特征向量eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)# 计算椭圆长短半轴chi2_val = chi2.ppf(confidence_level, df=2)alpha = np.sqrt(eigenvalues[0] * chi2_val)beta = np.sqrt(eigenvalues[1] * chi2_val)# 计算椭圆角度angle = np.degrees(np.arctan2(*eigenvectors[:, 0][::-1]))# 绘制椭圆ellipse = Ellipse(xy=mean, width=alpha * 2, height=beta * 2, angle=angle, fc='None', lw=2)ellipse_info = {'mean': mean,'alpha': alpha,'beta': beta,'angle': angle}return ellipse, ellipse_info
def raster_to_points(raster_path):# 读取栅格数据with rasterio.open(raster_path) as src:# 读取栅格数据的第一波段band1 = src.read(1)nodata = src.nodata# 获取栅格数据的行列号rows, cols = np.where(band1 != nodata)# 将行列号转换为地理坐标x_coords, y_coords = src.xy(rows, cols)# 将坐标组合成点数据points = np.column_stack((x_coords, y_coords))# 假设权重为栅格值weights = band1[rows, cols]return points, weights, src.crs
def create_ellipse_polygon(mean, alpha, beta, angle, num_points=100):# 生成椭圆上的点theta = np.linspace(0, 2 * np.pi, num_points)# 椭圆参数方程x = mean[0] + alpha * np.cos(theta) * np.cos(np.radians(angle)) - beta * np.sin(theta) * np.sin(np.radians(angle))y = mean[1] + alpha * np.cos(theta) * np.sin(np.radians(angle)) + beta * np.sin(theta) * np.cos(np.radians(angle))# 创建多边形return Polygon(zip(x, y))
def process_raster(raster_path, output_shp_path):# 创建输出目录（如果不存在）os.makedirs(os.path.dirname(output_shp_path), exist_ok=True)# 处理栅格数据points, weights, crs = raster_to_points(raster_path)weights = (weights - np.min(weights))# 计算椭圆_, ellipse_info = calculate_ellipse(points, weights, 0.63)mean = ellipse_info['mean']alpha = ellipse_info['alpha']beta = ellipse_info['beta']angle = ellipse_info['angle']# 保存为shpellipse_polygon = create_ellipse_polygon(mean, alpha, beta, angle, num_points=100)gdf = gpd.GeoDataFrame({'confidence': ['63%'],  # 置信度'center_x': [mean[0]],  # 中心点 X 坐标'center_y': [mean[1]],  # 中心点 Y 坐标'major_axis': [alpha],  # 长半轴'minor_axis': [beta],   # 短半轴'angle': [angle],       # 旋转角度'geometry': [ellipse_polygon]}, crs=crs)gdf.to_file(output_shp_path, driver='ESRI Shapefile')
if __name__ == '__main__':# 输入输出路径raster_path = 'input.tif'  # 输入栅格文件路径output_shp_path = 'output.shp'  # 输出shapefile路径# 处理数据process_raster(raster_path, output_shp_path)

结论

以大同市 2019 年的夜间灯光指数为例，程序运行结果与 Arcgis 地理处理中的标准差椭圆基本一致：（由于计算的浮点数问题存在较小偏差）

参考资料：
方向分布（标准差椭圆）算法
标准差椭圆算法实现