import os
import sys
import numpy as np
import geopandas as gpd
import rasterio
from rasterio.mask import mask
from rasterio.warp import calculate_default_transform, reproject, Resampling
import matplotlib.pyplot as plt
import matplotlib
import pandas as pd
from shapely.geometry import box, LineString, Point
import warnings
import glob
from datetime import datetime
import json
import tempfile
import shutil
warnings.filterwarnings('ignore')
# ==================== 配置部分 ====================
# 设置中文字体
matplotlib.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans']
matplotlib.rcParams['axes.unicode_minus'] = False
# 基础路径 - 请根据实际情况修改
BASE_DIR = r"C:\indian data"
# 数据目录
MODIS_DIR = os.path.join(BASE_DIR, "indian modis")
NIGHTLIGHT_DIR = os.path.join(BASE_DIR, "indian night")
OSM_DIR = os.path.join(BASE_DIR, "indian road")
BOUNDARY_DIR = os.path.join(BASE_DIR, "indian boundary")
# 输出目录
OUTPUT_DIR = r"C:\indian data\out"
os.makedirs(OUTPUT_DIR, exist_ok=True)
# 临时目录
OUTPUT_DIR = r"C:\indian data\linshi"
os.makedirs(OUTPUT_DIR, exist_ok=True)
# 分析年份
MODIS_YEARS = list(range(2015, 2025)) # 2015-2024
NIGHTLIGHT_YEARS = [2020, 2024]
OSM_YEARS = [2020, 2025]
# 目标投影(印度常用UTM 44N)
TARGET_CRS = "EPSG:32644"
# MODIS分类方案(多种可能的方案)
MODIS_CLASS_SCHEMES = {
'IGBP': {
0: "水体", 1: "常绿针叶林", 2: "常绿阔叶林", 3: "落叶针叶林",
4: "落叶阔叶林", 5: "混交林", 6: "密闭灌丛", 7: "开阔灌丛",
8: "木本稀树草原", 9: "稀树草原", 10: "草地", 11: "永久湿地",
12: "耕地", 13: "城市和建设用地", 14: "耕地/自然植被镶嵌",
15: "冰雪", 16: "贫瘠地"
}
}
# ==================== 工具函数 ====================
def copy_to_temp(src_path):
"""复制文件到临时英文路径,解决中文路径问题"""
if not os.path.exists(src_path):
print(f"文件不存在: {src_path}")
return None
try:
# 获取文件名并创建安全的英文名
filename = os.path.basename(src_path)
# 替换中文字符为下划线
safe_name = "".join([c if ord(c) < 128 else "_" for c in filename])
temp_path = os.path.join(TEMP_DIR, safe_name)
shutil.copy2(src_path, temp_path)
print(f"已复制文件到临时路径: {safe_name}")
return temp_path
except Exception as e:
print(f"复制文件失败: {e}")
return None
def find_files_by_year(directory, year, extensions):
"""根据年份查找文件"""
found_files = []
for ext in extensions:
patterns = [
f"*{year}*{ext}",
f"*{str(year)[2:]}*{ext}" # 最后两位数字
]
for pattern in patterns:
files = glob.glob(os.path.join(directory, pattern))
found_files.extend(files)
# 如果没有找到,尝试搜索所有文件
if not found_files:
for ext in extensions:
all_files = glob.glob(os.path.join(directory, f"*{ext}"))
# 检查文件名中是否包含年份
for file in all_files:
if str(year) in os.path.basename(file):
found_files.append(file)
return found_files
def detect_modis_class_scheme(filepath):
"""检测MODIS数据的分类方案"""
try:
with rasterio.open(filepath) as src:
data = src.read(1)
unique_values = np.unique(data)
print(f"数据唯一值: {unique_values[:20]}") # 显示前20个值
# 根据值的范围判断分类方案
max_val = np.max(unique_values)
if max_val <= 16:
urban_class = 13 # IGBP中的城市用地
scheme = 'IGBP'
elif max_val <= 14:
urban_class = 12 # UMD中的城市用地
scheme = 'UMD'
else:
# 尝试找到最可能的城市类别
urban_candidates = [v for v in unique_values if 10 <= v <= 20]
if urban_candidates:
urban_class = urban_candidates[0]
scheme = 'Custom'
else:
urban_class = 13 # 默认
scheme = 'Unknown'
print(f"检测到分类方案: {scheme}, 城市用地ID: {urban_class}")
return urban_class, scheme
except Exception as e:
print(f"检测分类方案失败: {e}")
return None, None
# ==================== 数据加载函数 ====================
def load_boundary_data():
"""加载印度边界数据"""
print("加载印度边界数据...")
# 查找边界文件
boundary_files = []
for root, dirs, files in os.walk(BOUNDARY_DIR):
for file in files:
if file.endswith(('.shp', '.geojson', '.gpkg')):
boundary_files.append(os.path.join(root, file))
if not boundary_files:
print("错误: 未找到边界文件")
return None
try:
boundary_path = boundary_files[0]
print(f"找到边界文件: {os.path.basename(boundary_path)}")
# 复制到临时目录避免中文路径问题
temp_boundary = copy_to_temp(boundary_path)
if temp_boundary is None:
print("错误: 无法复制边界文件")
return None
boundary = gpd.read_file(temp_boundary)
if boundary.crs is None:
boundary.set_crs("EPSG:4326", inplace=True)
print(f"边界数据加载完成,CRS: {boundary.crs}, 要素数量: {len(boundary)}")
# 统一到目标投影
if boundary.crs != TARGET_CRS:
boundary = boundary.to_crs(TARGET_CRS)
# 保存处理后的边界
output_path = os.path.join(OUTPUT_DIR, "india_boundary_processed.gpkg")
boundary.to_file(output_path, driver="GPKG")
print(f"边界数据已保存: {output_path}")
return boundary
except Exception as e:
print(f"加载边界文件失败: {e}")
return None
def find_and_load_modis_data(years, boundary):
"""查找并加载MODIS数据"""
print("查找并加载MODIS数据...")
modis_files = {}
urban_stats = {}
class_distributions = {}
for year in years:
print(f" 查找 {year} 年数据...")
# 查找文件
extensions = ['.tif', '.tiff', '.hdf', '.h5']
found_files = find_files_by_year(MODIS_DIR, year, extensions)
if not found_files:
print(f" 警告: 未找到{year}年的MODIS数据")
continue
# 使用第一个找到的文件
filepath = found_files[0]
print(f" 找到文件: {os.path.basename(filepath)}")
# 复制到临时目录
temp_file = copy_to_temp(filepath)
if temp_file is None:
print(f" 错误: 无法复制文件")
continue
try:
# 检测分类方案
urban_class, scheme = detect_modis_class_scheme(temp_file)
if urban_class is None:
print(f" 错误: 无法检测分类方案")
continue
# 处理数据
processed_path = process_raster_data(temp_file, boundary, f"modis_{year}")
if processed_path is None:
print(f" 错误: 无法处理数据")
continue
modis_files[year] = processed_path
# 分析城市用地
stats = analyze_modis_urban_area(processed_path, urban_class)
if stats:
urban_stats[year] = stats
print(f" {year}年城市面积: {stats['urban_area_km2']:,.0f} km² ({stats['urban_percentage']:.2f}%)")
# 记录分类分布
dist = get_class_distribution(processed_path)
class_distributions[year] = dist
except Exception as e:
print(f" 处理{year}年数据失败: {e}")
return modis_files, urban_stats, class_distributions
def find_and_load_nightlight_data(years, boundary):
"""查找并加载夜间灯光数据"""
print("查找并加载夜间灯光数据...")
nightlight_files = {}
nightlight_stats = {}
for year in years:
print(f" 查找 {year} 年数据...")
# 查找TIFF文件
extensions = ['.tif', '.tiff']
found_files = find_files_by_year(NIGHTLIGHT_DIR, year, extensions)
if not found_files:
print(f" 警告: 未找到{year}年的夜间灯光数据")
continue
# 使用第一个找到的文件
filepath = found_files[0]
print(f" 找到文件: {os.path.basename(filepath)}")
# 复制到临时目录
temp_file = copy_to_temp(filepath)
if temp_file is None:
print(f" 错误: 无法复制文件")
continue
try:
# 处理数据
processed_path = process_raster_data(temp_file, boundary, f"nightlight_{year}")
if processed_path is None:
print(f" 错误: 无法处理数据")
continue
nightlight_files[year] = processed_path
# 分析夜间灯光
stats = analyze_nightlight_data(processed_path)
if stats:
nightlight_stats[year] = stats
print(f" {year}年夜光均值: {stats['mean_value']:.2f}, 灯光面积: {stats['light_percentage']:.2f}%")
except Exception as e:
print(f" 处理{year}年数据失败: {e}")
return nightlight_files, nightlight_stats
def find_and_load_osm_data(years, boundary):
"""查找并加载OSM道路数据"""
print("查找并加载OSM道路数据...")
roads_data = {}
for year in years:
print(f" 查找 {year} 年数据...")
# 查找OSM文件
extensions = ['.pbf', '.osm', '.shp', '.geojson', '.gpkg']
found_files = []
for ext in extensions:
pattern = f"*{year}*{ext}"
files = glob.glob(os.path.join(OSM_DIR, pattern))
found_files.extend(files)
if not found_files:
# 如果没有年份标记的文件,尝试查找所有OSM文件
for ext in extensions:
files = glob.glob(os.path.join(OSM_DIR, f"*{ext}"))
found_files.extend(files)
if not found_files:
print(f" 警告: 未找到{year}年的OSM数据")
continue
# 使用第一个找到的文件
filepath = found_files[0]
print(f" 找到文件: {os.path.basename(filepath)}")
try:
# 加载道路数据
roads = load_osm_file(filepath, boundary)
if roads is not None and not roads.empty:
roads_data[year] = roads
print(f" {year}年道路数据: {len(roads)} 条道路")
else:
print(f" 错误: 无法加载道路数据")
except Exception as e:
print(f" 处理{year}年数据失败: {e}")
return roads_data
def load_osm_file(filepath, boundary):
"""加载OSM文件"""
try:
# 复制到临时目录
temp_file = copy_to_temp(filepath)
if temp_file is None:
return None
# 根据文件扩展名选择加载方法
ext = os.path.splitext(temp_file)[1].lower()
if ext in ['.shp', '.geojson', '.gpkg']:
# 使用geopandas直接读取
roads = gpd.read_file(temp_file)
elif ext in ['.pbf', '.osm']:
# 尝试使用osmnx读取
roads = load_osm_pbf(temp_file, boundary)
if roads is None:
return None
else:
print(f" 不支持的文件格式: {ext}")
return None
# 筛选主要道路
if 'highway' in roads.columns:
highway_types = ['motorway', 'trunk', 'primary', 'secondary', 'tertiary']
roads = roads[roads['highway'].isin(highway_types)]
print(f" 筛选后保留 {len(roads)} 条主要道路")
# 统一投影
if roads.crs != boundary.crs:
roads = roads.to_crs(boundary.crs)
# 计算道路长度
roads['length_km'] = roads.geometry.length / 1000
return roads
except Exception as e:
print(f" 加载OSM文件失败: {e}")
return None
def load_osm_pbf(filepath, boundary):
"""加载PBF格式的OSM数据"""
try:
# 尝试使用osmnx
import osmnx as ox
print(f" 使用osmnx加载PBF文件...")
# 获取边界范围
boundary_wgs84 = boundary.to_crs("EPSG:4326")
bounds = boundary_wgs84.total_bounds
# 使用边界下载数据
tags = {'highway': True}
roads = ox.features_from_bbox(
north=bounds[3], south=bounds[1],
east=bounds[2], west=bounds[0],
tags=tags
)
return roads
except ImportError:
print(f" osmnx未安装,无法读取PBF文件")
return None
except Exception as e:
print(f" 使用osmnx加载失败: {e}")
return None
# ==================== 数据处理函数 ====================
def process_raster_data(input_path, boundary, prefix):
"""处理栅格数据(重投影和裁剪)"""
try:
with rasterio.open(input_path) as src:
src_crs = src.crs
# 如果投影不同,需要重投影
if src_crs != boundary.crs:
output_path = os.path.join(TEMP_DIR, f"{prefix}_reprojected.tif")
# 计算目标投影参数
transform, width, height = calculate_default_transform(
src.crs, boundary.crs,
src.width, src.height,
*src.bounds
)
# 设置输出参数
kwargs = src.meta.copy()
kwargs.update({
'crs': boundary.crs,
'transform': transform,
'width': width,
'height': height
})
# 重投影
with rasterio.open(output_path, 'w', **kwargs) as dst:
for i in range(1, src.count + 1):
reproject(
source=rasterio.band(src, i),
destination=rasterio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=boundary.crs,
resampling=Resampling.bilinear
)
reprojected_path = output_path
else:
reprojected_path = input_path
# 裁剪到边界
clipped_path = os.path.join(TEMP_DIR, f"{prefix}_clipped.tif")
with rasterio.open(reprojected_path) as src2:
boundary_proj = boundary.to_crs(src2.crs)
out_image, out_transform = mask(
src2,
boundary_proj.geometry,
crop=True,
all_touched=True
)
out_meta = src2.meta.copy()
out_meta.update({
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform
})
with rasterio.open(clipped_path, "w", **out_meta) as dst:
dst.write(out_image)
return clipped_path
except Exception as e:
print(f" 处理栅格数据失败: {e}")
return None
def analyze_modis_urban_area(filepath, urban_class):
"""分析MODIS城市用地"""
try:
with rasterio.open(filepath) as src:
data = src.read(1)
# 统计城市用地
urban_mask = data == urban_class
urban_pixels = np.sum(urban_mask)
total_pixels = data.size
# 计算面积(假设500m分辨率)
pixel_area_km2 = 0.25 # 0.5km * 0.5km = 0.25 km²
urban_area_km2 = urban_pixels * pixel_area_km2
urban_percentage = (urban_pixels / total_pixels) * 100
# 统计其他可能与交通相关的类别
transport_related_classes = [12, 13, 16] # 耕地、城市、贫瘠地
transport_pixels = 0
for class_id in transport_related_classes:
if class_id != urban_class:
transport_pixels += np.sum(data == class_id)
transport_area_km2 = transport_pixels * pixel_area_km2
transport_percentage = (transport_pixels / total_pixels) * 100
return {
'urban_area_km2': urban_area_km2,
'urban_percentage': urban_percentage,
'urban_pixels': int(urban_pixels),
'transport_area_km2': transport_area_km2,
'transport_percentage': transport_percentage,
'total_pixels': int(total_pixels),
'urban_class_id': urban_class
}
except Exception as e:
print(f" 分析MODIS城市用地失败: {e}")
return None
def analyze_nightlight_data(filepath):
"""分析夜间灯光数据"""
try:
with rasterio.open(filepath) as src:
data = src.read(1)
# 处理nodata值
if src.nodata is not None:
data = np.where(data == src.nodata, 0, data)
# 确保数据非负
data = np.where(data < 0, 0, data)
# 计算统计信息
valid_data = data[data > 0]
if len(valid_data) > 0:
mean_value = np.mean(valid_data)
median_value = np.median(valid_data)
max_value = np.max(valid_data)
total_value = np.sum(valid_data)
# 计算有灯光的区域
if len(valid_data) > 100:
threshold = np.percentile(valid_data, 50) # 中位数作为阈值
else:
threshold = mean_value * 0.5
light_pixels = np.sum(data > threshold)
total_pixels = data.size
light_percentage = (light_pixels / total_pixels) * 100
else:
mean_value = 0
light_percentage = 0
total_value = 0
light_pixels = 0
return {
'mean_value': mean_value,
'total_value': total_value,
'light_percentage': light_percentage,
'light_pixels': int(light_pixels),
'total_pixels': int(total_pixels)
}
except Exception as e:
print(f" 分析夜间灯光数据失败: {e}")
return None
def get_class_distribution(filepath):
"""获取分类分布"""
try:
with rasterio.open(filepath) as src:
data = src.read(1)
unique_vals, counts = np.unique(data, return_counts=True)
return dict(zip(unique_vals, counts))
except Exception as e:
print(f" 获取分类分布失败: {e}")
return {}
# ==================== 分析函数 ====================
def analyze_urban_expansion(urban_stats):
"""分析城市扩张趋势"""
print("\n分析城市扩张趋势...")
if len(urban_stats) < 2:
print(" 数据不足,需要至少两年数据")
return None
years = sorted(urban_stats.keys())
# 创建DataFrame以便分析
df_data = []
for year in years:
stats = urban_stats[year]
df_data.append({
'year': year,
'urban_area_km2': stats['urban_area_km2'],
'urban_percentage': stats['urban_percentage'],
'transport_area_km2': stats['transport_area_km2'],
'transport_percentage': stats['transport_percentage']
})
df = pd.DataFrame(df_data)
print(" 年度城市扩张统计:")
print(df.to_string(index=False))
# 计算变化
changes = []
for i in range(1, len(years)):
year1, year2 = years[i-1], years[i]
area1 = urban_stats[year1]['urban_area_km2']
area2 = urban_stats[year2]['urban_area_km2']
# 避免除零
if area1 > 0:
change_km2 = area2 - area1
change_percent = ((area2 - area1) / area1) * 100
annual_rate = change_percent / (year2 - year1)
else:
change_km2 = area2
change_percent = 100 if area2 > 0 else 0
annual_rate = change_percent / (year2 - year1)
changes.append({
'period': f"{year1}-{year2}",
'area_change_km2': change_km2,
'percent_change': change_percent,
'annual_rate': annual_rate
})
print(f" {year1}-{year2}: 面积变化 {change_km2:+,.0f} km², "
f"变化率 {change_percent:+.2f}%, 年化率 {annual_rate:+.2f}%/年")
# 总体趋势
if len(years) >= 2:
first_year = years[0]
last_year = years[-1]
first_area = urban_stats[first_year]['urban_area_km2']
last_area = urban_stats[last_year]['urban_area_km2']
if first_area > 0:
total_change = last_area - first_area
total_percent_change = (total_change / first_area) * 100
avg_annual_rate = total_percent_change / (last_year - first_year)
else:
total_change = last_area
total_percent_change = 100 if last_area > 0 else 0
avg_annual_rate = total_percent_change / (last_year - first_year)
print(f"\n 总体变化 ({first_year}-{last_year}):")
print(f" 起始面积: {first_area:,.0f} km²")
print(f" 结束面积: {last_area:,.0f} km²")
print(f" 总面积变化: {total_change:+,.0f} km²")
print(f" 总变化率: {total_percent_change:+.2f}%")
print(f" 年均变化率: {avg_annual_rate:+.2f}%/年")
return {
'years': years,
'df': df,
'changes': changes,
'urban_stats': urban_stats
}
def analyze_road_infrastructure(roads_data, boundary):
"""分析道路基础设施"""
print("\n分析道路基础设施...")
if not roads_data:
print(" 无道路数据")
return None
all_stats = {}
for year, roads in roads_data.items():
if roads.empty:
print(f" {year}年: 无道路数据")
continue
# 基本统计
stats = {
'year': year,
'total_roads': len(roads),
'total_length_km': float(roads['length_km'].sum()),
'avg_length_km': float(roads['length_km'].mean()),
'max_length_km': float(roads['length_km'].max())
}
# 按道路类型统计
if 'highway' in roads.columns:
type_stats = roads.groupby('highway').agg({
'length_km': ['sum', 'count', 'mean']
}).round(1)
type_summary = {}
for htype, row in type_stats.iterrows():
type_summary[htype] = {
'count': int(row[('length_km', 'count')]),
'total_length': float(row[('length_km', 'sum')])
}
stats['type_summary'] = type_summary
all_stats[year] = stats
print(f" {year}年:")
print(f" 道路总数: {stats['total_roads']}")
print(f" 总长度: {stats['total_length_km']:.0f} km")
print(f" 平均长度: {stats['avg_length_km']:.1f} km")
if 'type_summary' in stats:
for htype, summary in stats['type_summary'].items():
print(f" {htype}: {summary['count']}条, 总长度 {summary['total_length']:.0f} km")
return all_stats
# ==================== 可视化函数 ====================
def create_urban_expansion_plot(urban_analysis, output_path):
"""创建城市扩张可视化图表"""
print("创建城市扩张可视化图表...")
if urban_analysis is None:
print(" 无城市扩张数据可可视化")
return None
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
years = urban_analysis['years']
urban_stats = urban_analysis['urban_stats']
urban_areas = [urban_stats[y]['urban_area_km2'] for y in years]
urban_percentages = [urban_stats[y]['urban_percentage'] for y in years]
transport_areas = [urban_stats[y]['transport_area_km2'] for y in years]
# 1. 城市面积趋势
ax = axes[0, 0]
bars = ax.bar(years, urban_areas, color='#FF6B6B', alpha=0.7)
ax.set_xlabel('年份', fontsize=12)
ax.set_ylabel('城市面积 (km²)', fontsize=12)
ax.set_title('印度城市面积变化趋势 (2015-2024)', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
# 添加数值标签
for bar, area in zip(bars, urban_areas):
height = bar.get_height()
if height > 0:
ax.text(bar.get_x() + bar.get_width()/2., height + max(urban_areas)*0.01,
f'{area:,.0f}', ha='center', va='bottom', fontsize=9)
# 2. 城市面积占比
ax = axes[0, 1]
ax.plot(years, urban_percentages, 'o-', linewidth=2, markersize=8, color='#4ECDC4')
ax.set_xlabel('年份', fontsize=12)
ax.set_ylabel('城市面积占比 (%)', fontsize=12)
ax.set_title('城市面积占比变化', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
# 添加数值标签
for x, y in zip(years, urban_percentages):
if y > 0:
ax.text(x, y + max(urban_percentages)*0.02, f'{y:.2f}%',
ha='center', va='bottom', fontsize=9)
# 3. 城市与交通用地对比
ax = axes[1, 0]
x = np.arange(len(years))
width = 0.35
bars1 = ax.bar(x - width/2, urban_areas, width, label='城市用地', color='#FF6B6B', alpha=0.7)
bars2 = ax.bar(x + width/2, transport_areas, width, label='交通相关用地', color='#45B7D1', alpha=0.7)
ax.set_xlabel('年份', fontsize=12)
ax.set_ylabel('面积 (km²)', fontsize=12)
ax.set_title('城市用地与交通相关用地对比', fontsize=14, fontweight='bold')
ax.set_xticks(x)
ax.set_xticklabels(years)
ax.legend()
ax.grid(True, alpha=0.3)
# 4. 增长率
ax = axes[1, 1]
if urban_analysis['changes']:
changes = urban_analysis['changes']
periods = [c['period'] for c in changes]
percent_changes = [c['percent_change'] for c in changes]
bars = ax.bar(periods, percent_changes,
color=['#96CEB4' if c > 0 else '#FF6B6B' for c in percent_changes],
alpha=0.7)
ax.set_xlabel('时期', fontsize=12)
ax.set_ylabel('变化率 (%)', fontsize=12)
ax.set_title('城市面积年增长率', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
ax.axhline(y=0, color='black', linestyle='-', linewidth=0.5)
# 添加数值标签
for bar, change in zip(bars, percent_changes):
height = bar.get_height()
if abs(height) > 0.1: # 忽略太小的变化
ax.text(bar.get_x() + bar.get_width()/2.,
height + (1 if height >= 0 else -1),
f'{change:+.1f}%',
ha='center', va='bottom' if height >= 0 else 'top', fontsize=9)
plt.suptitle('印度城市扩张与土地利用变化分析', fontsize=16, fontweight='bold', y=0.98)
plt.tight_layout()
plt.savefig(output_path, dpi=300, bbox_inches='tight')
plt.show()
print(f"图表已保存: {output_path}")
return output_path
def create_road_analysis_plot(road_stats, boundary, roads_data, output_path):
"""创建道路分析图表"""
print("创建道路分析图表...")
if not road_stats:
print(" 无道路数据可可视化")
return None
fig, axes = plt.subplots(1, 2, figsize=(16, 8))
# 1. 道路网络地图
ax = axes[0]
boundary.plot(ax=ax, color='lightgray', edgecolor='black', alpha=0.3)
# 绘制所有年份的道路
colors = ['red', 'blue', 'green', 'orange', 'purple']
for i, (year, roads) in enumerate(roads_data.items()):
if not roads.empty:
color = colors[i % len(colors)]
roads.plot(ax=ax, color=color, linewidth=1, alpha=0.7, label=f'{year}年')
ax.set_title('印度主要道路网络', fontsize=14, fontweight='bold')
ax.axis('off')
ax.legend(loc='upper right', fontsize=10)
# 2. 道路统计图表
ax = axes[1]
years = sorted(road_stats.keys())
total_lengths = [road_stats[y]['total_length_km'] for y in years]
total_roads = [road_stats[y]['total_roads'] for y in years]
x = np.arange(len(years))
width = 0.35
bars1 = ax.bar(x - width/2, total_lengths, width, label='总长度 (km)', color='#4CAF50', alpha=0.7)
ax.set_xlabel('年份', fontsize=12)
ax.set_ylabel('总长度 (km)', fontsize=12, color='#4CAF50')
ax.tick_params(axis='y', labelcolor='#4CAF50')
ax2 = ax.twinx()
bars2 = ax2.bar(x + width/2, total_roads, width, label='道路数量', color='#2196F3', alpha=0.7)
ax2.set_ylabel('道路数量', fontsize=12, color='#2196F3')
ax2.tick_params(axis='y', labelcolor='#2196F3')
ax.set_title('道路基础设施统计', fontsize=14, fontweight='bold')
ax.set_xticks(x)
ax.set_xticklabels(years)
# 合并图例
lines1, labels1 = ax.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax.legend(lines1 + lines2, labels1 + labels2, loc='upper left')
plt.suptitle('印度道路基础设施分析', fontsize=16, fontweight='bold', y=0.98)
plt.tight_layout()
plt.savefig(output_path, dpi=300, bbox_inches='tight')
plt.show()
print(f"图表已保存: {output_path}")
return output_path
def create_nightlight_analysis_plot(nightlight_stats, output_path):
"""创建夜间灯光分析图表"""
print("创建夜间灯光分析图表...")
if not nightlight_stats or len(nightlight_stats) < 2:
print(" 无足够的夜间灯光数据可可视化")
return None
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
years = sorted(nightlight_stats.keys())
light_percentages = [nightlight_stats[y]['light_percentage'] for y in years]
mean_values = [nightlight_stats[y]['mean_value'] for y in years]
light_pixels = [nightlight_stats[y]['light_pixels'] for y in years]
# 1. 灯光面积占比
ax = axes[0]
bars = ax.bar(years, light_percentages, color='#FFD166', alpha=0.7)
ax.set_xlabel('年份', fontsize=12)
ax.set_ylabel('灯光面积占比 (%)', fontsize=12)
ax.set_title('夜间灯光覆盖面积变化', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
for bar, percent in zip(bars, light_percentages):
height = bar.get_height()
if height > 0:
ax.text(bar.get_x() + bar.get_width()/2., height + max(light_percentages)*0.02,
f'{percent:.2f}%', ha='center', va='bottom', fontsize=10)
# 2. 灯光均值变化
ax = axes[1]
ax.plot(years, mean_values, 'o-', linewidth=2, markersize=8, color='#06D6A0')
ax.set_xlabel('年份', fontsize=12)
ax.set_ylabel('平均亮度值', fontsize=12)
ax.set_title('夜间灯光平均亮度变化', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
for x, y in zip(years, mean_values):
if y > 0:
ax.text(x, y + max(mean_values)*0.02, f'{y:.2f}',
ha='center', va='bottom', fontsize=10)
# 3. 灯光像素数量
ax = axes[2]
bars = ax.bar(years, light_pixels, color='#EF476F', alpha=0.7)
ax.set_xlabel('年份', fontsize=12)
ax.set_ylabel('灯光像素数量', fontsize=12)
ax.set_title('夜间灯光像素数量变化', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
for bar, pixels in zip(bars, light_pixels):
height = bar.get_height()
if height > 0:
ax.text(bar.get_x() + bar.get_width()/2., height + max(light_pixels)*0.01,
f'{pixels:,}', ha='center', va='bottom', fontsize=9)
plt.suptitle('印度夜间灯光变化分析', fontsize=16, fontweight='bold', y=0.98)
plt.tight_layout()
plt.savefig(output_path, dpi=300, bbox_inches='tight')
plt.show()
print(f"图表已保存: {output_path}")
return output_path
# ==================== 报告生成 ====================
def generate_analysis_report(urban_analysis, road_stats, nightlight_stats,
class_distributions, output_path):
"""生成分析报告"""
print("生成分析报告...")
with open(output_path, 'w', encoding='utf-8') as f:
f.write("=" * 80 + "\n")
f.write("印度交通基础设施变化监测分析报告\n")
f.write("=" * 80 + "\n\n")
f.write(f"报告生成时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n")
f.write(f"分析时期: 2015-2024年\n")
f.write(f"数据来源: MODIS土地覆盖数据、VIIRS夜间灯光数据、OSM道路数据\n\n")
# 1. 数据概览
f.write("1. 数据概览\n")
f.write("-" * 60 + "\n")
if urban_analysis:
f.write(f" MODIS数据年份: {', '.join(map(str, sorted(urban_analysis['years'])))}\n")
else:
f.write(f" MODIS数据: 无有效数据\n")
if nightlight_stats:
f.write(f" 夜间灯光数据年份: {', '.join(map(str, sorted(nightlight_stats.keys())))}\n")
else:
f.write(f" 夜间灯光数据: 无有效数据\n")
if road_stats:
f.write(f" OSM道路数据年份: {', '.join(map(str, sorted(road_stats.keys())))}\n")
else:
f.write(f" OSM道路数据: 无有效数据\n")
f.write("\n")
# 2. 城市扩张分析
f.write("2. 城市扩张分析\n")
f.write("-" * 60 + "\n")
if urban_analysis:
years = sorted(urban_analysis['years'])
f.write(" 2.1 年度城市面积统计:\n")
f.write(" 年份 城市面积(km²) 占比(%) 交通相关用地(km²)\n")
f.write(" -----------------------------------------------------\n")
for year in years:
stats = urban_analysis['urban_stats'][year]
f.write(f" {year} {stats['urban_area_km2']:>10,.0f} {stats['urban_percentage']:>6.2f} {stats['transport_area_km2']:>15,.0f}\n")
f.write("\n")
if urban_analysis['changes']:
f.write(" 2.2 城市扩张趋势:\n")
f.write(" 时期 面积变化(km²) 变化率(%) 年化率(%/年)\n")
f.write(" -------------------------------------------------\n")
for change in urban_analysis['changes']:
f.write(f" {change['period']} {change['area_change_km2']:>11,.0f} {change['percent_change']:>8.2f} {change['annual_rate']:>10.2f}\n")
f.write("\n")
else:
f.write(" 无有效的城市扩张数据\n\n")
# 3. 道路基础设施分析
f.write("3. 道路基础设施分析\n")
f.write("-" * 60 + "\n")
if road_stats:
for year, stats in sorted(road_stats.items()):
f.write(f" {year}年:\n")
f.write(f" 道路总数: {stats['total_roads']} 条\n")
f.write(f" 总长度: {stats['total_length_km']:,.0f} km\n")
f.write(f" 平均长度: {stats['avg_length_km']:.1f} km\n")
if 'type_summary' in stats:
f.write(f" 道路类型分布:\n")
for htype, summary in stats['type_summary'].items():
f.write(f" {htype}: {summary['count']}条, 总长度 {summary['total_length']:.0f} km\n")
f.write("\n")
else:
f.write(" 无有效的道路数据\n\n")
# 4. 夜间灯光分析
f.write("4. 夜间灯光分析\n")
f.write("-" * 60 + "\n")
if nightlight_stats:
f.write(" 4.1 夜间灯光统计:\n")
f.write(" 年份 平均亮度 灯光面积(%) 灯光像素数\n")
f.write(" -----------------------------------------\n")
for year in sorted(nightlight_stats.keys()):
stats = nightlight_stats[year]
f.write(f" {year} {stats['mean_value']:>8.2f} {stats['light_percentage']:>10.2f} {stats['light_pixels']:>12,}\n")
f.write("\n")
if len(nightlight_stats) >= 2:
years = sorted(nightlight_stats.keys())
first_year = years[0]
last_year = years[-1]
change_percent = (nightlight_stats[last_year]['light_percentage'] -
nightlight_stats[first_year]['light_percentage'])
change_value = (nightlight_stats[last_year]['mean_value'] -
nightlight_stats[first_year]['mean_value'])
f.write(f" 4.2 变化趋势 ({first_year}-{last_year}):\n")
f.write(f" 灯光面积变化: {change_percent:+.2f}%\n")
f.write(f" 平均亮度变化: {change_value:+.2f}\n")
f.write("\n")
else:
f.write(" 无有效的夜间灯光数据\n\n")
# 5. 分类分布信息
f.write("5. MODIS分类分布\n")
f.write("-" * 60 + "\n")
if class_distributions:
sample_year = list(class_distributions.keys())[0]
if class_distributions[sample_year]:
f.write(f" {sample_year}年分类分布 (前10个类别):\n")
total_pixels = sum(class_distributions[sample_year].values())
sorted_classes = sorted(class_distributions[sample_year].items(),
key=lambda x: x[1], reverse=True)[:10]
for class_id, count in sorted_classes:
percentage = (count / total_pixels) * 100
class_name = MODIS_CLASS_SCHEMES.get('IGBP', {}).get(class_id, f"类别{class_id}")
f.write(f" 类别 {class_id:2d} ({class_name}): {count:>10,} 像素 ({percentage:6.2f}%)\n")
f.write("\n")
else:
f.write(" 无分类分布数据\n\n")
# 6. 结论与建议
f.write("6. 结论与建议\n")
f.write("-" * 60 + "\n")
f.write(" 6.1 主要发现:\n")
if urban_analysis and urban_analysis['changes']:
changes = urban_analysis['changes']
if changes:
last_change = changes[-1]
f.write(f" • 城市扩张趋势: {last_change['period']}期间城市面积变化 {last_change['area_change_km2']:+,.0f} km²\n")
if road_stats:
total_roads = sum(stats['total_roads'] for stats in road_stats.values())
total_length = sum(stats['total_length_km'] for stats in road_stats.values())
f.write(f" • 道路基础设施: 共分析 {total_roads} 条主要道路,总长度 {total_length:,.0f} km\n")
if nightlight_stats and len(nightlight_stats) >= 2:
years = sorted(nightlight_stats.keys())
change = (nightlight_stats[years[-1]]['light_percentage'] -
nightlight_stats[years[0]]['light_percentage'])
f.write(f" • 夜间灯光变化: 灯光面积增加 {change:+.2f}%\n")
f.write("\n 6.2 结论:\n")
f.write(" • 基于多源遥感数据分析,可以监测印度交通基础设施变化\n")
f.write(" • 城市扩张、道路网络发展和夜间灯光变化呈现相关性\n")
f.write(" • 需要进一步结合高分辨率数据验证分析结果\n")
f.write("\n 6.3 建议:\n")
f.write(" 1. 定期更新多源遥感数据,保持监测连续性\n")
f.write(" 2. 结合社会经济数据,深入分析变化驱动因素\n")
f.write(" 3. 建立自动化的变化检测和预警系统\n")
f.write(" 4. 加强数据质量控制,确保分析结果可靠性\n")
f.write("\n" + "=" * 80 + "\n")
f.write("报告结束\n")
f.write("=" * 80 + "\n")
print(f"分析报告已保存: {output_path}")
# ==================== 主程序 ====================
def main():
"""主函数"""
print("=" * 80)
print("印度交通基础设施变化监测分析系统")
print("=" * 80)
try:
# 检查必要的库
print("\n检查必要的库...")
required_libs = ['geopandas', 'rasterio', 'matplotlib', 'pandas', 'numpy']
missing = []
for lib in required_libs:
try:
__import__(lib)
except ImportError:
missing.append(lib)
if missing:
print(f"错误: 缺少必要的库: {missing}")
print("请使用以下命令安装:")
print(f"pip install {' '.join(missing)}")
return
# 1. 加载边界数据
print("\n1. 加载边界数据...")
boundary = load_boundary_data()
if boundary is None:
print("错误: 无法加载边界数据,程序终止")
return
# 2. 加载和处理MODIS数据
print("\n2. 加载和处理MODIS数据...")
modis_files, urban_stats, class_distributions = find_and_load_modis_data(MODIS_YEARS, boundary)
if not modis_files:
print("警告: 没有找到有效的MODIS数据")
urban_analysis = None
else:
urban_analysis = analyze_urban_expansion(urban_stats)
# 3. 加载和处理夜间灯光数据
print("\n3. 加载和处理夜间灯光数据...")
nightlight_files, nightlight_stats = find_and_load_nightlight_data(NIGHTLIGHT_YEARS, boundary)
# 4. 加载和处理OSM道路数据
print("\n4. 加载和处理OSM道路数据...")
roads_data = find_and_load_osm_data(OSM_YEARS, boundary)
if roads_data:
road_stats = analyze_road_infrastructure(roads_data, boundary)
else:
road_stats = None
# 5. 创建可视化
print("\n5. 创建可视化...")
# 城市扩张图表
if urban_analysis:
urban_plot_path = os.path.join(OUTPUT_DIR, "urban_expansion_analysis.png")
create_urban_expansion_plot(urban_analysis, urban_plot_path)
# 道路分析图表
if road_stats and roads_data:
road_plot_path = os.path.join(OUTPUT_DIR, "road_infrastructure_analysis.png")
create_road_analysis_plot(road_stats, boundary, roads_data, road_plot_path)
# 夜间灯光图表
if nightlight_stats:
nightlight_plot_path = os.path.join(OUTPUT_DIR, "nightlight_analysis.png")
create_nightlight_analysis_plot(nightlight_stats, nightlight_plot_path)
# 6. 生成报告
print("\n6. 生成分析报告...")
report_path = os.path.join(OUTPUT_DIR, "analysis_report.txt")
generate_analysis_report(urban_analysis, road_stats, nightlight_stats,
class_distributions, report_path)
# 7. 保存中间数据
print("\n7. 保存中间数据...")
# 保存统计数据
if urban_stats:
urban_df = pd.DataFrame.from_dict(urban_stats, orient='index')
urban_csv = os.path.join(OUTPUT_DIR, "urban_statistics.csv")
urban_df.to_csv(urban_csv, index_label='year')
print(f" 城市统计数据: {urban_csv}")
if nightlight_stats:
nightlight_df = pd.DataFrame.from_dict(nightlight_stats, orient='index')
nightlight_csv = os.path.join(OUTPUT_DIR, "nightlight_statistics.csv")
nightlight_df.to_csv(nightlight_csv, index_label='year')
print(f" 夜间灯光统计数据: {nightlight_csv}")
if road_stats:
road_stats_json = os.path.join(OUTPUT_DIR, "road_statistics.json")
with open(road_stats_json, 'w') as f:
json.dump(road_stats, f, indent=2)
print(f" 道路统计数据: {road_stats_json}")
print("\n" + "=" * 80)
print("分析完成!")
print("=" * 80)
print(f"\n主要输出文件:")
print(f"1. 分析报告: {report_path}")
if 'urban_plot_path' in locals():
print(f"2. 城市扩张分析图: {urban_plot_path}")
if 'road_plot_path' in locals():
print(f"3. 道路基础设施分析图: {road_plot_path}")
if 'nightlight_plot_path' in locals():
print(f"4. 夜间灯光分析图: {nightlight_plot_path}")
print(f"\n数据处理统计:")
print(f" • 处理了 {len(modis_files) if 'modis_files' in locals() else 0} 年MODIS数据")
print(f" • 处理了 {len(nightlight_files) if 'nightlight_files' in locals() else 0} 年夜光数据")
print(f" • 分析了 {sum(len(r) for r in roads_data.values()) if roads_data else 0} 条道路")
print(f"\n输出目录: {OUTPUT_DIR}")
print(f"临时目录: {TEMP_DIR}")
except Exception as e:
print(f"\n程序执行出错: {str(e)}")
import traceback
traceback.print_exc()
print("\n故障排除建议:")
print("1. 检查数据文件是否存在且可读")
print("2. 确保有足够的磁盘空间")
print("3. 安装必要的Python库")
print("4. 检查文件路径中的中文字符")
# ==================== 运行程序 ====================
if __name__ == "__main__":
main() 基于以上代码进行调整,使得以下功能完整实现,结合GDAL、leafmap等库,基于MODIS数据完成印度交通基础设施变化的监测、分析和可视化。要求综合所学知识进行设计、处理、分析和可视化,并给出最终结果和结论。我有印度的2015年到2024年的modis数据,2020年和2024年的夜间灯光数据(tif格式),2020年到2025年的osm道路数据,osm道路数据帮我选取主要道路就行了,印度边界数据。这些数据的投影可能不一致。这个代码已经出现了无法加载边界数据的问题