没有网站做APP,尚层装饰,官方网站建设方案,惠州做网站建设价格Spatial Data Analysis#xff08;三#xff09;#xff1a;点模式分析
---- 1853年伦敦霍乱爆发
在此示例中#xff0c;我将演示如何使用 John Snow 博士的经典霍乱地图在 Python 中执行 KDE 分析和距离函数。
感谢 Robin Wilson 将所有数据数字化并将其转换为友好的 G…Spatial Data Analysis三点模式分析
---- 1853年伦敦霍乱爆发
在此示例中我将演示如何使用 John Snow 博士的经典霍乱地图在 Python 中执行 KDE 分析和距离函数。
感谢 Robin Wilson 将所有数据数字化并将其转换为友好的 GIS 格式。 原始数据从这里获得 http://blog.rtwilson.com/john-snows-known-cholera-analysis-data-in-modern-gis-formats/
具体来说我们需要这些包来进行分析
rasterio这是一个用于在 python 中处理栅格数据的包。 我们不使用栅格但在这里使用它的原因是显示 GeoTiff 图像作为底图该图像经过校正以与 GIS shapefile 对齐。seaborn这实际上是一个可视化包 但是它们确实具有可用于二维空间数据的 KDE 功能。pointpats这是一个用于进行 PPA 和最近邻分析的软件包。 我们需要用它来计算例如 G距离函数。
步骤1
导入所有需要的包
import geopandas as gpd
import rasterio
from rasterio.plot import show
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
from pointpats import distance_statistics as stats
from pointpats import PointPattern, PoissonPointProcess第2步
读入两个 shapefile
SnowGIS/Cholera_Deaths.shp 是一个包含所有死亡和位置的点形状文件SnowGIS/Pumps.shp 是包含所有泵的点形状文件
由于它们都是 shapefile一旦我们使用“gpd.read_file()”将它们读入 python我们就得到了 GeoDataFrame。
cases gpd.read_file(https://raw.githubusercontent.com/Ziqi-Li/GEO4162C/main/data/SnowGIS/Cholera_Deaths.shp)pump gpd.read_file(https://raw.githubusercontent.com/Ziqi-Li/GEO4162C/main/data/SnowGIS/Pumps.shp)在我们进行下一步操作之前最好先检查数据。 例如您需要执行以下操作
1.查看数据表了解GeoDataFrame中的属性2.检查几何图形看看是否可以绘制它们3.检查您读入的shapefile是否具有相同的crs
cases.head() #This returns you the first 5 rows in the GeoDataFrameIdCountgeometry003POINT (529308.741 181031.352)102POINT (529312.164 181025.172)201POINT (529314.382 181020.294)301POINT (529317.380 181014.259)404POINT (529320.675 181007.872) scriptconst buttonEl document.querySelector(#df-149da3d4-d8f5-4e15-82d9-0a5e55507ee2 button.colab-df-convert);buttonEl.style.display google.colab.kernel.accessAllowed ? block : none;async function convertToInteractive(key) {const element document.querySelector(#df-149da3d4-d8f5-4e15-82d9-0a5e55507ee2);const dataTable await google.colab.kernel.invokeFunction(convertToInteractive,[key], {});if (!dataTable) return;const docLinkHtml Like what you see? Visit the a target_blank hrefhttps://colab.research.google.com/notebooks/data_table.ipynbdata table notebook/a to learn more about interactive tables.;element.innerHTML ;dataTable[output_type] display_data;await google.colab.output.renderOutput(dataTable, element);const docLink document.createElement(div);docLink.innerHTML docLinkHtml;element.appendChild(docLink);}
/scriptsvg xmlns“http://www.w3.org/2000/svg” height24pxviewBox“0 0 24 24” width“24px”
cases.crs #This returns you the crsProjected CRS: EPSG:27700
Name: OSGB36 / British National Grid
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: United Kingdom (UK) - offshore to boundary of UKCS within 49°45N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.
- bounds: (-9.01, 49.75, 2.01, 61.01)
Coordinate Operation:
- name: British National Grid
- method: Transverse Mercator
Datum: Ordnance Survey of Great Britain 1936
- Ellipsoid: Airy 1830
- Prime Meridian: Greenwichcases.plot() #This returns you a map of all the casesAxes: pump.crsProjected CRS: EPSG:27700
Name: OSGB36 / British National Grid
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: United Kingdom (UK) - offshore to boundary of UKCS within 49°45N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.
- bounds: (-9.01, 49.75, 2.01, 61.01)
Coordinate Operation:
- name: British National Grid
- method: Transverse Mercator
Datum: Ordnance Survey of Great Britain 1936
- Ellipsoid: Airy 1830
- Prime Meridian: Greenwichpump.plot()Axes: 步骤 3
制作一个简单的地图覆盖带有泵的案例
ax cases.plot() #First we need to create an axispump.plot(axax,colororange) #Then plot the pumps using thes previously created axis as a parameter in the plot()Axes: ax cases.plot(markersize 20)
#Check out here for differernt shapes of the markers:
#https://matplotlib.org/api/markers_api.htmlpump.plot(axax,markersize500,marker*,colorred)Axes: 步骤4
使用 rasterio 读取 TIF 底图“SnowGIS/SnowMap.tif”
basemap rasterio.open(https://raw.githubusercontent.com/Ziqi-Li/GEO4162C/main/data/SnowGIS/SnowMap.tiff)show(basemap)Axes: 步骤 5
将地图叠加在一起
#First we create a figure with a size of (10,10)
f, ax plt.subplots(figsize(10,10))#Plot each layer, with the same axis ax
show(basemap,axax)cases.plot(axax, markersize 30, alpha0.8)
pump.plot(axax,markersize1000, marker*,colorred)
Axes: 步骤 6
执行核密度估计KDE
函数 sns.kdeplot() 有几个参数 GeoDataFrame cases x 和 y 坐标我们可以从 GeoDataFrame 中提取为 cases.geometry.x 和 cases.geometry.y bw_method 是 KDE 的带宽。 4.fill参数表示你想要计数图还是填充色图5.cmap参数表示您要使用的颜色 alpha参数表示透明度的级别。
基本上 4、5、6 是您可以更改的样式参数。
执行 KDE 的代码是
sns.kdeplot(datacases,xcases.geometry.x,ycases.geometry.y,bw_method0.5,fillTrue,cmapcoolwarm,alpha0.8)Axes: 让我们将 kde 覆盖在我们的地图上。 请注意您需要指定要在与其余地图相同的轴上绘制的 kdeplot。 所以我们需要在sns.kdeplot(axax,...)中添加axax。
f, ax plt.subplots(figsize(10,10))
show(basemap,axax)
cases.plot(axax, markersize 50)
pump.plot(axax,markersize1000,marker*,colorred)#4th layer
sns.kdeplot(axax, datacases,xcases.geometry.x,ycases.geometry.y,bw_method0.5,fillTrue,cmapcoolwarm,alpha0.6)Axes: 我们可以看到热点中心与中央泵对齐得很好
sns.kdeplot() 中的另一个重要参数是有一个 weights 参数您可以使用它根据每个位置的案例数量来权衡您的位置。
首先让我们回顾一下 GeoDataFrame 的情况我们发现有一列名为“Count”。然后让我们使用这个“Count”作为 KDE 中的权重。 为此请将“weightscases.Count”添加到现有的 kde 函数中
cases.head()IdCountgeometry003POINT (529308.741 181031.352)102POINT (529312.164 181025.172)201POINT (529314.382 181020.294)301POINT (529317.380 181014.259)404POINT (529320.675 181007.872) scriptconst buttonEl document.querySelector(#df-fc93c370-9178-4d6b-a535-1aca6b1bd379 button.colab-df-convert);buttonEl.style.display google.colab.kernel.accessAllowed ? block : none;async function convertToInteractive(key) {const element document.querySelector(#df-fc93c370-9178-4d6b-a535-1aca6b1bd379);const dataTable await google.colab.kernel.invokeFunction(convertToInteractive,[key], {});if (!dataTable) return;const docLinkHtml Like what you see? Visit the a target_blank hrefhttps://colab.research.google.com/notebooks/data_table.ipynbdata table notebook/a to learn more about interactive tables.;element.innerHTML ;dataTable[output_type] display_data;await google.colab.output.renderOutput(dataTable, element);const docLink document.createElement(div);docLink.innerHTML docLinkHtml;element.appendChild(docLink);}
/scriptsvg xmlns“http://www.w3.org/2000/svg” height24pxviewBox“0 0 24 24” width“24px”
f, ax plt.subplots(figsize(10,10))
show(basemap,axax,adjustNone)
cases.plot(axax, markersize 50,columnCount)
pump.plot(axax,markersize1000,marker*,colorred)sns.kdeplot(axax,datacases, xcases.geometry.x, ycases.geometry.y, bw_method0.5, fillTrue,cmapcoolwarm,alpha0.8,weightscases.Count)plt.tight_layout()plt.savefig(choloera.png,dpi600)你能找到什么 基本上 KDE 最热点正好指向中央泵 如果 John Snow 博士能够进行 KDE他的生活会变得更轻松吗
添加底图
在这里我添加了一小节关于使用“contextily”包将底图添加到绘图中的内容这在大多数情况下非常有用当我们没有另一个底图作为参考/提供必要的上下文时。
contextily 在 Google Colab 中不可用所以我们需要在这里安装它。
pip install -q contextilyimport contextily as cxf, ax plt.subplots(figsize(10,10))cases.plot(axax, markersize 50, figsize(9,9))
pump.plot(axax,markersize1000,marker*,colorred)
sns.kdeplot(axax,datacases, xcases.geometry.x, ycases.geometry.y, bw_method0.5, fillTrue,cmapcoolwarm,alpha0.8,weightscases.Count)#This is the new line of code:
#We need to put a crs into the function.
#The source defines the style of the base map: https://contextily.readthedocs.io/en/latest/providers_deepdive.html
#There are many options available.
cx.add_basemap(ax, sourcecx.providers.Stamen.TonerLite, crscases.crs) # we need to put a crs into the function.步骤 7
最近邻分析
下一节我们将使用距离函数和显着性检验进行最近邻分析。 这可以在很棒的“pointpats”包中轻松完成。
首先让我们从“cases”GeoDataFrame 中提取 x 和 y 坐标并将它们作为二维数组。
x cases.geometry.x.values
y cases.geometry.y.valuespoints np.array(list(zip(x,y)))其次让我们使用点创建一个“PointPattern()”类
pp PointPattern(points)/usr/local/lib/python3.10/dist-packages/libpysal/cg/shapes.py:1492: FutureWarning: Objects based on the Geometry class will deprecated and removed in a future version of libpysal.warnings.warn(dep_msg, FutureWarning)
/usr/local/lib/python3.10/dist-packages/libpysal/cg/shapes.py:1208: FutureWarning: Objects based on the Geometry class will deprecated and removed in a future version of libpysal.warnings.warn(dep_msg, FutureWarning)pppointpats.pointpattern.PointPattern at 0x7d59d1a5f010“knn”函数将返回每个点的最近邻居以及到该邻居的距离。
接下来让我们生成 100 个 CSR每个都是泊松过程作为我们的空基准。 请记住这将是我们进行显着性检验的置信区间。
CSRs PoissonPointProcess(pp.window, pp.n, 100, asPPTrue) # simulate CSR 100 times/usr/local/lib/python3.10/dist-packages/libpysal/cg/shapes.py:1923: FutureWarning: Objects based on the Geometry class will deprecated and removed in a future version of libpysal.warnings.warn(dep_msg, FutureWarning)
/usr/local/lib/python3.10/dist-packages/libpysal/cg/shapes.py:103: FutureWarning: Objects based on the Geometry class will deprecated and removed in a future version of libpysal.warnings.warn(dep_msg, FutureWarning)运行 G 距离 stats.Genv() 函数并将其可视化。
它显示 G 曲线高于置信包络线红色曲线表明我们在霍乱图中观察到聚集模式。
genv stats.Genv(pp, realizationsCSRs)genv.plot()您还可以轻松执行 K 或 F 距离函数。 同样K 曲线高于置信包络线红色曲线表明我们在霍乱地图中观察到聚集模式。
kenv stats.Kenv(pp, realizationsCSRs) # call Fenv for F function
kenv.plot()F 曲线低于置信包络线红色曲线表明我们在霍乱地图中观察到聚集模式。
import warnings
warnings.filterwarnings(ignore)
fenv stats.Fenv(pp, realizationsCSRs) # call Fenv for F function
fenv.plot()
