四个空间分析代码（M-K趋势检验、Sen-Slope趋势分析、变异系数法、多年均值），用于栅格数据输入，可输出栅格和图表。以及M-K和Sen-Slope方向比较图、研究区水系图（精度有限，但可优化）
·
输出栅格和输入栅格的地理坐标系一致。输出的图表是默认使用简易圆柱 (Plate Carrée) 投影：正方形经纬方格，1纬度长度=1经度长度的强制 1:1 坐标轴比例方格，建议根据研究区位置选择失真程度小的投影，以减少代码中默认的插值行为（用以填充经纬方格）。输出栅格如何定义投影，请询问AI，让AI对代码作出修改。
本文章代码均由AI生成，实践过有效。定义合理的投影之后可提高原先精度。
本文代码生成的图表对应像元值颜色的映射的选取存在一定主观性，为提高不同像元值间颜色的辨识度，可以改变颜色映射部分的代码。可以搜索标准色带Python代码，选取更有辨识度的色带进行颜色映射，把搜索到的代码给AI，让AI对本文代码原有的颜色映射模块这部分进行改进。
每张图片带有标题在图上方，在论文中直接使用并不规范，如有需要请用AI修改。
1、M-K趋势
import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
import rasterio
from rasterio.mask import mask
from scipy import stats
import glob
from tqdm import tqdm
import geopandas as gpd
from shapely.geometry import mapping
from concurrent.futures import ProcessPoolExecutor
from matplotlib.patches import Patch
# 设置非交互式后端，确保在服务器环境下也能生成图表
import matplotlib
matplotlib.use('Agg')
# 设置中文字体
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
# 进一步增大全局字体大小
plt.rcParams.update({
'font.size': 20, # 基础字体大小
'axes.titlesize': 28, # 标题字体大小
'axes.labelsize': 22, # 坐标轴标签字体大小
'xtick.labelsize': 20, # x轴刻度字体大小
'ytick.labelsize': 20, # y轴刻度字体大小
'legend.fontsize': 20, # 图例字体大小
'figure.titlesize': 32 # 图表标题字体大小
})
# 最小有效数据点数量
MIN_VALID_POINTS = 5
# 数据精度（保留小数位数）
DATA_PRECISION = 3
# 黑色标记值（用于所有值均为0的像元）
ZERO_VALUE_MARKER = 12
def calculate_mann_kendall(series):
"""计算Mann-Kendall趋势检验，增加数据检查和异常处理"""
# 排除NODATA值
valid_series = series[~np.isnan(series)]
# 检查有效数据点数量
if len(valid_series) < MIN_VALID_POINTS:
return np.nan, np.nan, np.nan
# 检查所有值是否相同（避免除以零）
if np.allclose(valid_series, valid_series[0], atol=1e-6):
return 0, 0, 1.0 # 无趋势，p值为1
# 计算S值
s = 0
for i in range(len(valid_series) - 1):
for j in range(i + 1, len(valid_series)):
s += np.sign(valid_series[j] - valid_series[i])
# 计算方差Var(S)
counts = np.unique(valid_series, return_counts=True)[1]
var_s = (len(valid_series) * (len(valid_series) - 1) * (2 * len(valid_series) + 5)) / 18
if len(counts) < len(valid_series): # 存在重复值
for c in counts:
var_s -= (c * (c - 1) * (2 * c + 5)) / 18
# 防止除以零
if var_s <= 0:
return 0, 0, 1.0 # 无趋势，p值为1
# 计算Z值
if s > 0:
z = (s - 1) / np.sqrt(var_s)
elif s < 0:
z = (s + 1) / np.sqrt(var_s)
else: # s == 0
z = 0
# 计算p值（双侧检验）
p = 2 * (1 - stats.norm.cdf(abs(z)))
return s, z, p
def classify_trend(z, p):
"""根据Z值和p值对趋势进行分类，处理NaN值"""
if np.isnan(z) or np.isnan(p):
return 0 # NODATA或无法计算趋势
if z == 0:
return 1 # 无趋势
# 显著性水平
if z > 0: # 增加趋势
if p < 0.001:
return 4 # 极显著增加 (p < 0.001)
elif p < 0.01:
return 5 # 高度显著增加 (p < 0.01)
elif p < 0.02:
return 6 # 显著增加 (p < 0.02)
elif p < 0.05:
return 7 # 一般显著增加 (p < 0.05)
else:
return 2 # 不显著增加
else: # 减少趋势
if p < 0.001:
return 11 # 极显著减少 (p < 0.001)
elif p < 0.01:
return 10 # 高度显著减少 (p < 0.01)
elif p < 0.02:
return 9 # 显著减少 (p < 0.02)
elif p < 0.05:
return 8 # 一般显著减少 (p < 0.05)
else:
return 3 # 不显著减少
def create_colormap():
"""创建自定义颜色映射，严格按照用户指定的颜色映射规律"""
colors = [
(0.7, 0.7, 0.7), # 1: 无趋势 - 灰色
(1.0, 1.0, 0.0), # 2: 不显著增加 - 黄色
(0.7, 1.0, 0.7), # 3: 不显著减少 - 淡绿色
(0.65, 0.33, 0.16), # 4: 极显著增加 (p < 0.001) - 褐色
(1.0, 0.0, 0.0), # 5: 高度显著增加 (p < 0.01) - 赤红色
(1.0, 0.5, 0.0), # 6: 显著增加 (p < 0.02) - 橙红色
(1.0, 0.84, 0.0), # 7: 一般显著增加 (p < 0.05) - 金色
(0.0, 1.0, 0.0), # 8: 一般显著减少 (p < 0.05) - 绿色
(0.0, 0.8, 0.8), # 9: 显著减少 (p < 0.02) - 天蓝色
(0.0, 0.0, 1.0), # 10: 高度显著减少 (p < 0.01) - 蓝色
(0.5, 0.0, 0.5), # 11: 极显著减少 (p < 0.001) - 紫罗兰色
]
cmap = LinearSegmentedColormap.from_list('trend_cmap', colors, N=11)
# 设置0（NODATA）为白色
cmap.set_under('white')
# 设置12（所有时期值均为0）为黑色，但不在色标中显示
cmap.set_over('black')
return cmap
def process_chunk(chunk_indices, all_data, zero_value_mask):
"""处理数据块的趋势分析，跳过零值标记的像元"""
height, width = all_data.shape[1:]
chunk_results = np.zeros((len(chunk_indices), width), dtype=np.uint8)
for idx, i in enumerate(chunk_indices):
for j in range(width):
# 如果是零值标记的像元，直接标记为12
if zero_value_mask[i, j]:
chunk_results[idx, j] = ZERO_VALUE_MARKER
continue
# 获取该位置的时间序列
series = all_data[:, i, j]
# 计算M-K统计量
s, z, p = calculate_mann_kendall(series)
# 分类趋势
trend_class = classify_trend(z, p)
chunk_results[idx, j] = trend_class
return chunk_indices, chunk_results
def calculate_class_percentages(trend_results, mask=None):
"""计算各趋势分类的面积占比，可选择在掩膜范围内计算"""
if mask is not None:
# 只考虑掩膜内的像素
masked_results = trend_results[mask]
else:
masked_results = trend_results.flatten()
total_pixels = len(masked_results)
if total_pixels == 0:
return {}, 0
# 统计每个分类的像素数
class_counts = np.bincount(masked_results, minlength=13)
# 计算在所有像素中的占比
percentages = {}
for class_id, count in enumerate(class_counts):
if count > 0:
percent = (count / total_pixels) * 100
percentages[class_id] = percent
return percentages, total_pixels, class_counts
def load_shapefile_mask(shapefile_path, raster_path):
"""加载SHP文件并创建掩膜"""
try:
# 读取SHP文件
if not os.path.exists(shapefile_path):
raise FileNotFoundError(f"SHP文件不存在: {shapefile_path}")
gdf = gpd.read_file(shapefile_path)
# 确保SHP文件和栅格数据投影一致
with rasterio.open(raster_path) as src:
raster_crs = src.crs
# 转换SHP文件投影
if gdf.crs != raster_crs:
print(f"警告: SHP文件投影与栅格数据不一致，正在转换...")
gdf = gdf.to_crs(raster_crs)
# 获取SHP文件的几何图形
geometries = [mapping(geom) for geom in gdf.geometry]
# 创建掩膜
out_image, out_transform = mask(src, geometries, crop=True, filled=False)
mask_array = ~out_image.mask[0] # 获取有效区域的掩膜
return mask_array
except Exception as e:
print(f"加载SHP文件时出错: {e}")
return None
def analyze_and_visualize(input_dir, output_dir, shapefile_path=None, debug_mode=False):
"""分析栅格数据并生成可视化结果，可选择在SHP文件范围内计算"""
try:
# 确保输出目录存在
os.makedirs(output_dir, exist_ok=True)
# 获取所有栅格文件
raster_files = glob.glob(os.path.join(input_dir, "*.tif"))
if not raster_files:
print("未找到栅格文件！")
return
print(f"找到 {len(raster_files)} 个栅格文件")
# 读取第一个栅格以获取元数据
first_raster_path = raster_files[0]
with rasterio.open(first_raster_path) as src:
meta = src.meta
nodata = src.nodata
height, width = src.shape
# 加载SHP文件掩膜（如果提供）
mask_array = None
if shapefile_path and os.path.exists(shapefile_path):
print(f"加载SHP文件掩膜: {shapefile_path}")
mask_array = load_shapefile_mask(shapefile_path, first_raster_path)
if mask_array is not None:
print(f"掩膜范围内有效像素占比: {np.mean(mask_array) * 100:.2f}%")
else:
print("警告: 无法加载SHP掩膜，将使用全区域计算")
else:
print("未提供SHP文件或文件不存在，将使用全区域计算")
# 一次性读取所有栅格数据
print("正在读取所有栅格数据...")
all_data = np.zeros((len(raster_files), height, width), dtype=np.float32)
for i, file in enumerate(tqdm(raster_files)):
try:
with rasterio.open(file) as src:
data = src.read(1)
# 将NODATA值转换为NaN
data[data == nodata] = np.nan
# 统一数据精度
data = np.round(data, decimals=DATA_PRECISION)
all_data[i] = data
except Exception as e:
print(f"读取文件 {file} 时出错: {e}")
all_data[i] = np.nan # 将出错的文件设为NODATA
# 识别所有24张栅格值均为0的像元
print("正在识别所有时期值均为0的像元...")
# 计算每个像元在所有时期中的和（忽略NaN）
pixel_sums = np.nansum(all_data, axis=0)
# 创建一个布尔掩码，标记所有24张栅格值均为0的像元
zero_value_mask = np.logical_and(np.isclose(pixel_sums, 0, atol=1e-6),
np.sum(~np.isnan(all_data), axis=0) == len(raster_files))
# 计算零值像元的比例
zero_value_percentage = np.mean(zero_value_mask) * 100
print(f"所有时期值均为0的像元占比: {zero_value_percentage:.2f}%")
# 检查并统计每个位置的有效数据点数量
valid_data_count = np.sum(~np.isnan(all_data), axis=0)
insufficient_data_mask = valid_data_count < MIN_VALID_POINTS
print(f"有效数据点不足({MIN_VALID_POINTS}个)的像素占比: {np.mean(insufficient_data_mask) * 100:.2f}%")
# 创建趋势结果数组
trend_results = np.zeros((height, width), dtype=np.uint8)
# 向量化处理每个栅格位置
print("正在进行M-K趋势分析...")
# 分块处理以提高内存效率
chunk_size = 100
total_chunks = (height + chunk_size - 1) // chunk_size
# 使用多进程加速
from concurrent.futures import ProcessPoolExecutor, as_completed
with ProcessPoolExecutor() as executor:
futures = []
for chunk_id in range(total_chunks):
start = chunk_id * chunk_size
end = min((chunk_id + 1) * chunk_size, height)
futures.append(executor.submit(process_chunk, range(start, end), all_data, zero_value_mask))
# 收集结果
for future in tqdm(as_completed(futures), total=total_chunks):
chunk_indices, chunk_results = future.result()
trend_results[chunk_indices] = chunk_results
# 将有效数据点不足的位置强制设为NODATA
trend_results[insufficient_data_mask] = 0
# 计算各分类占比（不排除零值像元）
print("正在计算趋势分类占比...")
percentages, total_pixels, class_counts = calculate_class_percentages(trend_results, mask_array)
# 统计SHP范围内的全零值像元数和占比
if mask_array is not None:
zero_count = np.sum(zero_value_mask & mask_array)
zero_percentage = (zero_count / np.sum(mask_array)) * 100
else:
zero_count = np.sum(zero_value_mask)
zero_percentage = (zero_count / trend_results.size) * 100
print(f"所有时期值均为0的像元数: {zero_count}")
print(f"所有时期值均为0的像元占比: {zero_percentage:.2f}%")
# 创建颜色映射（严格按照用户指定）
cmap = create_colormap()
# 可视化结果
plt.figure(figsize=(16, 12))
# 显示数据，将NODATA(0)显示为白色，值均为0的像元显示为黑色
im = plt.imshow(trend_results, cmap=cmap, interpolation='nearest', vmin=1, vmax=11)
# 添加标题
plt.title("M-K趋势分析结果", fontsize=32)
# 添加颜色条，严格按照用户指定的颜色映射和标注
cbar = plt.colorbar(im, fraction=0.05, pad=0.05, extend='neither') # 不显示扩展箭头
cbar.ax.tick_params(labelsize=20)
# 色标严格包含11个等级，且与用户提供的颜色映射规律一一对应
cbar.set_ticks(range(1, 12))
cbar.set_ticklabels([
'无趋势',
'不显著增加(p>0.05)',
'不显著减少(p>0.05)',
'极显著增加(p<0.001)',
'高度显著增加(p<0.01)',
'显著增加(p<0.02)',
'一般显著增加(p<0.05)',
'一般显著减少(p<0.05)',
'显著减少(p<0.02)',
'高度显著减少(p<0.01)',
'极显著减少(p<0.001)'
])
# 在右上角添加图例，说明黑色为所有时期值均为0的像元
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor='black', label='所有时期值均为0')]
plt.legend(handles=legend_elements, loc='upper right', frameon=True)
# 调整布局
plt.subplots_adjust(right=0.85, top=0.85)
# 移除坐标轴
plt.axis('off')
# 保存结果 - 命名为M-K修正2.png
output_path = os.path.join(output_dir, "M-K修正2.png")
plt.savefig(output_path, dpi=350, bbox_inches='tight')
print(f"结果已保存至: {output_path}")
# 保存趋势分类栅格
trend_meta = meta.copy()
trend_meta['dtype'] = 'uint8'
trend_meta['nodata'] = 0
# 栅格文件命名为mk_trend_classes_修正2.tif
trend_output_path = os.path.join(output_dir, "mk_trend_classes_修正2.tif")
with rasterio.open(trend_output_path, 'w', **trend_meta) as dst:
dst.write(trend_results, 1)
print(f"趋势分类栅格已保存至: {trend_output_path}")
# 保存百分比数据到txt文件 - 命名为M-K修正2_趋势分类占比.txt
txt_output_path = os.path.join(output_dir, "M-K修正2_趋势分类占比.txt")
with open(txt_output_path, 'w', encoding='utf-8') as f:
f.write("M-K趋势分析 - 各趋势分类面积占比\n")
f.write("=" * 60 + "\n")
if mask_array is not None:
f.write(f"分析范围: {shapefile_path}\n")
f.write(f"SHP范围内总像素数: {total_pixels}\n")
f.write(f"其中所有时期值均为0的像元数: {zero_count}\n")
f.write(f"全零值像元占比: {zero_percentage:.2f}%\n\n")
else:
f.write(f"分析范围: 全区域\n")
f.write(f"总像素数: {total_pixels}\n")
f.write(f"所有时期值均为0的像元数: {zero_count}\n")
f.write(f"全零值像元占比: {zero_percentage:.2f}%\n\n")
# 写入各分类占比和像元数
f.write("趋势分类统计（以SHP范围/全区域总像元为分母）:\n")
f.write("-" * 60 + "\n")
trend_labels = {
0: "NODATA",
1: "无趋势",
2: "不显著增加 (p>0.05)",
3: "不显著减少 (p>0.05)",
4: "极显著增加 (p<0.001)",
5: "高度显著增加 (p<0.01)",
6: "显著增加 (p<0.02)",
7: "一般显著增加 (p<0.05)",
8: "一般显著减少 (p<0.05)",
9: "显著减少 (p<0.02)",
10: "高度显著减少 (p<0.01)",
11: "极显著减少 (p<0.001)",
12: "所有时期值均为0"
}
for class_id in sorted(percentages.keys()):
label = trend_labels.get(class_id, f"未知类别 {class_id}")
count = class_counts[class_id] if class_id < len(class_counts) else 0
percent = percentages[class_id]
f.write(f"{label}: {count}个像元，占比{percent:.2f}%\n")
# 计算并写入增加和减少趋势的总占比和像元数
increasing_classes = [2, 4, 5, 6, 7]
decreasing_classes = [3, 8, 9, 10, 11]
total_increasing = sum(percentages.get(c, 0) for c in increasing_classes)
total_decreasing = sum(percentages.get(c, 0) for c in decreasing_classes)
total_increasing_count = sum(class_counts[c] if c < len(class_counts) else 0 for c in increasing_classes)
total_decreasing_count = sum(class_counts[c] if c < len(class_counts) else 0 for c in decreasing_classes)
f.write("\n")
f.write(f"总增加趋势: {total_increasing_count}个像元，占比{total_increasing:.2f}%\n")
f.write(f"总减少趋势: {total_decreasing_count}个像元，占比{total_decreasing:.2f}%\n")
# 写入显著等级的统计
f.write("\n")
f.write("各显著性等级统计:\n")
f.write("-" * 60 + "\n")
significant_classes = {
"极显著": [4, 11],
"高度显著": [5, 10],
"显著": [6, 9],
"一般显著": [7, 8],
"不显著": [2, 3]
}
for level, classes in significant_classes.items():
level_count = sum(class_counts[c] if c < len(class_counts) else 0 for c in classes)
level_percent = (level_count / total_pixels) * 100 if total_pixels > 0 else 0
f.write(f"{level}: {level_count}个像元，占比{level_percent:.2f}%\n")
print(f"趋势分类占比数据已保存至: {txt_output_path}")
except Exception as e:
print(f"执行分析时发生错误: {e}")
if __name__ == "__main__":
input_directory = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\M-K空间分析所用栅格"
output_directory = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\visualizations"
shapefile_path = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\黄河中游无定河流域" # 请替换为实际的SHP文件路径
# 启用调试模式以分析NODATA问题
debug_mode = True
analyze_and_visualize(input_directory, output_directory, shapefile_path, debug_mode)2、Sen-Slope趋势
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import geopandas as gpd
import rasterio
from rasterio.mask import mask
import logging
import numba
from datetime import datetime
from scipy import stats
# 设置中文显示
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]
plt.rcParams["axes.unicode_minus"] = False # 解决负号显示问题
# 配置日志记录
logging.basicConfig(
filename='sen_slope_analysis.log',
level=logging.ERROR,
format='%(asctime)s - %(levelname)s - %(message)s'
)
@numba.njit(parallel=True)
def calculate_sen_slope_numba(all_data, zero_mask, years):
"""使用Numba加速计算Sen's斜率"""
rows, cols = all_data.shape[1], all_data.shape[2]
slope_array = np.full((rows, cols), np.nan, dtype=np.float64)
for i in numba.prange(rows):
for j in range(cols):
if zero_mask[i, j]:
continue
pixel_values = all_data[:, i, j]
valid_count = 0
for t in range(len(pixel_values)):
if not np.isnan(pixel_values[t]):
valid_count += 1
if valid_count < 3:
continue
slopes = []
for m in range(len(pixel_values)):
if np.isnan(pixel_values[m]):
continue
for n in range(m + 1, len(pixel_values)):
if np.isnan(pixel_values[n]):
continue
s = (pixel_values[n] - pixel_values[m]) / (years[n] - years[m])
slopes.append(s)
if slopes:
slopes_array = np.array(slopes)
slope = np.median(slopes_array)
slope_array[i, j] = slope
return slope_array
def extract_year_from_filename(filename):
"""从文件名中提取年份信息"""
try:
# 假设文件名格式为 'MOD16A3GF.A2000001.ET_500m.tif.tif.tif'
parts = filename.split('.')
if len(parts) < 3:
raise ValueError(f"文件名格式不符合预期: {filename}")
a_part = parts[1] # 例如 'A2000001'
if not a_part.startswith('A'):
raise ValueError(f"文件名格式不符合预期: {filename}")
year_str = a_part[1:5] # 提取年份部分
return int(year_str)
except Exception as e:
logging.error(f"从文件名 {filename} 提取年份时出错: {str(e)}")
raise
def calculate_sen_slope(raster_files, shapefile, zero_mask, shapefile_mask):
"""计算Sen's斜率"""
# 修改：接收mask_raster_with_shapefile的四个返回值
first_raster, _, meta, _ = mask_raster_with_shapefile(raster_files[0], shapefile)
rows, cols = first_raster.shape[1], first_raster.shape[2]
# 从文件名中提取年份并排序
years = []
file_year_pairs = []
for raster_path in raster_files:
filename = os.path.basename(raster_path)
year = extract_year_from_filename(filename)
years.append(year)
file_year_pairs.append((year, raster_path))
# 按年份排序
file_year_pairs.sort(key=lambda x: x[0])
sorted_years = [pair[0] for pair in file_year_pairs]
sorted_files = [pair[1] for pair in file_year_pairs]
print(f"已识别并排序的年份: {sorted_years}")
# 准备用于计算斜率的年份索引（以年为单位）
year_indices = np.array([y - sorted_years[0] for y in sorted_years], dtype=np.float64)
print(f"开始计算Sen's斜率，共{len(sorted_files)}个时期")
all_data = np.full((len(sorted_files), rows, cols), np.nan, dtype=np.float64)
for t, raster_path in enumerate(sorted_files):
try:
with rasterio.open(raster_path) as src:
data = src.read(1, masked=True)
data = np.ma.filled(data, np.nan)
all_data[t] = data
except Exception as e:
logging.error(f"读取栅格文件 {raster_path} 时出错: {str(e)}")
# 只在shapefile范围内处理数据
all_data[:, ~shapefile_mask] = np.nan
all_data[:, zero_mask] = np.nan
slope_array = calculate_sen_slope_numba(all_data, zero_mask, year_indices)
# 只在shapefile范围内保留计算结果
slope_array[~shapefile_mask] = np.nan
print(f"Sen's斜率计算完成")
classified_array = classify_slope(slope_array, zero_mask)
return slope_array, meta, classified_array
def classify_slope(slope_array, zero_mask):
"""对斜率进行分类"""
classified_array = np.full_like(slope_array, np.nan, dtype=np.float32)
non_zero_mask = ~zero_mask
# 根据新的分类规则进行分类
classified_array[non_zero_mask & (slope_array >= -4) & (slope_array < -2)] = -2
classified_array[non_zero_mask & (slope_array >= -2) & (slope_array < 0)] = -1
classified_array[non_zero_mask & (slope_array == 0)] = 0
classified_array[non_zero_mask & (slope_array > 0) & (slope_array < 2)] = 1
classified_array[non_zero_mask & (slope_array >= 2) & (slope_array < 4)] = 2
classified_array[non_zero_mask & (slope_array >= 4) & (slope_array < 6)] = 3
classified_array[non_zero_mask & (slope_array >= 6) & (slope_array < 8)] = 4
classified_array[non_zero_mask & (slope_array >= 8) & (slope_array < 10)] = 5
classified_array[non_zero_mask & (slope_array >= 10) & (slope_array < 12)] = 6
classified_array[non_zero_mask & (slope_array >= 12) & (slope_array < 14)] = 7
classified_array[non_zero_mask & (slope_array >= 14) & (slope_array < 16)] = 8
classified_array[non_zero_mask & (slope_array >= 16) & (slope_array <= 18)] = 9
classified_array[zero_mask] = np.nan
return classified_array
def create_colorbar():
"""创建自定义色标"""
# 调整颜色映射和标签
colors = [
'#E6E6FA', # 薰衣草色 (-4 to -2)
'#FF00FF', # 品红色 (-2 to 0)
'#FF0000', # 红色 (0 to 2)
'#8B4513', # 赭色 (2 to 4)
'#FFA500', # 橙子色 (4 to 6)
'#F0E68C', # 浅金色 (6 to 8)
'#BDB76B', # 暗黄褐色 (8 to 10)
'#32CD32', # 柠檬绿色 (10 to 12)
'#00FFFF', # 青色 (12 to 14)
'#00BFFF', # 深天蓝色 (14 to 16)
'#9370DB', # 中紫色 (16 to 18)
]
# 调整色标边界
bounds = [-4, -2, 0, 2, 4, 6, 8, 10, 12, 14, 16, 18]
# 创建自定义颜色映射
cmap = mcolors.LinearSegmentedColormap.from_list('custom_cmap', colors, N=len(colors))
cmap.set_bad(color='white') # 设置NaN值为白色
norm = mcolors.BoundaryNorm(bounds, cmap.N)
return cmap, norm, bounds
def calculate_statistics(slope_array, zero_mask, classified_array, shapefile_mask):
"""计算统计数据，分母只包含shp文件范围内的像元"""
print(f"[统计计算] shapefile范围内像元数: {np.sum(shapefile_mask)}")
print(f"[统计计算] 零值像元掩码形状: {zero_mask.shape}, 唯一值: {np.unique(zero_mask)}")
print(f"[统计计算] 斜率数组形状: {slope_array.shape}, NaN数量: {np.sum(np.isnan(slope_array))}")
# 只在shapefile范围内计算统计数据
slope_in_shape = slope_array[shapefile_mask]
zero_in_shape = zero_mask[shapefile_mask]
classified_in_shape = classified_array[shapefile_mask]
# 计算不同区域的像元数
total_pixels = np.sum(shapefile_mask)
nan_pixels = np.sum(np.isnan(slope_in_shape))
valid_pixels = total_pixels - nan_pixels
# 计算零值像元在shapefile范围内的数量
total_zero_pixels = np.sum(zero_in_shape)
# 计算零值像元在有效区域中的数量
valid_mask = ~np.isnan(slope_in_shape)
valid_zero_pixels = np.sum(zero_in_shape & valid_mask)
# 计算零值像元在NaN区域中的数量
nan_zero_pixels = np.sum(zero_in_shape & ~valid_mask)
print(f"[统计计算] shapefile范围内总像元数: {total_pixels}")
print(f"[统计计算] shapefile范围内NaN像元数: {nan_pixels}")
print(f"[统计计算] shapefile范围内有效像元数(非NaN): {valid_pixels}")
print(f"[统计计算] shapefile范围内24年均为0值的像元数: {total_zero_pixels}")
print(f"[统计计算] shapefile范围内24年均为0值的像元数(有效区域): {valid_zero_pixels}")
print(f"[统计计算] shapefile范围内24年均为0值的像元数(NaN区域): {nan_zero_pixels}")
# 使用shapefile范围内的所有像元作为分母（包括NaN值）
denominator = total_pixels
# 计算各区间像元数
category_counts = {
-2: np.sum(classified_in_shape == -2),
-1: np.sum(classified_in_shape == -1),
0: np.sum(classified_in_shape == 0),
1: np.sum(classified_in_shape == 1),
2: np.sum(classified_in_shape == 2),
3: np.sum(classified_in_shape == 3),
4: np.sum(classified_in_shape == 4),
5: np.sum(classified_in_shape == 5),
6: np.sum(classified_in_shape == 6),
7: np.sum(classified_in_shape == 7),
8: np.sum(classified_in_shape == 8),
9: np.sum(classified_in_shape == 9),
}
# 计算各区间占比，使用shapefile范围内的所有像元作为分母
category_percentages = {cat: (count / denominator) * 100 for cat, count in category_counts.items()}
# 计算零值像元占比
zero_percentage = (total_zero_pixels / denominator) * 100 if denominator > 0 else 0
# 计算占比总和（排除NaN值）
total_percentage = sum(category_percentages.values())
# 计算斜率最大值和最小值
valid_slope = slope_in_shape[valid_mask & ~zero_in_shape]
min_slope = np.min(valid_slope) if len(valid_slope) > 0 else np.nan
max_slope = np.max(valid_slope) if len(valid_slope) > 0 else np.nan
# 计算中位数和众数
slope_median = np.median(valid_slope) if len(valid_slope) > 0 else np.nan
slope_mode = stats.mode(valid_slope, keepdims=True)[0][0] if len(valid_slope) > 0 else np.nan
return {
'total_pixels': total_pixels,
'nan_pixels': nan_pixels,
'valid_pixels': valid_pixels,
'total_zero_pixels': total_zero_pixels,
'valid_zero_pixels': valid_zero_pixels,
'nan_zero_pixels': nan_zero_pixels,
'denominator': denominator,
'zero_percentage': zero_percentage,
'category_counts': category_counts,
'category_percentages': category_percentages,
'total_percentage': total_percentage,
'min_slope': min_slope,
'max_slope': max_slope,
'median_slope': slope_median,
'mode_slope': slope_mode,
}
def save_statistics_to_txt(statistics, output_path):
"""将统计数据保存到txt文件"""
with open(output_path, 'w', encoding='utf-8') as f:
f.write(f"Sen's斜率分析统计结果\n")
f.write(f"分析日期: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n")
f.write("=" * 70 + "\n")
f.write(f"shapefile范围内总像元数: {statistics['total_pixels']}\n")
f.write(f"其中NaN像元数: {statistics['nan_pixels']}\n")
f.write(f"有效像元数(非NaN): {statistics['valid_pixels']}\n")
f.write(f"shapefile范围内24年均为0值的像元数: {statistics['total_zero_pixels']}\n")
f.write(f"其中，在有效区域中的零值像元数: {statistics['valid_zero_pixels']}\n")
f.write(f"其中，在NaN区域中的零值像元数: {statistics['nan_zero_pixels']}\n")
f.write("=" * 70 + "\n")
f.write(f"统计占比计算使用的分母: {statistics['denominator']} (shapefile范围内包含24年全为0的NaN值像元)\n")
f.write(f"零值像元占比: {statistics['zero_percentage']:.4f}%\n")
f.write("=" * 70 + "\n")
f.write(f"Sen's斜率最小值: {statistics['min_slope']:.4f}\n")
f.write(f"Sen's斜率最大值: {statistics['max_slope']:.4f}\n")
f.write(f"Sen's斜率中位数: {statistics['median_slope']:.4f}\n")
f.write(f"Sen's斜率众数: {statistics['mode_slope']:.4f}\n")
f.write("=" * 70 + "\n")
f.write("各斜率区间占比统计（分母包含shapefile范围内的NaN值）:\n")
category_labels = {
-2: '[-4, -2)',
-1: '[-2, 0)',
0: '[0, 0]',
1: '(0, 2)',
2: '[2, 4)',
3: '[4, 6)',
4: '[6, 8)',
5: '[8, 10)',
6: '[10, 12)',
7: '[12, 14)',
8: '[14, 16)',
9: '[16, 18]'
}
for cat in sorted(category_labels.keys()):
count = statistics['category_counts'][cat]
percent = statistics['category_percentages'][cat]
f.write(f"{category_labels[cat]:<10}: {count:>10d} 个像元 ({percent:.4f}%)\n")
f.write("=" * 70 + "\n")
f.write(f"所有区间占比总和: {statistics['total_percentage']:.4f}%\n")
f.write("=" * 70 + "\n")
def read_raster_files(folder_path):
"""读取文件夹中的所有栅格文件"""
raster_files = []
for file in os.listdir(folder_path):
if file.endswith(('.tif', '.TIF')):
raster_files.append(os.path.join(folder_path, file))
return sorted(raster_files)
def read_shapefile(shp_path):
"""读取shapefile文件"""
for file in os.listdir(shp_path):
if file.endswith('.shp'):
return gpd.read_file(os.path.join(shp_path, file))
raise FileNotFoundError("未找到shapefile文件")
def mask_raster_with_shapefile(raster_path, shapefile, dtype=np.float32):
"""使用shapefile对栅格进行掩膜处理"""
try:
with rasterio.open(raster_path) as src:
# 获取原始栅格的元数据
original_meta = src.meta
# 使用shapefile裁剪栅格
geoms = shapefile.geometry.values
out_image, out_transform = mask(src, geoms, crop=True, nodata=np.nan)
out_meta = src.meta
out_image = out_image.astype(dtype)
# 记录裁剪后的栅格尺寸
rows, cols = out_image.shape[1], out_image.shape[2]
# 创建一个掩码，表示shapefile范围内的区域
# 任何在裁剪后栅格中非NaN的位置都被认为是shapefile范围内
shapefile_mask = ~np.isnan(out_image[0])
# 计算shapefile范围内的像元数
shapefile_pixels = np.sum(shapefile_mask)
print(f"[掩膜处理] 原始栅格尺寸: {original_meta['width']} x {original_meta['height']}")
print(f"[掩膜处理] 裁剪后栅格尺寸: {cols} x {rows}")
print(f"[掩膜处理] shapefile范围内的像元数: {shapefile_pixels}")
print(
f"[掩膜处理] shapefile范围外的像元数: {original_meta['width'] * original_meta['height'] - shapefile_pixels}")
out_meta.update({
"driver": "GTiff",
"height": rows,
"width": cols,
"transform": out_transform,
"dtype": dtype
})
return out_image, out_transform, out_meta, shapefile_mask
except Exception as e:
logging.error(f"处理栅格文件 {raster_path} 时出错: {str(e)}")
raise
def identify_zero_pixels(raster_files, shapefile, tolerance=1e-6, output_debug_folder=None):
"""识别所有时期均为0值的像元，增加详细调试信息和可视化"""
print(f"开始识别24年均为0值的像元，共{len(raster_files)}个时期")
# 确保调试输出文件夹存在
if output_debug_folder:
os.makedirs(output_debug_folder, exist_ok=True)
# 获取第一个栅格的尺寸和掩码，同时获取shapefile_mask
first_raster, _, _, shapefile_mask = mask_raster_with_shapefile(raster_files[0], shapefile)
rows, cols = first_raster.shape[1], first_raster.shape[2]
# 创建一个全True的掩码，表示所有像元初始都被认为是0值
zero_mask = np.ones((rows, cols), dtype=bool)
# 记录每个栅格的基本统计信息
grid_stats = []
# 可视化设置
if output_debug_folder:
plt.figure(figsize=(15, 10))
for i, raster_path in enumerate(raster_files):
try:
with rasterio.open(raster_path) as src:
# 使用shapefile裁剪栅格
geoms = shapefile.geometry.values
out_image, _ = mask(src, geoms, crop=True, nodata=np.nan)
raster_data = out_image[0] # 获取第一个波段
# 统计基本信息
file_stats = {
'filename': os.path.basename(raster_path),
'min_value': np.nanmin(raster_data),
'max_value': np.nanmax(raster_data),
'mean_value': np.nanmean(raster_data),
'nan_count': np.sum(np.isnan(raster_data)),
'zero_count_exact': np.sum(raster_data == 0),
'zero_count_tolerance': np.sum(np.abs(raster_data) < tolerance)
}
# 识别当前栅格中的0值像元（使用容差）
current_zero = np.abs(raster_data) < tolerance
# 只在shapefile范围内更新zero_mask
current_zero[~shapefile_mask] = False
# 保存当前时期的零值像元图用于调试
if output_debug_folder:
plt.clf()
plt.imshow(current_zero, cmap='binary')
plt.title(f"时期 {i + 1}: {os.path.basename(raster_path)}\n"
f"零值像元数: {np.sum(current_zero)}")
plt.colorbar(label='是否为零值')
plt.savefig(os.path.join(output_debug_folder, f"zero_pixels_period_{i + 1}.png"), dpi=300)
# 更新zero_mask，只保留所有时期都为0的像元
zero_mask = np.logical_and(zero_mask, current_zero)
# 更新文件统计信息
file_stats['cumulative_zero_count'] = np.sum(zero_mask)
grid_stats.append(file_stats)
print(f"处理文件 {i + 1}/{len(raster_files)}: {os.path.basename(raster_path)}")
print(f" 最小值: {file_stats['min_value']:.6f}, 最大值: {file_stats['max_value']:.6f}")
print(f" 精确0值像元数: {file_stats['zero_count_exact']}")
print(f" 容差范围内0值像元数: {file_stats['zero_count_tolerance']}")
print(f" 累积24年均为0值的像元数: {file_stats['cumulative_zero_count']}")
except Exception as e:
logging.error(f"识别零值像元时处理文件 {raster_path} 出错: {str(e)}")
raise
# 只在shapefile范围内保留零值像元
zero_mask[~shapefile_mask] = False
# 保存最终的零值像元图
if output_debug_folder and np.sum(zero_mask) > 0:
plt.clf()
plt.imshow(zero_mask, cmap='binary')
plt.title(f"24年均为0值的像元分布\n总数: {np.sum(zero_mask)}")
plt.colorbar(label='是否为零值')
plt.savefig(os.path.join(output_debug_folder, "final_zero_pixels.png"), dpi=300)
# 输出总体统计信息
total_zero_pixels = np.sum(zero_mask)
print("=" * 50)
print(f"最终24年均为0值的像元数: {total_zero_pixels}")
print(f"占shapefile范围内总面积比例: {total_zero_pixels / np.sum(shapefile_mask) * 100:.4f}%")
# 如果没有找到零值像元，输出详细统计信息帮助调试
if total_zero_pixels == 0:
print("\n警告: 未找到24年均为0值的像元")
print("各栅格文件统计信息:")
for stats in grid_stats:
print(f"\n文件: {stats['filename']}")
print(f" 最小值: {stats['min_value']:.6f}, 最大值: {stats['max_value']:.6f}")
print(f" NaN数量: {stats['nan_count']}")
print(f" 精确0值像元数: {stats['zero_count_exact']}")
print(f" 容差范围内0值像元数: {stats['zero_count_tolerance']}")
# 检查是否有接近零的值
near_zero = np.logical_and(raster_data > 0, raster_data < 1e-3)
near_zero_count = np.sum(near_zero)
if near_zero_count > 0:
print(f" 接近零的值(0-1e-3)数量: {near_zero_count}")
print(f" 最接近零的非零值: {np.nanmin(raster_data[raster_data > 0]):.10f}")
return zero_mask, shapefile_mask
def find_key_statistics_locations(slope_array, zero_mask, shapefile_mask, statistics):
"""找到斜率极大值、极小值、中位数和众数像元的位置"""
print("正在查找关键统计值的位置...")
# 只考虑shapefile范围内且非零值的有效像元
valid_mask = shapefile_mask & ~zero_mask
valid_slope = slope_array[valid_mask]
# 找到极大值和极小值的位置
max_slope = statistics['max_slope']
min_slope = statistics['min_slope']
median_slope = statistics['median_slope']
mode_slope = statistics['mode_slope']
# 创建标记数组
max_locations = np.zeros_like(slope_array, dtype=bool)
min_locations = np.zeros_like(slope_array, dtype=bool)
median_locations = np.zeros_like(slope_array, dtype=bool)
mode_locations = np.zeros_like(slope_array, dtype=bool)
# 由于浮点数精度问题，使用容差比较
tolerance = 1e-6
# 找到极大值位置
max_locations[valid_mask] = np.abs(slope_array[valid_mask] - max_slope) < tolerance
# 找到极小值位置
min_locations[valid_mask] = np.abs(slope_array[valid_mask] - min_slope) < tolerance
# 找到中位数位置
median_locations[valid_mask] = np.abs(slope_array[valid_mask] - median_slope) < tolerance
# 找到众数位置
mode_locations[valid_mask] = np.abs(slope_array[valid_mask] - mode_slope) < tolerance
print(f"找到 {np.sum(max_locations)} 个极大值像元")
print(f"找到 {np.sum(min_locations)} 个极小值像元")
print(f"找到 {np.sum(median_locations)} 个中位数像元")
print(f"找到 {np.sum(mode_locations)} 个众数像元")
return max_locations, min_locations, median_locations, mode_locations
def main():
"""主函数，执行完整的Sen's斜率分析流程"""
# 设置路径
raster_folder = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\M-K空间分析所用栅格"
shp_folder = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\黄河中游无定河流域"
output_folder = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\slope趋势分析"
reference_raster_path = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\results\bijiaoM-K修正1.tif"
# 确保输出文件夹存在
os.makedirs(output_folder, exist_ok=True)
print("开始Sen's斜率分析...")
# 读取参考栅格尺寸
print("正在读取参考栅格尺寸...")
try:
with rasterio.open(reference_raster_path) as ref_src:
reference_width = ref_src.width
reference_height = ref_src.height
reference_transform = ref_src.transform
print(f"参考栅格尺寸: {reference_width} x {reference_height}")
except Exception as e:
logging.error(f"读取参考栅格时出错: {str(e)}")
print(f"警告: 无法读取参考栅格，将使用默认尺寸")
reference_width = None
reference_height = None
reference_transform = None
# 读取数据
print("正在读取栅格数据...")
raster_files = read_raster_files(raster_folder)
if not raster_files:
print("未找到栅格文件！")
return
print(f"找到 {len(raster_files)} 个栅格文件")
print("正在读取shapefile...")
shapefile = read_shapefile(shp_folder)
print("正在识别24年均为0值的像元...")
debug_folder = os.path.join(output_folder, "zero_pixels_debug")
zero_mask, shapefile_mask = identify_zero_pixels(raster_files, shapefile, tolerance=1e-6,
output_debug_folder=debug_folder)
print(f"[主函数] 识别出的24年均为0值的像元数: {np.sum(zero_mask)}")
print(f"[主函数] shapefile范围内像元数: {np.sum(shapefile_mask)}")
# 计算Sen's斜率
print("正在计算Sen's斜率...")
slope_array, meta, classified_array = calculate_sen_slope(raster_files, shapefile, zero_mask, shapefile_mask)
print(f"[主函数] 计算完成后零值像元数: {np.sum(zero_mask)}")
# 如果读取到参考栅格尺寸，调整输出栅格尺寸
if reference_width and reference_height:
print(f"调整输出栅格尺寸为: {reference_width} x {reference_height}")
# 创建一个新的空数组，大小为参考栅格尺寸
new_slope_array = np.full((reference_height, reference_width), np.nan, dtype=np.float64)
new_classified_array = np.full((reference_height, reference_width), np.nan, dtype=np.float32)
new_zero_mask = np.full((reference_height, reference_width), False, dtype=bool)
new_shapefile_mask = np.full((reference_height, reference_width), False, dtype=bool)
# 计算可以复制的区域
min_rows = min(slope_array.shape[0], reference_height)
min_cols = min(slope_array.shape[1], reference_width)
# 将计算结果复制到新数组中
new_slope_array[:min_rows, :min_cols] = slope_array[:min_rows, :min_cols]
new_classified_array[:min_rows, :min_cols] = classified_array[:min_rows, :min_cols]
new_zero_mask[:min_rows, :min_cols] = zero_mask[:min_rows, :min_cols]
new_shapefile_mask[:min_rows, :min_cols] = shapefile_mask[:min_rows, :min_cols]
# 更新元数据
meta.update({
'width': reference_width,
'height': reference_height,
'transform': reference_transform
})
slope_array = new_slope_array
classified_array = new_classified_array
zero_mask = new_zero_mask
shapefile_mask = new_shapefile_mask # 确保更新后的shapefile_mask被使用
print(f"[主函数] 调整尺寸后零值像元数: {np.sum(zero_mask)}")
print(f"[主函数] 调整尺寸后shapefile范围内像元数: {np.sum(shapefile_mask)}")
# 计算统计数据
print("正在计算统计数据...")
statistics = calculate_statistics(slope_array, zero_mask, classified_array, shapefile_mask)
print(f"[主函数] 统计数据中的零值像元数: {statistics['total_zero_pixels']}")
# 保存统计数据
stats_file = os.path.join(output_folder, "slopeXX6.txt")
save_statistics_to_txt(statistics, stats_file)
print(f"统计数据已保存至: {stats_file}")
# 查找关键统计值的位置
max_locations, min_locations, median_locations, mode_locations = find_key_statistics_locations(
slope_array, zero_mask, shapefile_mask, statistics
)
# 创建主图表
print("正在创建主图表...")
cmap, norm, bounds = create_colorbar()
plt.figure(figsize=(18, 16))
# 创建一个背景层，显示shapefile范围
background = np.ma.masked_where(shapefile_mask, np.ones_like(slope_array))
plt.imshow(background, cmap='gray', alpha=0.2)
# 显示斜率图
im = plt.imshow(slope_array, cmap=cmap, norm=norm, interpolation='nearest')
# 显示黑色的全0像元
zero_display = np.ma.masked_where(~zero_mask, zero_mask)
plt.imshow(zero_display, cmap='Greys', alpha=0.8, vmin=0, vmax=1)
# 添加颜色条
cbar = plt.colorbar(im, ticks=bounds, pad=0.03, shrink=0.8)
# 调整色标标题位置和字体
cbar.set_label('Sen\'s斜率', fontsize=30, labelpad=5, rotation=90)
cbar_label = cbar.ax.get_yaxis().get_label()
cbar_label.set_verticalalignment('bottom')
cbar_label.set_position((1.0, 0.5))
cbar_label.set_horizontalalignment('center')
# 修改色标刻度标签
tick_labels = ['-4', '-2', '0(不含)', '2', '4', '6', '8', '10', '12', '14', '16', '18']
cbar.ax.set_yticklabels(tick_labels, fontsize=24)
# 去掉横纵坐标轴
plt.axis('off')
# 添加右上角图例
from matplotlib.patches import Patch
legend_elements = [
Patch(facecolor='#808080', label='24年均为0值'),
Patch(facecolor='#D3D3D3', label='shapefile范围外')
]
plt.legend(handles=legend_elements, loc='upper right', fontsize=24,
frameon=True, facecolor='white', edgecolor='black', framealpha=0.8,
handlelength=1.5, handleheight=1.5, borderpad=0.5, labelspacing=0.5)
plt.title('Sen\'s斜率分析结果', fontsize=36, pad=30)
plt.tight_layout()
# 保存主图表
chart_file = os.path.join(output_folder, "slopeXX6.png")
plt.savefig(chart_file, dpi=300, bbox_inches='tight')
print(f"主图表已保存至: {chart_file}")
# 创建标注关键统计值的图表
print("正在创建标注关键统计值的图表...")
plt.figure(figsize=(18, 16))
# 显示背景层和斜率图
plt.imshow(background, cmap='gray', alpha=0.2)
# 重新创建图像对象，避免引用上一个图表的图像
im2 = plt.imshow(slope_array, cmap=cmap, norm=norm, interpolation='nearest', alpha=0.7)
# 显示黑色的全0像元
plt.imshow(zero_display, cmap='Greys', alpha=0.8, vmin=0, vmax=1)
# 标注极大值位置 (DarkMagenta)
max_display = np.ma.masked_where(~max_locations, np.ones_like(slope_array))
plt.imshow(max_display, cmap=mcolors.ListedColormap(['DarkMagenta']), alpha=0.8)
# 标注极小值位置 (DarkCyan)
min_display = np.ma.masked_where(~min_locations, np.ones_like(slope_array))
plt.imshow(min_display, cmap=mcolors.ListedColormap(['DarkCyan']), alpha=0.8)
# 标注中位数位置 (DarkBlue)
median_display = np.ma.masked_where(~median_locations, np.ones_like(slope_array))
plt.imshow(median_display, cmap=mcolors.ListedColormap(['DarkBlue']), alpha=0.8)
# 标注众数位置 (LightGreen)
mode_display = np.ma.masked_where(~mode_locations, np.ones_like(slope_array))
plt.imshow(mode_display, cmap=mcolors.ListedColormap(['LightGreen']), alpha=0.8)
# 添加颜色条，使用新创建的图像对象
cbar = plt.colorbar(im2, ticks=bounds, pad=0.03, shrink=0.8)
cbar.set_label('Sen\'s斜率', fontsize=30, labelpad=5, rotation=90)
cbar.ax.set_yticklabels(tick_labels, fontsize=24)
# 去掉横纵坐标轴
plt.axis('off')
# 添加图例
legend_elements = [
Patch(facecolor='#808080', label='24年均为0值'),
Patch(facecolor='#D3D3D3', label='shapefile范围外'),
Patch(facecolor='DarkMagenta',
label=f'极大值点 ({np.sum(max_locations)}个, 值: {statistics["max_slope"]:.4f})'),
Patch(facecolor='DarkCyan', label=f'极小值点 ({np.sum(min_locations)}个, 值: {statistics["min_slope"]:.4f})'),
Patch(facecolor='DarkBlue',
label=f'中位数点 ({np.sum(median_locations)}个, 值: {statistics["median_slope"]:.4f})'),
Patch(facecolor='LightGreen', label=f'众数点 ({np.sum(mode_locations)}个, 值: {statistics["mode_slope"]:.4f})'),
]
plt.legend(handles=legend_elements, loc='upper right', fontsize=20,
frameon=True, facecolor='white', edgecolor='black', framealpha=0.8,
handlelength=1.5, handleheight=1.5, borderpad=0.5, labelspacing=0.5)
# 在像元旁边添加标签
print("正在添加标签...")
# 定义标签样式
label_styles = {
'max': {'color': 'DarkMagenta', 'prefix': 'Max: '},
'min': {'color': 'DarkCyan', 'prefix': 'Min: '},
'median': {'color': 'DarkBlue', 'prefix': 'Median: '},
'mode': {'color': 'LightGreen', 'prefix': 'Mode: '}
}
# 收集所有需要标注的点
all_points = []
# 添加最大值点
max_y, max_x = np.nonzero(max_locations)
for x, y in zip(max_x, max_y):
all_points.append({'x': x, 'y': y, 'value': statistics['max_slope'], 'type': 'max'})
# 添加最小值点
min_y, min_x = np.nonzero(min_locations)
for x, y in zip(min_x, min_y):
all_points.append({'x': x, 'y': y, 'value': statistics['min_slope'], 'type': 'min'})
# 添加中位数点
median_y, median_x = np.nonzero(median_locations)
for x, y in zip(median_x, median_y):
all_points.append({'x': x, 'y': y, 'value': statistics['median_slope'], 'type': 'median'})
# 添加众数点
mode_y, mode_x = np.nonzero(mode_locations)
for x, y in zip(mode_x, mode_y):
all_points.append({'x': x, 'y': y, 'value': statistics['mode_slope'], 'type': 'mode'})
# 简单的标签位置优化算法
# 定义标签的大小和最小间距
label_width = 100
label_height = 30
min_distance = label_width * 1.5
# 已经放置的标签位置
placed_labels = []
# 尝试为每个点找到合适的位置
for point in all_points:
# 定义8个可能的位置：上、下、左、右、左上、右上、左下、右下
possible_positions = [
{'dx': 0, 'dy': -label_height - 10}, # 上
{'dx': 0, 'dy': 10}, # 下
{'dx': -label_width - 10, 'dy': 0}, # 左
{'dx': 10, 'dy': 0}, # 右
{'dx': -label_width - 10, 'dy': -label_height - 10}, # 左上
{'dx': 10, 'dy': -label_height - 10}, # 右上
{'dx': -label_width - 10, 'dy': 10}, # 左下
{'dx': 10, 'dy': 10} # 右下
]
# 随机打乱位置顺序，增加多样性
np.random.shuffle(possible_positions)
# 找到最佳位置
best_position = None
for pos in possible_positions:
label_x = point['x'] + pos['dx']
label_y = point['y'] + pos['dy']
# 检查是否与已放置的标签冲突
conflict = False
for placed in placed_labels:
dist = np.sqrt((label_x - placed['x']) ** 2 + (label_y - placed['y']) ** 2)
if dist < min_distance:
conflict = True
break
if not conflict:
best_position = pos
break
# 如果没有找到不冲突的位置，选择第一个位置
if best_position is None:
best_position = possible_positions[0]
# 计算最终标签位置
label_x = point['x'] + best_position['dx']
label_y = point['y'] + best_position['dy']
# 记录已放置的标签
placed_labels.append({'x': label_x, 'y': label_y})
# 添加标签
style = label_styles[point['type']]
plt.annotate(
f"{style['prefix']}{point['value']:.4f}",
(point['x'], point['y']),
xytext=(label_x, label_y),
textcoords='data',
arrowprops=dict(arrowstyle="->", color=style['color']),
bbox=dict(boxstyle="round,pad=0.3", fc=style['color'], ec='black', alpha=0.7),
color='white',
fontsize=12,
ha='center',
va='center'
)
plt.title('Sen\'s斜率分析 - 关键统计值位置', fontsize=36, pad=30)
plt.tight_layout()
# 保存标注图表
marked_chart_file = os.path.join(output_folder, "sloppeXX6增加图.png")
plt.savefig(marked_chart_file, dpi=300, bbox_inches='tight')
print(f"标注关键统计值的图表已保存至: {marked_chart_file}")
# 保存分类后的栅格
print("正在保存分类后的栅格...")
classified_file = os.path.join(output_folder, "classified_slopeXX6.tif")
meta.update({"dtype": "float32", "nodata": np.nan})
with rasterio.open(classified_file, 'w', **meta) as dst:
dst.write(classified_array, 1)
print("分析完成！")
if __name__ == "__main__":
main() 3、变异系数法
import os
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
import rasterio
from rasterio.mask import mask
from matplotlib.colors import LinearSegmentedColormap
import pandas as pd
# 设置中文显示
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]
plt.rcParams["axes.unicode_minus"] = False # 解决负号显示问题
plt.rcParams["font.size"] = 32 # 增大基础字体大小
def calculate_cv(input_dir, shp_path, output_dir):
"""
计算指定shp文件范围内栅格数据的变异系数
参数:
input_dir: 栅格文件所在目录
shp_path: 矢量边界文件路径
output_dir: 结果输出目录
"""
# 确保输出目录存在
os.makedirs(output_dir, exist_ok=True)
# 读取shp文件
gdf = gpd.read_file(shp_path)
# 确保shp文件有CRS
if gdf.crs is None:
print("警告: shp文件没有定义CRS，假设为WGS84 (EPSG:4326)")
gdf = gdf.set_crs(epsg=4326)
# 存储所有栅格数据的数组
raster_arrays = []
profile = None
# 读取所有栅格文件
raster_files = [f for f in os.listdir(input_dir) if f.endswith(('.tif', '.tiff'))]
if not raster_files:
print("未找到栅格文件")
return
# 统计成功和失败的栅格数量
success_count = 0
failed_count = 0
# 优化：只在需要时转换shp的CRS
shp_crs_transformed = False
gdf_transformed = None
for raster_file in raster_files:
raster_path = os.path.join(input_dir, raster_file)
try:
with rasterio.open(raster_path) as src:
# 获取栅格的CRS
raster_crs = src.crs
# 检查shp和栅格的CRS是否一致
if gdf.crs != raster_crs:
if not shp_crs_transformed:
print(
f"警告: {raster_file} 的CRS ({raster_crs}) 与shp的CRS ({gdf.crs}) 不匹配，正在转换shp的CRS")
gdf_transformed = gdf.to_crs(raster_crs)
shp_crs_transformed = True
else:
if gdf_transformed is None:
gdf_transformed = gdf
# 裁剪栅格到shp范围
out_image, out_transform = mask(src, gdf_transformed.geometry, crop=True, nodata=np.nan)
out_meta = src.meta
# 更新元数据
out_meta.update({
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform,
"nodata": np.nan
})
# 检查裁剪结果是否有效
valid_pixels = np.sum(~np.isnan(out_image))
if valid_pixels == 0:
print(f"警告: {raster_file} 裁剪后全部为NaN值，可能是由于shp和栅格没有重叠区域")
failed_count += 1
continue
# 存储栅格数据
raster_arrays.append(out_image[0]) # 假设是单波段栅格
# 保存第一个栅格的元数据用于输出
if profile is None:
profile = out_meta
success_count += 1
print(f"成功处理栅格: {raster_file}，有效像素: {valid_pixels}")
except Exception as e:
print(f"处理栅格 {raster_file} 时出错: {e}")
failed_count += 1
continue
# 输出处理统计信息
print(f"栅格处理完成: 成功 {success_count} 个, 失败 {failed_count} 个")
if not raster_arrays:
print("没有成功读取任何栅格数据")
return
# 转换为numpy数组
raster_data = np.array(raster_arrays)
# 创建掩码：识别所有时间步长中值都为0的像元
all_zero_mask = np.all(raster_data == 0, axis=0)
print(f"识别到所有时间步长中值都为0的像元数量: {np.sum(all_zero_mask)}")
# 创建NaN值掩码（shp范围外的区域）
nan_mask = np.isnan(raster_data[0])
# 计算0值像元占总面积的比例（只考虑shp范围内的像元）
total_pixels = np.sum(~nan_mask)
zero_percentage = np.sum(all_zero_mask) / total_pixels * 100
print(f"0值像元占shp范围内总面积的比例: {zero_percentage:.2f}%")
# 创建掩码，排除0值像元和NaN值像元
valid_mask = ~all_zero_mask & ~nan_mask
# 检查是否有足够的有效数据点
valid_pixels_per_location = np.sum(~np.isnan(raster_data), axis=0)
min_valid_pixels = 2 # 至少需要两个有效数据点来计算标准差
valid_locations = valid_pixels_per_location >= min_valid_pixels
# 组合掩码，确保0值像元不参与CV计算
final_mask = valid_mask & valid_locations
if not np.any(final_mask):
print(f"错误: 没有位置有足够的有效数据点 (至少需要{min_valid_pixels}个) 来计算变异系数")
return
# 计算均值和标准差，仅考虑有足够有效数据且非0的位置
masked_raster = np.ma.masked_array(raster_data, mask=np.isnan(raster_data) | ~final_mask)
mean_data = np.ma.mean(masked_raster, axis=0).filled(np.nan)
std_data = np.ma.std(masked_raster, axis=0, ddof=1).filled(np.nan) # 样本标准差
# 计算变异系数，避免除以零
cv_data = np.zeros_like(mean_data)
valid_divisions = (mean_data != 0) & (std_data != 0) & final_mask
cv_data[valid_divisions] = std_data[valid_divisions] / mean_data[valid_divisions]
# 保存CV栅格
cv_output_path = os.path.join(output_dir, "cv修正1_result.tif")
with rasterio.open(cv_output_path, 'w', **profile) as dst:
dst.write(cv_data, 1)
# 保存0值掩码
zero_mask_path = os.path.join(output_dir, "cv修正1_zero_mask.tif")
with rasterio.open(zero_mask_path, 'w', **profile) as dst:
dst.write(all_zero_mask.astype(np.uint8), 1)
# 生成图表
plot_cv(cv_data, all_zero_mask, nan_mask, profile, gdf.to_crs(profile['crs']), output_dir)
# 生成统计数据txt
generate_stats(cv_data, all_zero_mask, nan_mask, output_dir, zero_percentage)
print("处理完成!")
def plot_cv(cv_data, zero_mask, nan_mask, profile, gdf, output_dir):
"""
生成CV空间分布图，标记0值像元为黑色，shp范围外为白色
参数:
cv_data: CV值数组
zero_mask: 全0值像元掩码
nan_mask: NaN值（shp范围外）掩码
profile: 栅格元数据
gdf: 矢量边界数据
output_dir: 输出目录
"""
# 创建自定义颜色映射，从CV=0.05开始的渐变
colors = [(0.87, 0.0, 0.0), # CV=0.05: 赤红
(1.0, 0.64, 0.0), # 橙色
(1.0, 1.0, 0.0), # 黄色
(0.5, 0.9, 0.5), # 黄绿色
(0.0, 0.75, 0.75), # 水色
(0.0, 0.7, 1.0), # 深天蓝色
(0.8, 0.6, 0.8), # 蓟色
(0.75, 0.75, 0.75)] # CV=0.35: 银色
# 创建颜色映射，将颜色映射到0.05-0.35区间
cmap = LinearSegmentedColormap.from_list('custom_cmap', colors, N=256)
cmap.set_bad(color='white') # 设置NaN值显示为白色
# 计算栅格数据的宽高比
height, width = cv_data.shape
aspect_ratio = width / height
# 创建图表，根据栅格比例调整大小
fig_width = 24 # 固定宽度
fig_height = fig_width / aspect_ratio # 根据比例计算高度
fig = plt.figure(figsize=(fig_width, fig_height))
# 使用transform属性手动计算边界
transform = profile['transform']
left = transform.c
top = transform.f
right = transform.c + width * transform.a
bottom = transform.f + height * transform.e
# 创建单一显示数组：
# - shp范围外的区域(nan_mask)设为NaN（显示为白色）
# - 0值像元设为-1（后续通过cmap的under属性映射为黑色）
# - 其他值保持原有CV值
display_data = np.copy(cv_data)
display_data[nan_mask] = np.nan
display_data[zero_mask] = -1 # 特殊值，用于映射为黑色
# 创建栅格图
ax = fig.add_axes([0.1, 0.15, 0.7, 0.7]) # 调整位置
# 设置颜色映射的下限为黑色
cmap.set_under(color='black')
# 绘制图像，将颜色映射到0.05-0.35区间
im = ax.imshow(display_data, cmap=cmap, extent=[left, right, bottom, top],
vmin=0.05, vmax=0.35, aspect='auto')
# 确保图表中只显示裁剪区域内的内容
ax.set_xlim(left, right)
ax.set_ylim(bottom, top)
# 去掉横纵坐标轴和图框
ax.axis('off')
# 设置标题
fig.suptitle("CV空间分布图", fontsize=48, y=0.92)
# 添加色标
cbar_ax = fig.add_axes([0.85, 0.15, 0.03, 0.7]) # 调整色标位置和宽度
# 明确指定色标范围和刻度
cbar = fig.colorbar(im, cax=cbar_ax, ticks=[0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35])
# 调整色标标签朝向向外
cbar.set_label('变异系数', rotation=90, labelpad=20, fontsize=36)
cbar.ax.tick_params(labelsize=32) # 增大色标刻度字体
# 保存图表
plt.savefig(os.path.join(output_dir, "cv空间修正1.png"), dpi=300, bbox_inches='tight')
plt.close()
def generate_stats(cv_data, zero_mask, nan_mask, output_dir, zero_percentage):
"""
生成CV值统计数据txt和条状图，排除0值像元和shp范围外的区域，只分六类并排除占比为0的类别
参数:
cv_data: CV值数组
zero_mask: 全0值像元掩码
nan_mask: NaN值（shp范围外）掩码
output_dir: 输出目录
zero_percentage: 0值像元占比
"""
# 创建排除0值像元、NaN值像元的有效CV数据
valid_mask = ~zero_mask & ~nan_mask
valid_cv = cv_data[valid_mask & ~np.isnan(cv_data)]
if valid_cv.size == 0:
print("没有有效数据用于统计")
return
# 计算最大最小值
min_cv = np.min(valid_cv)
max_cv = np.max(valid_cv)
# 确定六类区间，确保覆盖整个范围
interval = (max_cv - min_cv) / 6
bins = [min_cv + i * interval for i in range(7)]
bin_labels = [f"{bins[i]:.3f}-{bins[i + 1]:.3f}" for i in range(6)]
# 计算各区间占比
hist, _ = np.histogram(valid_cv, bins=bins)
percentages = hist / valid_cv.size * 100
# 创建DataFrame，排除占比为0的类别
stats_df = pd.DataFrame({
'CV区间': bin_labels,
'像素数量': hist,
'占比(%)': percentages
})
# 过滤掉占比为0的类别
stats_df = stats_df[stats_df['占比(%)'] > 0]
# 如果所有类别都被过滤掉，添加一个默认类别
if stats_df.empty:
stats_df = pd.DataFrame({
'CV区间': ['全部范围'],
'像素数量': [valid_cv.size],
'占比(%)': [100.0]
})
# 计算不含0的占比总和
non_zero_percentage_sum = np.sum(stats_df['占比(%)'])
# 添加统计摘要
summary = pd.DataFrame({
'统计项': ['最小值', '最大值', '平均值', '中位数', '总像素数', '0值像素占比', '不含0的占比总和'],
'值': [min_cv, max_cv, np.mean(valid_cv), np.median(valid_cv), valid_cv.size, f"{zero_percentage:.2f}%",
f"{non_zero_percentage_sum:.2f}%"]
})
# 保存到txt
output_path = os.path.join(output_dir, "cv统计修正1.txt")
with open(output_path, 'w', encoding='utf-8') as f:
f.write("CV值区间统计 (不包含0值像元)\n" + "=" * 40 + "\n")
f.write(stats_df.to_string(index=False) + "\n\n")
f.write("统计摘要\n" + "=" * 40 + "\n")
f.write(summary.to_string(index=False))
print(f"统计数据已保存到 {output_path}")
# 生成条状图
plt.figure(figsize=(24, 16))
# 绘制直方图，使用与主图表完全相同的颜色映射
fig, ax = plt.subplots(figsize=(24, 16))
# 创建与空间分布图相同的颜色映射
colors = [(0.87, 0.0, 0.0), # 赤红
(1.0, 0.64, 0.0), # 橙色
(1.0, 1.0, 0.0), # 黄色
(0.5, 0.9, 0.5), # 黄绿色
(0.0, 0.75, 0.75), # 水色
(0.0, 0.7, 1.0), # 深天蓝色
(0.8, 0.6, 0.8), # 蓟色
(0.75, 0.75, 0.75)] # 银色
# 计算每个区间的中心值
bin_centers = [(bins[i] + bins[i + 1]) / 2 for i in range(6)]
# 为每个区间分配对应的颜色
# 首先将颜色映射到0.05-0.35的范围
color_map = LinearSegmentedColormap.from_list('custom_cmap', colors, N=256)
norm = plt.Normalize(0.05, 0.35)
# 为每个区间选择颜色
bar_colors = []
for center in bin_centers:
# 将中心值归一化到0-1范围
normalized_value = (center - 0.05) / (0.35 - 0.05)
normalized_value = max(0, min(1, normalized_value)) # 确保在0-1范围内
bar_colors.append(color_map(normalized_value))
# 绘制条形图，使用与主图表完全对应的颜色
bars = ax.bar(stats_df['CV区间'], stats_df['占比(%)'], color=[bar_colors[i] for i in stats_df.index],
edgecolor='black')
# 添加数据标签，保留三位小数
for bar in bars:
height = bar.get_height()
ax.text(bar.get_x() + bar.get_width() / 2., height + 0.05,
f'{height:.3f}%', ha='center', va='bottom', fontsize=36)
# 设置图表标题和坐标轴标签
ax.set_title('CV值分布直方图 (不包含0值像元)', fontsize=52)
ax.set_xlabel('CV区间', fontsize=40)
ax.set_ylabel('占比 (%)', fontsize=40)
# 设置x轴标签旋转
ax.tick_params(axis='x', rotation=45, labelsize=36)
ax.tick_params(axis='y', labelsize=36)
# 将注释放在右上角
annotation = f'注:\n0值像元占研究区\n总面积的 {zero_percentage:.2f}%\n不含0的占比总和为 {non_zero_percentage_sum:.2f}%'
ax.text(0.98, 0.98, annotation, ha='right', va='top', fontsize=36,
transform=ax.transAxes, bbox={"facecolor": "white", "alpha": 0.8, "pad": 10},
linespacing=1.5)
# 调整布局
plt.tight_layout()
# 保存图表
plt.savefig(os.path.join(output_dir, "cv分布修正1.png"), dpi=200, bbox_inches='tight')
plt.close()
if __name__ == "__main__":
# 输入输出路径
input_dir = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\M-K空间分析所用栅格"
shp_path = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\黄河中游无定河流域"
output_dir = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\变异系数空间分析"
# 执行分析
calculate_cv(input_dir, shp_path, output_dir)4、多年均值
import os
import glob
import time
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
import rasterio
from rasterio.mask import mask
from matplotlib.colors import LinearSegmentedColormap
# 设置中文字体和全局字体大小
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]
plt.rcParams["axes.unicode_minus"] = False # 解决负号显示问题
plt.rcParams["font.size"] = 28 # 增大字体
def process_raster_data(raster_folder, shp_path, output_folder):
"""
处理栅格数据，计算多年均值，生成图表和统计报告
参数:
raster_folder (str): 栅格数据文件夹路径
shp_path (str): 矢量边界文件路径
output_folder (str): 输出结果文件夹路径
"""
# 确保输出文件夹存在
os.makedirs(output_folder, exist_ok=True)
print("开始处理栅格数据...")
start_time = time.time()
# 读取矢量边界
print("读取矢量边界文件...")
gdf = gpd.read_file(shp_path)
# 获取所有栅格文件路径
raster_files = glob.glob(os.path.join(raster_folder, "*.tif"))
if not raster_files:
print("未找到栅格文件！")
return
print(f"找到 {len(raster_files)} 个栅格文件")
# 存储所有栅格数据的数组
all_data = []
profile = None
# 循环处理每个栅格文件
for i, raster_file in enumerate(raster_files):
try:
print(f"处理栅格文件 {i + 1}/{len(raster_files)}: {os.path.basename(raster_file)}")
# 读取栅格数据
with rasterio.open(raster_file) as src:
# 裁剪栅格数据到矢量边界内
out_image, out_transform = mask(src, gdf.geometry, crop=True)
out_meta = src.meta
# 保存第一个栅格的配置信息
if profile is None:
profile = out_meta.copy()
# 获取栅格数据（去掉波段维度）
data = out_image[0]
# 数据类型转换和异常值处理
# 使用float64避免计算过程中溢出
data = data.astype(np.float64)
# 检查是否存在无限值
if np.isinf(data).any():
print(f"警告: 文件 {raster_file} 包含无限值，已替换为 NaN")
data[np.isinf(data)] = np.nan
# 跳过全是无效值的栅格
if np.all(np.isnan(data)):
print(f"警告: 文件 {raster_file} 全是无效值，已跳过")
continue
all_data.append(data)
except Exception as e:
print(f"处理文件 {raster_file} 时出错: {e}")
if not all_data:
print("没有有效的栅格数据！")
return
print(f"栅格数据处理完成，耗时: {time.time() - start_time:.2f} 秒")
# 识别24年栅格均为0值的像元
if len(all_data) == 24: # 确保有24年数据
print("识别24年栅格均为0值的像元...")
all_data_array = np.array(all_data)
# 创建全0掩码（所有年份该位置都为0）
all_zero_mask = np.all(all_data_array == 0, axis=0)
print(f"识别到 {np.sum(all_zero_mask)} 个像元在所有年份中均为0值")
else:
print(f"数据年份数为 {len(all_data)}，不等于24。将不进行全0值像元的特殊处理")
all_zero_mask = np.zeros_like(all_data[0], dtype=bool)
# 创建有效数据掩码（排除全0值但包含shp范围外的像元）
valid_data_mask = np.ones_like(all_zero_mask, dtype=bool)
valid_data_mask[all_zero_mask] = False
# 创建shp范围内的掩码
shp_mask = ~np.all(all_data_array == profile['nodata'], axis=0)
# 计算多年均值（仅针对有效数据且非全0值）
print("计算多年均值...")
# 使用float64进行计算以避免溢出
masked_data = np.ma.array(all_data, mask=[~valid_data_mask for _ in range(len(all_data))], dtype=np.float64)
mean_data = np.ma.mean(masked_data, axis=0).filled(np.nan)
# 检查计算结果是否包含inf或NaN
if np.isinf(mean_data).any() or np.isnan(mean_data).all():
print("警告: 多年均值计算结果包含无限值或全部为NaN，尝试调整计算方法...")
# 尝试逐像元计算均值，避免一次性计算导致的溢出
mean_data = np.zeros_like(all_data[0], dtype=np.float64)
count = np.zeros_like(all_data[0], dtype=np.int32)
for data in all_data:
valid_pixels = ~np.isnan(data) & valid_data_mask
mean_data[valid_pixels] += data[valid_pixels]
count[valid_pixels] += 1
# 避免除以0
valid_count = count > 0
mean_data[valid_count] /= count[valid_count]
mean_data[~valid_count] = np.nan
# 再次检查结果
if np.isinf(mean_data).any() or np.isnan(mean_data).all():
print("错误: 调整计算方法后仍无法得到有效结果！")
return
# 版本号
version = "新新1"
# 保存多年均值栅格
print("保存多年均值栅格...")
mean_output_path = os.path.join(output_folder, f"多年均值{version}.tif")
profile.update({"driver": "GTiff", "height": mean_data.shape[0],
"width": mean_data.shape[1], "transform": out_transform})
# 保存为float32以节省空间
with rasterio.open(mean_output_path, "w", **profile) as dst:
dst.write(mean_data.astype(np.float32), 1)
# 自定义颜色映射
colors = [
(1.0, 0.75, 0.8), # 粉色 (200)
(0.58, 0.29, 0.58), # 中紫色
(0.29, 0.59, 0.82), # 宝蓝色
(0.25, 0.88, 0.82), # 绿松石色
(0.4, 0.9, 0.4), # 草坪绿
(1.0, 1.0, 0.0), # 黄色
(1.0, 0.39, 0.28), # 番茄色
(0.8, 0.0, 0.0) # 赤红色 (500)
]
max_value = 500 # 颜色映射的最大值
positions = [0.0, 0.15, 0.3, 0.5, 0.6, 0.75, 0.85, 1.0]
# 创建自定义颜色映射
custom_cmap = LinearSegmentedColormap.from_list('custom_cmap', list(zip(positions, colors)))
custom_cmap.set_bad('white') # shp范围外显示为白色
custom_cmap.set_under('black') # 小于vmin的值显示为黑色
# 生成图表（多年均值新新1）
print("生成空间分布图...")
# 调整图表尺寸，取消超高分辨率
fig, ax = plt.subplots(figsize=(20, 18)) # 减小图表尺寸
# 调整子图布局
plt.subplots_adjust(left=0.10, right=0.85, top=0.88, bottom=0.10)
# 创建用于显示的数组
mean_data_display = mean_data.copy()
# 将shp范围外的区域设为NaN，确保显示为白色
mean_data_display[~shp_mask] = np.nan
# 将24年均为0值的像元设为199（略小于vmin），确保显示为黑色
mean_data_display[all_zero_mask] = 199
# 整体缩小显示，但保持完整数据
scale_factor = 0.94
# 创建一个新的数组，用于存储缩小后的数据
height, width = mean_data.shape
new_height = int(height * scale_factor)
new_width = int(width * scale_factor)
# 创建一个更大的白色背景
display_data = np.full_like(mean_data, np.nan)
# 计算居中放置的位置
y_start = (height - new_height) // 2
x_start = (width - new_width) // 2
# 使用scipy.ndimage.zoom进行缩放
try:
from scipy.ndimage import zoom
# 对数据进行缩放
zoomed_data = zoom(mean_data_display, scale_factor, order=0)
# 确保缩放后的数据尺寸与目标位置匹配
zoomed_height, zoomed_width = zoomed_data.shape
display_height = min(new_height, zoomed_height)
display_width = min(new_width, zoomed_width)
display_data[y_start:y_start + display_height, x_start:x_start + display_width] = zoomed_data[:display_height,
:display_width]
except ImportError:
# 如果没有scipy，使用简单的切片方法
print("警告: 未找到scipy库，将使用简单方法进行缩放")
y_step = max(1, int(1 / scale_factor))
x_step = max(1, int(1 / scale_factor))
sliced_height = (height // y_step)
sliced_width = (width // x_step)
display_height = min(new_height, sliced_height)
display_width = min(new_width, sliced_width)
display_data[y_start:y_start + display_height, x_start:x_start + display_width] = mean_data_display[
:display_height * y_step:y_step,
:display_width * x_step:x_step]
# 设置颜色映射和值范围
im = ax.imshow(display_data, cmap=custom_cmap, vmin=200, vmax=max_value)
# 自定义颜色条
cbar = plt.colorbar(im, ax=ax, label='均值', pad=0.03, shrink=0.7)
cbar.ax.set_frame_on(True)
cbar.outline.set_linewidth(2)
cbar.outline.set_color('black')
# 调整字体大小
cbar.ax.tick_params(labelsize=28)
cbar.set_label('均值', fontsize=32, labelpad=15)
# 在颜色条上添加特殊值标记
cbar.set_ticks([200, 250, 300, 350, 400, 450, 500])
# 设置标题
ax.set_title('多年均值空间分布', fontsize=40, pad=20)
# 移除边框和轴线
ax.axis('off')
for spine in ax.spines.values():
spine.set_visible(False)
plt.tight_layout()
# 保存图像，降低DPI以提高性能
chart_path = os.path.join(output_folder, f"多年均值{version}.png")
plt.savefig(chart_path, dpi=300, bbox_inches='tight')
plt.close()
print(f"空间分布图已保存到: {chart_path}")
# 仅保留一个版本的图像，移除其他高分辨率版本
# 计算统计数据（仅针对shp范围内且非0值的像元）
print("计算统计数据...")
stats_data = mean_data[shp_mask & valid_data_mask]
# 计算总像元数（shp范围内除去24年均为0值的像元数）
total_valid_pixels = np.sum(shp_mask & valid_data_mask)
# 计算0值像元数和占比（仅针对shp范围内）
zero_pixels = np.sum(shp_mask & all_zero_mask)
total_shp_pixels = np.sum(shp_mask)
zero_percentage = zero_pixels / total_shp_pixels * 100
# 检查数据是否有效
if stats_data.size == 0:
print("处理后没有有效数据！")
return
# 确定合理的间隔
min_val = np.min(stats_data)
max_val = np.max(stats_data)
# 根据数据范围确定分箱数量和区间
if max_val - min_val > 20:
num_bins = 10
else:
num_bins = min(10, int(max_val - min_val) + 1)
# 确保区间分割合理
bin_edges = np.linspace(min_val, max_val, num_bins + 1)
# 计算直方图
hist, _ = np.histogram(stats_data, bins=bin_edges)
# 计算每个区间的像元占比（基于shp范围内总像元数，不含24年均为0值的像元），保留三位小数
percentages = hist / total_valid_pixels * 100
# 创建数据框
df = pd.DataFrame({
'值区间': [f'{bin_edges[i]:.2f}-{bin_edges[i + 1]:.2f}' for i in range(len(bin_edges) - 1)],
'像元数量': hist,
'占比(%)': percentages
})
# 添加统计摘要
stats = pd.DataFrame({
'统计量': ['最小值', '最大值', '平均值', '中位数', '标准差'],
'值': [np.min(stats_data), np.max(stats_data), np.mean(stats_data),
np.median(stats_data), np.std(stats_data)]
})
# 生成占比图（直方图）
print("生成直方图...")
try:
# 调整图表尺寸，使其更宽更高
plt.figure(figsize=(32, 18)) # 宽度从25增加到32，高度保持18
# 创建子图，设置边框
ax = plt.subplot(111)
# 计算每个柱子的中心值，用于颜色映射
bin_centers = [(bin_edges[i] + bin_edges[i + 1]) / 2 for i in range(len(bin_edges) - 1)]
# 归一化中心值到0-1范围，用于颜色映射
norm = plt.Normalize(vmin=200, vmax=max_value)
# 为每个柱子分配颜色
colors = [custom_cmap(norm(center)) for center in bin_centers]
# 绘制直方图，使用占比作为纵坐标
bars = plt.bar(range(len(hist)), percentages, color=colors, edgecolor='black', linewidth=1.5)
# 更新标题和标签，增大字体
plt.title('均值分布直方图', fontsize=56, pad=30) # 增大标题字体
plt.xlabel('均值区间', fontsize=50, labelpad=25) # 增大x轴标签字体
plt.ylabel('占比(%)', fontsize=50, labelpad=25) # 增大y轴标签字体
# 设置x轴刻度和标签，调整字体大小和旋转角度
plt.xticks(range(len(hist)), [f'{bin_edges[i]:.0f}-{bin_edges[i + 1]:.0f}' for i in range(len(bin_edges) - 1)],
rotation=30, ha='center', fontsize=42) # 增大x轴刻度字体
plt.yticks(fontsize=42) # 增大y轴刻度字体
# 计算最大柱子高度的10%作为注释的垂直偏移量
offset = max(percentages) * 0.08
# 添加数值标签，保留三位小数，将标签放在柱子上方但更靠近柱子顶部，增大字体
for bar in bars:
height = bar.get_height()
plt.text(bar.get_x() + bar.get_width() / 2., height + offset,
f'{height:.3f}%', ha='center', va='bottom', fontsize=46) # 增大标签字体
# 设置子图边框
for spine in ax.spines.values():
spine.set_visible(True)
spine.set_linewidth(2.5) # 增加边框线宽
# 调整y轴上限，确保有足够空间放置标签
plt.ylim(0, max(percentages) * 1.15) # 调整y轴上限，留出15%的空间
# 调整子图位置，留出更多空间
plt.subplots_adjust(left=0.08, right=0.98, top=0.88, bottom=0.15)
plt.tight_layout()
# 保存直方图
hist_path = os.path.join(output_folder, f"均值分布直方图{version}.png")
plt.savefig(hist_path, dpi=300, bbox_inches='tight')
plt.close()
print(f"直方图已成功保存到: {hist_path}")
except Exception as e:
print(f"生成直方图时出错: {e}")
print("直方图生成失败，但其他处理继续进行...")
# 保存到文本文件
print("保存统计信息到文本文件...")
txt_path = os.path.join(output_folder, f"栅格统计信息{version}.txt")
with open(txt_path, 'w', encoding='utf-8') as f:
f.write("=" * 50 + "\n")
f.write(f"栅格数据统计分析报告（仅shp范围内）{version}\n")
f.write("=" * 50 + "\n\n")
f.write("数据摘要:\n")
f.write(f"总栅格文件数: {len(raster_files)}\n")
f.write(f"有效栅格文件数: {len(all_data)}\n")
f.write(f"shp分析区域总像元数: {total_shp_pixels}\n")
f.write(f"shp范围内24年均为0值的像元数: {zero_pixels}\n")
f.write(f"24年均为0值像元占shp文件范围内总像元数的比值: {zero_percentage:.3f}%\n\n")
f.write("统计量（排除shp范围内24年均为0值的像元）:\n")
for _, row in stats.iterrows():
f.write(f"{row['统计量']}: {row['值']:.4f}\n")
f.write("\n")
f.write("各均值区间像元分布（以shp文件范围内除去24年均为0值的像元数为总数）:\n")
f.write("{:<15} {:<15} {:<15}\n".format("均值区间", "像元数量", "占比(%)"))
for _, row in df.iterrows():
f.write("{:<15} {:<15} {:<15.3f}\n".format(
row['值区间'], row['像元数量'], row['占比(%)']))
# 添加占比总和
total_percentage = np.sum(percentages)
f.write("-" * 45 + "\n")
f.write("{:<15} {:<15} {:<15.3f}\n".format(
"总计", total_valid_pixels, total_percentage))
print(f"处理完成！结果已保存到: {output_folder}")
print(f"总耗时: {time.time() - start_time:.2f} 秒")
if __name__ == "__main__":
# 输入参数
raster_folder = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\M-K空间分析所用栅格"
shp_path = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\黄河中游无定河流域"
output_folder = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\多年均值"
# 执行处理
process_raster_data(raster_folder, shp_path, output_folder) 5、M-K趋势与Sen-Slope趋势计算后的趋势方向（增加或减少）的重合比较
该代码需要用到前两段M-K趋势与Sen-Slope趋势代码生成的栅格用于比较
import os
import numpy as np
import rasterio
from rasterio.mask import mask
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.colors import LinearSegmentedColormap
import geopandas as gpd
from shapely.geometry import mapping
def compare_rasters(mk_raster_path, slope_raster_path, shp_path, output_dir="."):
"""比较两个栅格文件并生成多种比较图表"""
try:
# 设置中文字体和全局字体大小
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]
plt.rcParams.update({
'font.size': 16, # 基础字体大小
'axes.titlesize': 18, # 标题字体大小
'axes.labelsize': 16, # 坐标轴标签字体大小
'xtick.labelsize': 14, # x轴刻度字体大小
'ytick.labelsize': 14, # y轴刻度字体大小
'legend.fontsize': 16, # 图例字体大小
'figure.titlesize': 20 # 图表标题字体大小
})
# 确保输出目录存在
os.makedirs(output_dir, exist_ok=True)
# 读取栅格数据
with rasterio.open(mk_raster_path) as mk_src, rasterio.open(slope_raster_path) as slope_src:
# 确保两个栅格具有相同的形状和地理参考
if mk_src.shape != slope_src.shape or mk_src.transform != slope_src.transform:
raise ValueError("两个栅格的形状或地理参考不匹配")
# 读取栅格数据
mk_data = mk_src.read(1)
slope_data = slope_src.read(1)
# 创建有效数据掩码（排除无效值）
valid_mask = np.logical_and(~np.isnan(mk_data), ~np.isnan(slope_data))
# 计算正趋势重合、负趋势重合、零趋势重合和不重合的比例
pos_trend = np.logical_and(mk_data > 0, slope_data > 0)
neg_trend = np.logical_and(mk_data < 0, slope_data < 0)
zero_trend = np.logical_and(mk_data == 0, slope_data == 0)
# 修复：分两步使用logical_and组合多个条件
non_match = np.logical_and(~pos_trend, ~neg_trend)
non_match = np.logical_and(non_match, ~zero_trend)
non_match = np.logical_and(non_match, valid_mask)
total_valid = np.count_nonzero(valid_mask)
percent_pos = (np.count_nonzero(pos_trend) / total_valid) * 100
percent_neg = (np.count_nonzero(neg_trend) / total_valid) * 100
percent_zero = (np.count_nonzero(zero_trend) / total_valid) * 100
percent_non_match = (np.count_nonzero(non_match) / total_valid) * 100
# 创建比较结果栅格
comparison = np.zeros_like(mk_data, dtype=np.int16)
comparison[pos_trend] = 1 # 正趋势重合
comparison[neg_trend] = 2 # 负趋势重合
comparison[zero_trend] = 3 # 零趋势重合
comparison[non_match] = 0 # 不重合
comparison[~valid_mask] = -9999 # 无效值
# 新增：计算M-K像元值为5和Slope像元值为5-9的重合区域
mk_value_5 = mk_data == 5
slope_value_5_9 = np.logical_and(slope_data >= 5, slope_data <= 9)
overlap_mask = np.logical_and(mk_value_5, slope_value_5_9)
# 计算区域面积（以像元数表示）
mk_5_area = np.count_nonzero(mk_value_5)
slope_5_9_area = np.count_nonzero(slope_value_5_9)
overlap_area = np.count_nonzero(overlap_mask)
# 计算M-K=5与Slope=5-9重叠区域占Slope=5-9区域的百分比
if slope_5_9_area > 0:
overlap_percent_of_slope = (overlap_area / slope_5_9_area) * 100
else:
overlap_percent_of_slope = 0.0
# 计算重合百分比（原逻辑保留）
if total_valid > 0:
overlap_percentage = (overlap_area / total_valid) * 100
else:
overlap_percentage = 0.0
# 创建重合比较结果栅格
overlap_comparison = np.zeros_like(mk_data, dtype=np.int16)
overlap_comparison[overlap_mask] = 1 # 重合区域
overlap_comparison[~valid_mask] = -9999 # 无效值
# 如果提供了shapefile，使用它来掩膜结果
if shp_path and os.path.exists(shp_path):
# 读取shapefile
gdf = gpd.read_file(shp_path)
# 确保shapefile是WGS84坐标系
gdf = gdf.to_crs(mk_src.crs)
# 创建掩膜
geom = gdf.geometry.values[0]
# 将几何对象转换为GeoJSON格式
geom_geojson = [mapping(geom)]
# 掩膜栅格
out_image, out_transform = mask(mk_src, geom_geojson, crop=True)
mask_data = out_image[0]
# 获取掩膜的有效区域
mask_valid = mask_data != mk_src.nodata
# 将掩膜外的区域设置为无效值
comparison[~mask_valid] = -9999
overlap_comparison[~mask_valid] = -9999
# 创建图表1：原始趋势比较
chart_path = os.path.join(output_dir, "M-Kslope比较新1.png")
# 创建新布局的图表 (2行2列)
fig = plt.figure(figsize=(15, 15))
# 1. M-K趋势栅格
ax1 = plt.subplot(221)
# 修改颜色映射，M-K无趋势改为黄色
cmap_mk = LinearSegmentedColormap.from_list(
'mk_cmap',
[(0, 'blue'), (0.5, 'yellow'), (1, 'red')]
)
im1 = ax1.imshow(mk_data, cmap=cmap_mk, vmin=-1, vmax=1)
ax1.set_title('M-K趋势栅格')
ax1.axis('off') # 去掉坐标轴
# 创建M-K趋势图例 - 调整位置更靠近显示区域
mk_legend_labels = {
'red': 'M-K增加',
'blue': 'M-K减少',
'yellow': 'M-K无趋势'
}
mk_patches = [mpatches.Patch(color=color, label=label) for color, label in mk_legend_labels.items()]
ax1.legend(handles=mk_patches, bbox_to_anchor=(0.5, -0.05), loc='upper center', ncol=3)
# 添加M-K趋势注释，放在显示区域下面，更靠近色标
plt.figtext(0.25, 0.49, "M-K增加和减少趋势判断标准是Z值正负", ha="center",
bbox={"facecolor": "white", "alpha": 0.8, "pad": 5})
# 2. Slope趋势栅格
ax2 = plt.subplot(222)
# 修改颜色映射，无趋势为黄色，但Slope没有斜率为0的像元，所以图例中不显示
cmap_slope = LinearSegmentedColormap.from_list(
'slope_cmap',
[(0, 'blue'), (0.5, 'yellow'), (1, 'red')]
)
im2 = ax2.imshow(slope_data, cmap=cmap_slope, vmin=-1, vmax=1)
ax2.set_title('Slope趋势栅格')
ax2.axis('off') # 去掉坐标轴
# 创建Slope趋势图例 - 调整位置更靠近显示区域，不包含无趋势
slope_legend_labels = {
'red': 'Slope增加',
'blue': 'Slope减少'
}
slope_patches = [mpatches.Patch(color=color, label=label) for color, label in slope_legend_labels.items()]
ax2.legend(handles=slope_patches, bbox_to_anchor=(0.5, -0.05), loc='upper center', ncol=2)
# 添加Slope趋势注释，放在显示区域下面，更靠近色标
plt.figtext(0.75, 0.49, "Slope增加和减少趋势判断标准是斜率值正负", ha="center",
bbox={"facecolor": "white", "alpha": 0.8, "pad": 5})
# 3. 比较结果栅格 - 放到第二行，使用223和224合并
ax3 = plt.subplot2grid((2, 2), (1, 0), colspan=2)
# 修改颜色映射，正趋势重合改为蓝绿色，负趋势重合改为水色
cmap_comparison = LinearSegmentedColormap.from_list(
'comparison_cmap',
[(0 / 3, 'magenta'), (1 / 3, '#2CA02C'), (2 / 3, '#1F77B4'), (3 / 3, '#FDE725')]
)
# 设置无效值(-9999)的颜色为透明
cmap_comparison.set_under('k', alpha=0)
im3 = ax3.imshow(comparison, cmap=cmap_comparison, vmin=0, vmax=3, clim=(0, 3))
ax3.set_title('趋势比较结果')
ax3.axis('off') # 去掉坐标轴
# 创建图例 - 调整位置更靠近显示区域，包含占比信息，保留三位小数
legend_labels = {
1: f'正趋势重合 ({percent_pos:.3f}%)',
2: f'负趋势重合 ({percent_neg:.3f}%)',
3: f'零趋势重合 ({percent_zero:.3f}%)',
0: f'不重合 ({percent_non_match:.3f}%)'
}
patches = [mpatches.Patch(color=cmap_comparison(i / 3), label=legend_labels[i]) for i in range(4)]
ax3.legend(handles=patches, bbox_to_anchor=(0.5, -0.05), loc='upper center', ncol=4)
# 添加注释
plt.figtext(0.5, 0.01, "注：不包含之前的24年均为0值的像元", ha="center",
bbox={"facecolor": "white", "alpha": 0.8, "pad": 5})
plt.tight_layout(rect=[0, 0.03, 1, 0.95]) # 调整布局，为底部的注释留出空间
# 保存主图表
plt.savefig(chart_path, dpi=300, bbox_inches='tight')
plt.close()
# 创建图表2：M-K像元值为5和Slope像元值为5-9的重合比较图
overlap_chart_path = os.path.join(output_dir, "M-K像元值为5和Slope像元值为5-9的重合比较图.png")
fig = plt.figure(figsize=(12, 10))
ax = plt.subplot(111)
# 创建自定义颜色映射，重合区域为绿色，背景为透明
cmap_overlap = LinearSegmentedColormap.from_list(
'overlap_cmap',
[(0, 'none'), (1, '#2CA02C')] # 0为透明，1为绿色
)
cmap_overlap.set_under('k', alpha=0) # 设置无效值为透明
im = ax.imshow(overlap_comparison, cmap=cmap_overlap, vmin=0, vmax=1, clim=(0, 1))
ax.set_title('M-K像元值为5和Slope像元值为5-9的重合区域')
ax.axis('off') # 去掉坐标轴
# 创建图例，更新为显示重叠区域占Slope=5-9区域的百分比
legend_labels = {
'#2CA02C': f'重合区域 ({overlap_percentage:.3f}%)',
'gray': f'重合区域占Slope=5-9区域的: {overlap_percent_of_slope:.3f}%'
}
patches = [mpatches.Patch(color=color, label=label) for color, label in legend_labels.items()]
ax.legend(handles=patches, bbox_to_anchor=(0.5, -0.05), loc='upper center', ncol=2)
# 添加注释
plt.figtext(0.5, 0.01, "注：M-K值=5且Slope值在5-9之间的像元", ha="center",
bbox={"facecolor": "white", "alpha": 0.8, "pad": 5})
plt.tight_layout(rect=[0, 0.03, 1, 0.95]) # 调整布局，为底部的注释留出空间
# 保存重合比较图
plt.savefig(overlap_chart_path, dpi=300, bbox_inches='tight')
plt.close()
# 打印详细统计信息
print(f"\n统计信息:")
print(f"1. M-K值为5的像元数: {mk_5_area}")
print(f"2. Slope值为5-9的像元数: {slope_5_9_area}")
print(f"3. 重合区域像元数: {overlap_area}")
print(f"4. 重合区域占Slope值5-9区域的百分比: {overlap_percent_of_slope:.3f}%")
print(f"5. 重合区域占总有效区域的百分比: {overlap_percentage:.3f}%")
print(f"\n图表已保存至: {chart_path}")
print(f"重合比较图已保存至: {overlap_chart_path}")
return chart_path, overlap_chart_path
except Exception as e:
print(f"比较栅格时出错: {e}")
return None, None
def main():
"""主函数，用于演示模块功能"""
# 设置实际数据路径
mk_raster_path = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\results\bijiaoM-K修正1.tif"
slope_raster_path = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\slope趋势分析\classified_slopeXX6.tif"
# 更新为您提供的shp文件路径
shp_path = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\黄河中游无定河流域\黄河中游无定河流域.shp"
output_dir = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\M-K空间分析\比较分析结果"
# 比较两个栅格
chart_path, overlap_chart_path = compare_rasters(mk_raster_path, slope_raster_path, shp_path, output_dir)
if chart_path and overlap_chart_path:
print(f"比较图表已生成: {chart_path}")
print(f"重合比较图已生成: {overlap_chart_path}")
if __name__ == "__main__":
main() 6、研究区水系图
需要流域水系shapfile文件一系列（后缀.shp.shx.dbf等）【可由研究区DEM文件在ArcGIS中提取得到】，以及研究区边界shapfile文件一系列（后缀.shp.shx.dbf等）。
这些shapfile文件的地理坐标系需一致，如不一样就不能如本代码这样直接叠加，修改方法可问AI。
import geopandas as gpd
import matplotlib.pyplot as plt
import os
import numpy as np
from matplotlib.ticker import MultipleLocator
# 设置中文字体
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]
# 文件路径设置
shp1_path = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\研究区对照描述水系图\quhoushiliang.shp"
shp2_path = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\研究区对照描述水系图\无定河流域.shp"
output_dir = r"D:\毕设文件\分析方法：小波分析和一列数据MK检验\研究区对照描述水系图\对照水系图输出文件"
output_path = os.path.join(output_dir, "叠加水系图1.png")
# 确保输出目录存在
os.makedirs(output_dir, exist_ok=True)
try:
# 读取第一个Shapefile (水系)
print("正在读取第一个Shapefile (水系)...")
shp1 = gpd.read_file(shp1_path)
print("第一个图层CRS:", shp1.crs)
# 读取第二个Shapefile (流域)
print("正在读取第二个Shapefile (流域)...")
shp2 = gpd.read_file(shp2_path)
print("第二个图层CRS:", shp2.crs)
# 强制设置为WGS 84 (Plate Carrée投影)
if shp1.crs is None:
shp1 = shp1.set_crs(epsg=4326)
if shp2.crs is None:
shp2 = shp2.set_crs(epsg=4326)
# 转换为Plate Carrée投影 (EPSG:4326)
shp1 = shp1.to_crs(epsg=4326)
shp2 = shp2.to_crs(epsg=4326)
print("转换后两个图层的CRS:", shp1.crs)
# 计算共同边界框
xmin1, ymin1, xmax1, ymax1 = shp1.total_bounds
xmin2, ymin2, xmax2, ymax2 = shp2.total_bounds
xmin = min(xmin1, xmin2)
ymin = min(ymin1, ymin2)
xmax = max(xmax1, xmax2)
ymax = max(ymax1, ymax2)
# 计算中心点经纬度
center_lat = (ymin + ymax) / 2
# 计算经纬度范围
lon_range = xmax - xmin
lat_range = ymax - ymin
# 计算正确的纵横比（考虑地球形状）
aspect_ratio = lon_range / lat_range * np.cos(np.radians(center_lat))
# 设置图形尺寸（以英寸为单位）
fig_width = 12 # 宽度固定
fig_height = fig_width / aspect_ratio # 根据比例计算高度
# 创建图形
print("正在创建叠加图...")
fig, ax = plt.subplots(figsize=(fig_width, fig_height))
# 叠加第一个Shapefile (水系)
shp1.plot(ax=ax, color='darkblue', linewidth=1.5, label='水系')
# 叠加第二个Shapefile (流域边界)
shp2.boundary.plot(ax=ax, color='darkred', linewidth=2, linestyle='--', label='流域边界')
shp2.plot(ax=ax, color='none', edgecolor='darkred', alpha=0.1)
# 设置显示范围（添加适当边距）
margin_x = lon_range * 0.1
margin_y = lat_range * 0.1
ax.set_xlim(xmin - margin_x, xmax + margin_x)
ax.set_ylim(ymin - margin_y, ymax + margin_y)
# 设置坐标轴比例为1:1，确保网格为正方形
ax.set_aspect('equal')
# 添加经纬度网格线
gl = ax.grid(True, linestyle='--', alpha=0.7, color='gray')
ax.set_axisbelow(True) # 将网格线置于底层
# 设置网格间隔（根据数据范围自动调整）
if lon_range > 10:
lon_interval = 1
elif lon_range > 2:
lon_interval = 0.5
else:
lon_interval = 0.1
if lat_range > 10:
lat_interval = 1
elif lat_range > 2:
lat_interval = 0.5
else:
lat_interval = 0.1
# 设置坐标轴刻度
ax.xaxis.set_major_locator(MultipleLocator(lon_interval))
ax.yaxis.set_major_locator(MultipleLocator(lat_interval))
# 添加标题和图例（修改标题为"研究区水系图"）
plt.title('研究区水系图', fontsize=16)
plt.legend(loc='upper right', fontsize=12)
plt.xlabel('经度', fontsize=12)
plt.ylabel('纬度', fontsize=12)
# 保存图形（使用相同的DPI确保比例正确）
print(f"正在保存图形到: {output_path}")
plt.tight_layout()
plt.savefig(output_path, dpi=300, bbox_inches='tight')
print("叠加图创建完成!")
plt.show()
except Exception as e:
print(f"发生错误: {e}")
DAMO开发者矩阵，由阿里巴巴达摩院和中国互联网协会联合发起，致力于探讨最前沿的技术趋势与应用成果，搭建高质量的交流与分享平台，推动技术创新与产业应用链接，围绕“人工智能与新型计算”构建开放共享的开发者生态。
更多推荐
1
0- 0
已为社区贡献1条内容
所有评论(0)