目录
热点分析 (Getis-Ord Gi*)概述
Getis-Ord Gi* 统计量是一种空间统计方法，广泛用于 识别和分析空间聚集现象，包括城市热岛效应。该方法有助于确定某一特征（如气温）在空间上的集中程度。高值聚集区域是热点，低值聚集区域是冷点。
什么是热点分析？
热点分析 (Getis-Ord Gi)* 是一种用于空间数据分析的技术，主要用于识别地理空间数据中值的聚集模式，可以帮助我们理解哪些区域存在高值或低值的聚集，这些聚集通常被称为“热点”或“冷点”，Gi* 统计量为数据集中的每个要素（例如地图上的点或区域）计算一个z得分。这个z得分可以用来判断在该位置附近是否存在显著的高值或低值聚集。
Getis-Ord Gi 原理*
热点分析 (Getis-Ord Gi*) 和高/低聚类分析 (Getis-Ord General G) 有哪些区别？
1、高/低聚类分析 (Getis-Ord General G)
全局统计量：Getis-Ord General G 是一个全局统计量，它用于检测整个研究区域内是否存在全局性的高值或低值聚类。也就是说，它提供了一个单一的度量值来说明数据集中是否存在全局的空间自相关性。
输出：General G 统计量给出一个单一的数值，用来衡量整个数据集中高值或低值是否呈现聚集的趋势。
应用场景：适合于初步检查数据中是否存在空间自相关性，或者在整个研究区域内是否有显著的高值或低值聚类。
2、热点分析 (Getis-Ord Gi)*
局部统计量：Getis-Ord Gi* 是一个局部统计量，它为数据集中的每个要素计算一个统计值。这意味着对于每一个点或区域，Gi* 统计量都会给出一个关于其周围环境是否形成热点或冷点的评估。
输出：Gi* 统计量为每个要素提供了 z 得分和 p 值，这些值用于判断某个特定位置是否为热点或冷点，以及这种聚集是否具有统计显著性。
应用场景：适用于详细地分析热点或冷点的位置以及它们的空间分布模式，帮助识别特定区域内的聚集特征。
3、关系总结
- 相似之处：两者都基于相同的基础理论——Getis-Ord 空间自相关模型，并且都使用了空间权重矩阵来定义要素间的空间关系。
- 不同之处：General G 提供了一个整体视角，而 Gi* 则提供了局部视角，使我们可以详细了解哪些地方形成了热点或冷点。
基于Python实现热点分析 (Getis-Ord Gi*)
准备工作：库包安装
Python代码实现
Python代码实现步骤如下：
1、导入栅格数据（.tif文件）。
2、计算邻接矩阵，确定每个栅格单元的邻域。
3、计算Getis-Ord Gi*统计量。
4、标准化结果并绘制热点图。
绘制图形如下：
Python代码如下：
本案例在指定的 shapefile（如 rural_shapefile）范围内计算统计量并绘制图形，需要通过地理信息工具（如 geopandas）加载 shapefile，然后基于该范围对湿度数据进行约束计算。
import os
import numpy as np
import geopandas as gpd
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from shapely.geometry import Point
from scipy.stats import zscore
# 设置全局字体为 Times New Roman
plt.rcParams["font.family"] = "Times New Roman"
# 数据文件夹路径
data_dir = r"C:\Database\3_CDMet\Relative humidity"
variable_name = "rhu" # 相对湿度变量名
output_directory = "C:/Database/3_CDMet/Figures"
# SHP 文件路径
Chinashp_file_path = "C:/PycharmProjects/GBA_Data_Process/shpFile/China/中国.shp"
# 加载中国边界数据
gdf_China = gpd.read_file(Chinashp_file_path)
def read_daily_rhu(data_dir, year, day, lat_range=None, lon_range=None):
"""
从 NetCDF 文件中读取指定年份和某一天的湿度数据，并应用缩放和无效值处理。
可根据给定的经纬度范围裁剪数据。
"""
file_path = os.path.join(data_dir, f"CDMet_rhu_{year}.nc")
if not os.path.exists(file_path):
raise FileNotFoundError(f"文件未找到: {file_path}")
with Dataset(file_path, mode="r") as nc_file:
rhu_var = nc_file.variables[variable_name]
scale_factor = getattr(rhu_var, "scale_factor", 1.0)
add_offset = getattr(rhu_var, "add_offset", 0.0)
fill_value = getattr(rhu_var, "_FillValue", -32767.0)
lat = nc_file.variables["lat"][:]
lon = nc_file.variables["lon"][:]
daily_data = rhu_var[day, :, :].astype(np.float32)
daily_data[daily_data == fill_value] = np.nan
# 裁剪经纬度范围
if lat_range:
lat_min, lat_max = lat_range
lat_idx = np.where((lat >= lat_min) & (lat <= lat_max))[0]
lat = lat[lat_idx]
daily_data = daily_data[lat_idx, :]
if lon_range:
lon_min, lon_max = lon_range
lon_idx = np.where((lon >= lon_min) & (lon <= lon_max))[0]
lon = lon[lon_idx]
daily_data = daily_data[:, lon_idx]
return daily_data, lat, lon
def apply_shapefile_mask(data, lat, lon, shapefile):
"""
使用 shapefile 约束湿度数据，将 shapefile 范围外的值设置为 NaN。
"""
# 加载 shapefile
gdf_boundary = gpd.read_file(shapefile)
boundary = gdf_boundary.unary_union # 合并为一个整体多边形
# 构造网格点坐标
lon_grid, lat_grid = np.meshgrid(lon, lat)
points = np.column_stack((lon_grid.ravel(), lat_grid.ravel()))
points_geom = [Point(pt) for pt in points]
# 检查点是否在 shapefile 范围内
mask = np.array([boundary.contains(pt) for pt in points_geom])
mask = mask.reshape(data.shape)
# 将 shapefile 范围外的值设置为 NaN
data[~mask] = np.nan
return data
def calculate_getis_ord_gi_star(data, lat, lon, dist_threshold=2.0):
"""
逐点计算 Getis-Ord Gi* 统计量，避免计算完整的距离矩阵。
"""
# 构造网格坐标
lon_grid, lat_grid = np.meshgrid(lon, lat)
points = np.column_stack((lat_grid.ravel(), lon_grid.ravel()))
data_flat = data.ravel()
# 排除 NaN 值
valid_mask = ~np.isnan(data_flat)
data_flat = data_flat[valid_mask]
points = points[valid_mask]
# 全局统计量
n = len(data_flat)
mean_x = np.mean(data_flat)
std_x = np.std(data_flat)
# 初始化 Gi* 统计量
Gi_star = np.full(n, np.nan)
for i, point in enumerate(points):
# 计算邻居：基于经纬度距离
distances = np.sqrt((points[:, 0] - point[0]) ** 2 + (points[:, 1] - point[1]) ** 2)
neighbors = distances <= dist_threshold
# 局部加权计算
local_sum = np.sum(data_flat[neighbors])
weight_sum = np.sum(neighbors)
weight_sq_sum = np.sum(neighbors ** 2)
numerator = local_sum - mean_x * weight_sum
denominator = std_x * np.sqrt((n * weight_sq_sum - weight_sum ** 2) / (n - 1))
Gi_star[i] = numerator / denominator
# 将 Gi* 统计量还原为二维网格
Gi_star_grid = np.full(lat_grid.shape, np.nan)
Gi_star_grid.ravel()[valid_mask] = Gi_star
return Gi_star_grid
def plot_hotspot(data, lat, lon, lat_range, lon_range, title):
"""
绘制 Getis-Ord Gi* 热点分析结果图。
"""
# 确保输入的 经纬度范围有效
min_lat, max_lat = lat_range
min_lon, max_lon = lon_range
fig, ax = plt.subplots(figsize=(10, 6))
hotspot_plot = ax.pcolormesh(lon, lat, data, cmap="RdBu", shading="auto")
cbar = plt.colorbar(hotspot_plot, ax=ax, fraction=0.05, pad=0.04)
cbar.set_label("Getis-Ord Gi*", fontsize=14)
ax.set_title(title, fontsize=16)
ax.set_xlabel("Longitude (°E)", fontsize=14)
ax.set_ylabel("Latitude (°N)", fontsize=14)
ax.tick_params(axis="both", which="major", labelsize=13)
gdf_China.boundary.plot(ax=ax, color="black", linewidth=1, label="China Boundary")
ax.legend()
# 设置 x 和 y 轴的范围
ax.set_xlim(min_lon, max_lon)
ax.set_ylim(min_lat, max_lat)
output_file = os.path.join(output_directory, f"{title}.png")
plt.savefig(output_file, dpi=600, bbox_inches="tight")
print(f"Exported: {output_file}")
plt.show()
plt.close()
# 主程序
if __name__ == "__main__":
# 输入文件路径
FID_obj = 33112 # 指定 FID 编号(GBA)
output_directory = f"C:/Database/3_CDMet/Result_{FID_obj}"
urban_shapefile = f"C:/Database/10 GUB-Global/GUB_Global_2015/Polygon_FID_{FID_obj}.shp" # 城市边界文件路径
rural_shapefile = f"C:/Database/10 GUB-Global/GUB_Global_2015/Rural_boundary_FID_{FID_obj}.shp" # 乡村边界文件路径
# 指定年份和天数
year = 2020
day = 0 # 0 表示 1 月 1 日
# 指定经纬度范围
lat_range = (20, 25)
lon_range = (110, 120)
# 读取湿度数据
daily_data, lat, lon = read_daily_rhu(data_dir, year, day, lat_range, lon_range)
print(f"数据形状: {daily_data.shape}")
print(f"纬度范围: {lat.min()} 至 {lat.max()}")
print(f"经度范围: {lon.min()} 至 {lon.max()}")
print(f"数据范围: 最小值={np.nanmin(daily_data)}, 最大值={np.nanmax(daily_data)}")
# 应用 shapefile 约束
daily_data = apply_shapefile_mask(daily_data, lat, lon, rural_shapefile)
# 计算不同阈值下的 Getis-Ord Gi* 统计量
Gi_star_1km = calculate_getis_ord_gi_star(daily_data, lat, lon, dist_threshold=1.0)
Gi_star_2km = calculate_getis_ord_gi_star(daily_data, lat, lon, dist_threshold=2.0)
Gi_star_5km = calculate_getis_ord_gi_star(daily_data, lat, lon, dist_threshold=5.0)
# 绘制不同阈值下的热点图
plot_hotspot(Gi_star_1km, lat, lon, lat_range, lon_range, "UDI (1km threshold)")
plot_hotspot(Gi_star_2km, lat, lon, lat_range, lon_range, "UDI (2km threshold)")
plot_hotspot(Gi_star_5km, lat, lon, lat_range, lon_range, "UDI (5km threshold)")
如何确定dist_threshold阈值？
dist_threshold 的确定原理：dist_threshold 决定了分析中点的邻域范围，其选择需根据以下原则：
1、空间分辨率：
dist_threshold 应与数据的空间分辨率匹配。例如，如果数据网格的经纬度间隔是 1°，则设置 dist_threshold = 1（或更大）可以捕获足够的邻居点。
2、研究目标：
- 如果分析局部聚集性（如城市热岛效应），选择较小的 dist_threshold。
- 如果分析大范围区域趋势，选择较大的 dist_threshold。
3、数据特性：
数据稀疏时，设置较大的 dist_threshold 以确保每个点有足够的邻居参与计算。
参考
1、CSDN博客-ArcGIS热点分析 (Getis-Ord Gi*)——基于地级市尺度的七普人口普查数据的热点与冷点分析
扫一扫
4725
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
