50 years, 50 colors 二分匹配图

本文探讨了一个有趣的算法问题——如何用有限次数的机会踩破所有相同颜色的气球。通过构建一个n*n的矩阵来模拟放置各色气球的情况,并采用匹配算法来找出最少的操作次数。如果无法完成任务,则输出不能被全部踩破的气球颜色。

所有点的x,y之间连一条边,求最小点覆盖

题意:一个n*n矩阵,每个格子放一个气球,气球有颜色。
一个人一次可以选择一种颜色的气球,再选择一行或者一列,
把该种颜色的气球踩破;

你有k次机会,看否把某种颜色的气球全部踩破。若有些颜色的气球不能被踩破
,按从小大的顺序输出这些气球。否则,输出-1.

#include<bits/stdc++.h>
using namespace std;
const int maxn=200;
int n,k;
set<int>sett;
set<int> ::iterator it;
int mapp[maxn][maxn];
int answer[maxn];
int cx[maxn],cy[maxn],inpath[maxn];
int findpath(int u,int xx)
{
    for(int i=1;i<=n;i++)
    {
        int v=i;
        if(!inpath[v]&&mapp[u][v]==xx)
        {
            inpath[v]=1;
            if(cy[v]==-1||findpath(cy[v],xx))
            {
                cx[u]=v;
                cy[v]=u;
                return 1;
            }
        }
    }
    return 0;
}
int maxmatch(int xx)
{
    int res=0;
    for(int i=0;i<=n;i++)
    {
        cx[i]=-1;
        cy[i]=-1;
    }
    for(int i=1;i<=n;i++)
    {
        if(cx[i]==-1)
        {

            memset(inpath,0,sizeof(inpath));
            res+=findpath(i,xx);
            //printf("i:%d %d\n",i,res);
        }
    }
    return res;
}
int main ()
{
    while(~scanf("%d%d",&n,&k))
    {
        if(n==0&&k==0)
        {
            break;
        }
        sett.clear();
        for(int i=1;i<=n;i++)
        {
            for(int j=1;j<=n;j++)
            {
                scanf("%d",&mapp[i][j]);
                sett.insert(mapp[i][j]);
            }
        }
        int sum=0;
        int cnt=0;
        for(it=sett.begin();it!=sett.end();it++)
        {
            int ans=maxmatch(*it);
            //printf("i:%d %d\n",i,ans);
            if(ans>k)
            {
                sum++;
                answer[++cnt]=*it;
            }
        }
        if(sum==0)
        {
            printf("-1\n");
        }
        else
        {
            sort(answer+1,answer+1+cnt);
            for(int i=1;i<cnt;i++)
            {
                printf("%d ",answer[i]);
            }
            printf("%d\n",answer[cnt]);
        }
    }
}
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道路数据帮我选取主要道路就行了,印度边界数据。这些数据的投影可能不一致。这个代码已经出现了无法加载边界数据的问题
12-17
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值