import numpy as np from osgeo import gdal import math import os from typing import Optional, Tuple, Dict, Any class SlopeAspectGenerator: """ 坡度和坡向生成器 支持生成单独的坡度、坡向文件，或包含两者的合并文件 """ def __init__(self, dem_path: str): """ 初始化坡度和坡向生成器 Args: dem_path: 输入DEM文件路径 """ self.dem_path = dem_path self.dataset = gdal.Open(dem_path, gdal.GA_ReadOnly) if not self.dataset: raise FileNotFoundError(f"无法打开DEM文件: {dem_path}") # 获取地理信息 self.geotransform = self.dataset.GetGeoTransform() self.projection = self.dataset.GetProjection() self.band = self.dataset.GetRasterBand(1) self.data = self.band.ReadAsArray() # 获取栅格信息 self.cols = self.dataset.RasterXSize self.rows = self.dataset.RasterYSize self.no_data = self.band.GetNoDataValue() # 如果DEM没有无效值设置，使用-9999 if self.no_data is None: self.no_data = -9999.0 print(f"警告: DEM未设置无效值，将使用 {self.no_data} 作为无效值") # 像元大小（考虑非正方形像元） self.cell_size_x = abs(self.geotransform[1]) self.cell_size_y = abs(self.geotransform[5]) print(f"DEM信息: {self.rows}行 x {self.cols}列") print(f"X方向像元大小: {self.cell_size_x:.4f}") print(f"Y方向像元大小: {self.cell_size_y:.4f}") print(f"无效值: {self.no_data}") def calculate_slope(self, z_factor: float = 1.0, algorithm: str = 'Horn') -> np.ndarray: """ 计算坡度 Args: z_factor: 垂直单位转换因子 algorithm: 算法选项 'Horn' 或 'ZevenbergenThorne' Returns: 坡度数组（单位：度） """ print(f"正在计算坡度 (算法: {algorithm})...") # 创建输出数组 slope = np.full((self.rows, self.cols), self.no_data, dtype=np.float32) # 创建填充后的数据（用于边界处理） padded_data = np.pad( self.data, pad_width=1, mode='constant', constant_values=self.no_data ) # 创建有效掩码 valid_mask = (self.data != self.no_data) padded_valid = np.pad(valid_mask, pad_width=1, mode='constant', constant_values=False) for i in range(self.rows): for j in range(self.cols): if not valid_mask[i, j]: continue i_pad, j_pad = i + 1, j + 1 # 获取3x3窗口 window = padded_data[i_pad-1:i_pad+2, j_pad-1:j_pad+2] window_valid = padded_valid[i_pad-1:i_pad+2, j_pad-1:j_pad+2] # 如果窗口内有无效值，跳过 if not np.all(window_valid): continue if algorithm.upper() == 'ZEVERBERGENTHORNE': # Zevenbergen & Thorne 算法 dz_dx = (window[1, 2] - window[1, 0]) / (2 * self.cell_size_x) dz_dy = (window[0, 1] - window[2, 1]) / (2 * self.cell_size_y) else: # Horn 算法（默认） dz_dx = ((window[0, 2] + 2 * window[1, 2] + window[2, 2]) - (window[0, 0] + 2 * window[1, 0] + window[2, 0])) / (8 * self.cell_size_x) dz_dy = ((window[2, 0] + 2 * window[2, 1] + window[2, 2]) - (window[0, 0] + 2 * window[0, 1] + window[0, 2])) / (8 * self.cell_size_y) # 计算坡度（弧度转角度） slope_rad = math.atan(z_factor * math.sqrt(dz_dx**2 + dz_dy**2)) slope_deg = math.degrees(slope_rad) slope[i, j] = slope_deg # 边界处理：将边缘设为无效值 slope[0, :] = self.no_data slope[-1, :] = self.no_data slope[:, 0] = self.no_data slope[:, -1] = self.no_data valid_slope = slope[slope != self.no_data] if len(valid_slope) > 0: print(f"坡度范围: {np.min(valid_slope):.2f} - {np.max(valid_slope):.2f} 度") return slope def calculate_aspect(self, algorithm: str = 'Horn') -> np.ndarray: """ 计算坡向 Args: algorithm: 算法选项 'Horn' 或 'ZevenbergenThorne' Returns: 坡向数组（单位：度，0-360，北为0度，顺时针增加） """ print(f"正在计算坡向 (算法: {algorithm})...") # 创建输出数组 aspect = np.full((self.rows, self.cols), -1.0, dtype=np.float32) # 平地用-1表示 # 创建填充后的数据 padded_data = np.pad( self.data, pad_width=1, mode='constant', constant_values=self.no_data ) # 创建有效掩码 valid_mask = (self.data != self.no_data) padded_valid = np.pad(valid_mask, pad_width=1, mode='constant', constant_values=False) for i in range(self.rows): for j in range(self.cols): if not valid_mask[i, j]: continue i_pad, j_pad = i + 1, j + 1 # 获取3x3窗口 window = padded_data[i_pad-1:i_pad+2, j_pad-1:j_pad+2] window_valid = padded_valid[i_pad-1:i_pad+2, j_pad-1:j_pad+2] # 如果窗口内有无效值，跳过 if not np.all(window_valid): continue if algorithm.upper() == 'ZEVERBERGENTHORNE': # Zevenbergen & Thorne 算法 dz_dx = (window[1, 2] - window[1, 0]) / (2 * self.cell_size_x) dz_dy = (window[0, 1] - window[2, 1]) / (2 * self.cell_size_y) else: # Horn 算法（默认） dz_dx = ((window[0, 2] + 2 * window[1, 2] + window[2, 2]) - (window[0, 0] + 2 * window[1, 0] + window[2, 0])) / (8 * self.cell_size_x) dz_dy = ((window[2, 0] + 2 * window[2, 1] + window[2, 2]) - (window[0, 0] + 2 * window[0, 1] + window[0, 2])) / (8 * self.cell_size_y) # 计算坡向 if abs(dz_dx) < 1e-10 and abs(dz_dy) < 1e-10: # 平地 aspect[i, j] = -1.0 else: aspect_rad = math.atan2(dz_dy, -dz_dx) # atan2(y, x) aspect_deg = math.degrees(aspect_rad) # 转换为0-360度（北为0度，顺时针） if aspect_deg < 0: aspect_deg += 360.0 aspect[i, j] = aspect_deg # 边界处理：将边缘设为-1 aspect[0, :] = -1.0 aspect[-1, :] = -1.0 aspect[:, 0] = -1.0 aspect[:, -1] = -1.0 valid_aspect = aspect[aspect != -1.0] if len(valid_aspect) > 0: print(f"坡向范围: {np.min(valid_aspect):.2f} - {np.max(valid_aspect):.2f} 度") return aspect def calculate_slope_aspect(self, z_factor: float = 1.0, algorithm: str = 'Horn') -> Tuple[np.ndarray, np.ndarray]: """ 同时计算坡度和坡向（优化版本，减少重复计算） Args: z_factor: 垂直单位转换因子 algorithm: 算法选项 Returns: (slope_array, aspect_array) """ print(f"同时计算坡度和坡向 (算法: {algorithm})...") # 创建输出数组 slope = np.full((self.rows, self.cols), self.no_data, dtype=np.float32) aspect = np.full((self.rows, self.cols), -1.0, dtype=np.float32) # 创建填充后的数据 padded_data = np.pad( self.data, pad_width=1, mode='constant', constant_values=self.no_data ) # 创建有效掩码 valid_mask = (self.data != self.no_data) padded_valid = np.pad(valid_mask, pad_width=1, mode='constant', constant_values=False) for i in range(self.rows): for j in range(self.cols): if not valid_mask[i, j]: continue i_pad, j_pad = i + 1, j + 1 # 获取3x3窗口 window = padded_data[i_pad-1:i_pad+2, j_pad-1:j_pad+2] window_valid = padded_valid[i_pad-1:i_pad+2, j_pad-1:j_pad+2] # 如果窗口内有无效值，跳过 if not np.all(window_valid): continue # 计算导数（一次计算，两个都用） if algorithm.upper() == 'ZEVERBERGENTHORNE': dz_dx = (window[1, 2] - window[1, 0]) / (2 * self.cell_size_x) dz_dy = (window[0, 1] - window[2, 1]) / (2 * self.cell_size_y) else: dz_dx = ((window[0, 2] + 2 * window[1, 2] + window[2, 2]) - (window[0, 0] + 2 * window[1, 0] + window[2, 0])) / (8 * self.cell_size_x) dz_dy = ((window[2, 0] + 2 * window[2, 1] + window[2, 2]) - (window[0, 0] + 2 * window[0, 1] + window[0, 2])) / (8 * self.cell_size_y) # 计算坡度 slope_rad = math.atan(z_factor * math.sqrt(dz_dx**2 + dz_dy**2)) slope_deg = math.degrees(slope_rad) slope[i, j] = slope_deg # 计算坡向 if abs(dz_dx) < 1e-10 and abs(dz_dy) < 1e-10: aspect[i, j] = -1.0 else: aspect_rad = math.atan2(dz_dy, -dz_dx) aspect_deg = math.degrees(aspect_rad) if aspect_deg < 0: aspect_deg += 360.0 aspect[i, j] = aspect_deg # 边界处理 slope[0, :] = self.no_data slope[-1, :] = self.no_data slope[:, 0] = self.no_data slope[:, -1] = self.no_data aspect[0, :] = -1.0 aspect[-1, :] = -1.0 aspect[:, 0] = -1.0 aspect[:, -1] = -1.0 # 显示统计信息 valid_slope = slope[slope != self.no_data] valid_aspect = aspect[aspect != -1.0] if len(valid_slope) > 0: print(f"坡度范围: {np.min(valid_slope):.2f} - {np.max(valid_slope):.2f} 度") if len(valid_aspect) > 0: print(f"坡向范围: {np.min(valid_aspect):.2f} - {np.max(valid_aspect):.2f} 度") return slope, aspect def save_raster(self, data: np.ndarray, output_path: str, band_name: str = '', data_type: int = gdal.GDT_Float32, no_data_value: Optional[float] = None) -> None: """ 保存栅格数据为TIFF文件 Args: data: 栅格数据数组 output_path: 输出文件路径 band_name: 波段名称 data_type: GDAL数据类型 no_data_value: 无效值 """ driver = gdal.GetDriverByName('GTiff') out_dataset = driver.Create( output_path, self.cols, self.rows, 1, # 单波段 data_type, options=['COMPRESS=LZW', 'PREDICTOR=2', 'TILED=YES'] ) out_dataset.SetGeoTransform(self.geotransform) out_dataset.SetProjection(self.projection) out_band = out_dataset.GetRasterBand(1) out_band.WriteArray(data) if no_data_value is not None: out_band.SetNoDataValue(no_data_value) elif band_name == 'Slope': out_band.SetNoDataValue(self.no_data) elif band_name == 'Aspect': out_band.SetNoDataValue(-1.0) if band_name: out_band.SetDescription(band_name) # 设置统计信息（加速QGIS等软件加载） out_band.ComputeStatistics(False) out_dataset.FlushCache() out_dataset = None print(f"文件已保存: {output_path}") if band_name: print(f" 波段: {band_name}") def save_combined_raster(self, slope_data: np.ndarray, aspect_data: np.ndarray, output_path: str) -> None: """ 保存包含坡度和坡向两个波段的TIFF文件 Args: slope_data: 坡度数据 aspect_data: 坡向数据 output_path: 输出文件路径 """ driver = gdal.GetDriverByName('GTiff') # 创建包含2个波段的数据集 out_dataset = driver.Create( output_path, self.cols, self.rows, 2, # 2个波段 gdal.GDT_Float32, options=['COMPRESS=LZW', 'PREDICTOR=2', 'TILED=YES'] ) # 设置地理信息 out_dataset.SetGeoTransform(self.geotransform) out_dataset.SetProjection(self.projection) # 写入坡度波段（波段1） slope_band = out_dataset.GetRasterBand(1) slope_band.WriteArray(slope_data) slope_band.SetNoDataValue(self.no_data) slope_band.SetDescription('Slope') slope_band.SetMetadataItem('UNITS', 'degrees') slope_band.SetMetadataItem('RANGE', '0-90') slope_band.ComputeStatistics(False) # 写入坡向波段（波段2） aspect_band = out_dataset.GetRasterBand(2) aspect_band.WriteArray(aspect_data) aspect_band.SetNoDataValue(-1.0) aspect_band.SetDescription('Aspect') aspect_band.SetMetadataItem('UNITS', 'degrees') aspect_band.SetMetadataItem('RANGE', '0-360') aspect_band.SetMetadataItem('DIRECTION', 'Clockwise from North') aspect_band.ComputeStatistics(False) # 设置数据集元数据 out_dataset.SetMetadata({ 'SOURCE_DEM': os.path.basename(self.dem_path), 'PROCESSING_METHOD': 'SlopeAspectGenerator', 'CELL_SIZE_X': str(self.cell_size_x), 'CELL_SIZE_Y': str(self.cell_size_y) }) out_dataset.FlushCache() out_dataset = None print(f"合并文件已保存: {output_path}") print(f" 波段1: Slope (坡度, 单位:度)") print(f" 波段2: Aspect (坡向, 单位:度)") def generate_slope_tif(self, output_path: Optional[str] = None, z_factor: float = 1.0, algorithm: str = 'Horn') -> np.ndarray: """ 生成坡度TIFF文件 Args: output_path: 输出文件路径（可选） z_factor: 垂直单位转换因子 algorithm: 算法选项 Returns: 坡度数据数组 """ if output_path is None: base_name = os.path.splitext(self.dem_path)[0] output_path = f"{base_name}_slope.tif" print(f"\n生成坡度文件...") slope = self.calculate_slope(z_factor, algorithm) self.save_raster(slope, output_path, 'Slope', no_data_value=self.no_data) return slope def generate_aspect_tif(self, output_path: Optional[str] = None, algorithm: str = 'Horn') -> np.ndarray: """ 生成坡向TIFF文件 Args: output_path: 输出文件路径（可选） algorithm: 算法选项 Returns: 坡向数据数组 """ if output_path is None: base_name = os.path.splitext(self.dem_path)[0] output_path = f"{base_name}_aspect.tif" print(f"\n生成坡向文件...") aspect = self.calculate_aspect(algorithm) self.save_raster(aspect, output_path, 'Aspect', no_data_value=-1.0) return aspect def generate_slope_aspect_separate(self, slope_output: Optional[str] = None, aspect_output: Optional[str] = None, z_factor: float = 1.0, algorithm: str = 'Horn') -> Tuple[np.ndarray, np.ndarray]: """ 分别生成坡度和坡向TIFF文件 Args: slope_output: 坡度输出路径（可选） aspect_output: 坡向输出路径（可选） z_factor: 垂直单位转换因子 algorithm: 算法选项 Returns: (slope_data, aspect_data) """ # 设置输出路径 base_name = os.path.splitext(self.dem_path)[0] if slope_output is None: slope_output = f"{base_name}_slope.tif" if aspect_output is None: aspect_output = f"{base_name}_aspect.tif" print(f"\n分别生成坡度和坡向文件...") # 同时计算坡度和坡向（优化版本） slope, aspect = self.calculate_slope_aspect(z_factor, algorithm) # 保存文件 self.save_raster(slope, slope_output, 'Slope', no_data_value=self.no_data) self.save_raster(aspect, aspect_output, 'Aspect', no_data_value=-1.0) return slope, aspect def generate_combined_tif(self, output_path: Optional[str] = None, z_factor: float = 1.0, algorithm: str = 'Horn') -> Tuple[np.ndarray, np.ndarray]: """ 生成包含坡度和坡向的合并TIFF文件 Args: output_path: 输出文件路径（可选） z_factor: 垂直单位转换因子 algorithm: 算法选项 Returns: (slope_data, aspect_data) """ if output_path is None: base_name = os.path.splitext(self.dem_path)[0] output_path = f"{base_name}_slope_aspect.tif" print(f"\n生成合并文件...") # 同时计算坡度和坡向 slope, aspect = self.calculate_slope_aspect(z_factor, algorithm) # 保存合并文件 self.save_combined_raster(slope, aspect, output_path) return slope, aspect def generate_all(self, output_dir: Optional[str] = None, z_factor: float = 1.0, algorithm: str = 'Horn') -> Dict[str, str]: """ 生成所有类型的文件：单独的坡度、单独的坡向、合并文件、统计信息 Args: output_dir: 输出目录（可选） z_factor: 垂直单位转换因子 algorithm: 算法选项 Returns: 包含所有输出文件路径的字典 """ # 设置输出路径 base_name = os.path.basename(self.dem_path) file_base = os.path.splitext(base_name)[0] if output_dir: if not os.path.exists(output_dir): os.makedirs(output_dir) output_prefix = os.path.join(output_dir, file_base) else: output_prefix = file_base # 计算坡度和坡向（只计算一次） print(f"\n计算坡度和坡向...") slope, aspect = self.calculate_slope_aspect(z_factor, algorithm) # 生成各种文件 results = {} # 1. 单独的坡度文件 slope_file = f"{output_prefix}_slope.tif" self.save_raster(slope, slope_file, 'Slope', no_data_value=self.no_data) results['slope'] = slope_file # 2. 单独的坡向文件 aspect_file = f"{output_prefix}_aspect.tif" self.save_raster(aspect, aspect_file, 'Aspect', no_data_value=-1.0) results['aspect'] = aspect_file # 3. 合并文件 combined_file = f"{output_prefix}_slope_aspect.tif" self.save_combined_raster(slope, aspect, combined_file) results['combined'] = combined_file # 4. 统计信息文件（文本） stats_file = f"{output_prefix}_statistics.txt" self.save_statistics(slope, aspect, stats_file) results['statistics'] = stats_file print(f"\n所有文件生成完成！") for key, path in results.items(): print(f" {key}: {os.path.basename(path)}") return results def save_statistics(self, slope_data: np.ndarray, aspect_data: np.ndarray, output_path: str) -> None: """ 保存统计信息到文本文件 Args: slope_data: 坡度数据 aspect_data: 坡向数据 output_path: 输出文件路径 """ valid_slope = slope_data[slope_data != self.no_data] valid_aspect = aspect_data[aspect_data != -1.0] with open(output_path, 'w', encoding='utf-8') as f: f.write("坡度和坡向统计信息\n") f.write("=" * 50 + "\n") f.write(f"输入DEM: {os.path.basename(self.dem_path)}\n") f.write(f"栅格大小: {self.rows} x {self.cols}\n") f.write(f"处理时间: {np.datetime64('now')}\n\n") f.write("坡度统计:\n") f.write("-" * 30 + "\n") if len(valid_slope) > 0: f.write(f"最小值: {np.min(valid_slope):.4f} 度\n") f.write(f"最大值: {np.max(valid_slope):.4f} 度\n") f.write(f"平均值: {np.mean(valid_slope):.4f} 度\n") f.write(f"标准差: {np.std(valid_slope):.4f} 度\n") f.write(f"中位数: {np.median(valid_slope):.4f} 度\n") f.write(f"有效像元数: {len(valid_slope)} / {self.rows * self.cols}\n") # 坡度分级统计 f.write("\n坡度分级统计:\n") slope_classes = ['平地(0-2)', '缓坡(2-5)', '斜坡(5-15)', '陡坡(15-30)', '急陡坡(30-45)', '峭壁(45-90)'] bins = [0, 2, 5, 15, 30, 45, 90] hist, _ = np.histogram(valid_slope, bins=bins) total = sum(hist) for i in range(len(slope_classes)): if hist[i] > 0: percentage = (hist[i] / total) * 100 f.write(f"{slope_classes[i]}: {hist[i]} ({percentage:.2f}%)\n") else: f.write("无有效数据\n") f.write("\n坡向统计:\n") f.write("-" * 30 + "\n") if len(valid_aspect) > 0: f.write(f"最小值: {np.min(valid_aspect):.4f} 度\n") f.write(f"最大值: {np.max(valid_aspect):.4f} 度\n") f.write(f"平均值: {np.mean(valid_aspect):.4f} 度\n") f.write(f"标准差: {np.std(valid_aspect):.4f} 度\n") f.write(f"中位数: {np.median(valid_aspect):.4f} 度\n") f.write(f"有效像元数: {len(valid_aspect)} / {self.rows * self.cols}\n") # 坡向分类统计 f.write("\n坡向分类统计:\n") directions = ['北(0-22.5)', '东北(22.5-67.5)', '东(67.5-112.5)', '东南(112.5-157.5)', '南(157.5-202.5)', '西南(202.5-247.5)', '西(247.5-292.5)', '西北(292.5-337.5)', '北(337.5-360)'] bins = [0, 22.5, 67.5, 112.5, 157.5, 202.5, 247.5, 292.5, 337.5, 360] hist, _ = np.histogram(valid_aspect, bins=bins) total = sum(hist) for i in range(len(directions)): if hist[i] > 0: percentage = (hist[i] / total) * 100 f.write(f"{directions[i]}: {hist[i]} ({percentage:.2f}%)\n") else: f.write("无有效数据\n") print(f"统计信息已保存: {output_path}") def get_dem_info(self) -> Dict[str, Any]: """ 获取DEM文件的详细信息 Returns: 包含DEM信息的字典 """ info = { 'file_path': self.dem_path, 'file_name': os.path.basename(self.dem_path), 'rows': self.rows, 'cols': self.cols, 'cell_size_x': self.cell_size_x, 'cell_size_y': self.cell_size_y, 'no_data_value': self.no_data, 'projection': self.projection } # 获取高程统计 valid_data = self.data[self.data != self.no_data] if len(valid_data) > 0: info['elevation_min'] = float(np.min(valid_data)) info['elevation_max'] = float(np.max(valid_data)) info['elevation_mean'] = float(np.mean(valid_data)) info['elevation_std'] = float(np.std(valid_data)) info['valid_cells'] = int(len(valid_data)) info['total_cells'] = self.rows * self.cols info['valid_percentage'] = (len(valid_data) / (self.rows * self.cols)) * 100 return info def get_slope_statistics(self, slope_data: np.ndarray) -> Dict[str, Any]: """ 获取坡度统计信息 Args: slope_data: 坡度数据数组 Returns: 坡度统计信息字典 """ valid_slope = slope_data[slope_data != self.no_data] stats = { 'min': float(np.min(valid_slope)) if len(valid_slope) > 0 else None, 'max': float(np.max(valid_slope)) if len(valid_slope) > 0 else None, 'mean': float(np.mean(valid_slope)) if len(valid_slope) > 0 else None, 'std': float(np.std(valid_slope)) if len(valid_slope) > 0 else None, 'median': float(np.median(valid_slope)) if len(valid_slope) > 0 else None, 'valid_cells': int(len(valid_slope)), 'total_cells': self.rows * self.cols } return stats def get_aspect_statistics(self, aspect_data: np.ndarray) -> Dict[str, Any]: """ 获取坡向统计信息 Args: aspect_data: 坡向数据数组 Returns: 坡向统计信息字典 """ valid_aspect = aspect_data[aspect_data != -1.0] stats = { 'min': float(np.min(valid_aspect)) if len(valid_aspect) > 0 else None, 'max': float(np.max(valid_aspect)) if len(valid_aspect) > 0 else None, 'mean': float(np.mean(valid_aspect)) if len(valid_aspect) > 0 else None, 'std': float(np.std(valid_aspect)) if len(valid_aspect) > 0 else None, 'median': float(np.median(valid_aspect)) if len(valid_aspect) > 0 else None, 'valid_cells': int(len(valid_aspect)), 'total_cells': self.rows * self.cols } # 坡向分类统计 if len(valid_aspect) > 0: bins = [0, 22.5, 67.5, 112.5, 157.5, 202.5, 247.5, 292.5, 337.5, 360] directions = ['北', '东北', '东', '东南', '南', '西南', '西', '西北', '北'] hist, _ = np.histogram(valid_aspect, bins=bins) direction_stats = {} for i in range(len(directions)): direction_stats[directions[i]] = { 'count': int(hist[i]), 'percentage': (hist[i] / len(valid_aspect)) * 100 if len(valid_aspect) > 0 else 0 } stats['direction_distribution'] = direction_stats return stats def close(self): """关闭数据集，释放资源""" if hasattr(self, 'dataset') and self.dataset: self.dataset = None if hasattr(self, 'band') and self.band: self.band = None print("资源已释放") def __enter__(self): """支持with语句""" return self def __exit__(self, exc_type, exc_val, exc_tb): """退出with语句时自动关闭""" self.close() # 使用示例和辅助函数 def create_slope_aspect(dem_path: str, output_type: str = 'combined', output_path: Optional[str] = None, z_factor: float = 1.0, algorithm: str = 'Horn') -> Dict[str, str]: """ 创建坡度坡向的便捷函数 Args: dem_path: DEM文件路径 output_type: 输出类型 'slope', 'aspect', 'separate', 'combined', 'all' output_path: 输出路径（对于单个文件）或输出目录 z_factor: 垂直单位转换因子 algorithm: 算法选项 Returns: 包含输出文件路径的字典 """ print(f"处理DEM文件: {dem_path}") with SlopeAspectGenerator(dem_path) as generator: # 获取DEM信息 dem_info = generator.get_dem_info() print(f"DEM信息: {dem_info['rows']}x{dem_info['cols']}, " f"高程范围: {dem_info.get('elevation_min', 'N/A'):.2f}-" f"{dem_info.get('elevation_max', 'N/A'):.2f}") results = {} if output_type == 'slope': # 只生成坡度 slope_file = output_path or f"{os.path.splitext(dem_path)[0]}_slope.tif" slope = generator.generate_slope_tif(slope_file, z_factor, algorithm) results['slope'] = slope_file # 获取统计信息 slope_stats = generator.get_slope_statistics(slope) print(f"坡度统计: {slope_stats['min']:.2f}度 - {slope_stats['max']:.2f}度") elif output_type == 'aspect': # 只生成坡向 aspect_file = output_path or f"{os.path.splitext(dem_path)[0]}_aspect.tif" aspect = generator.generate_aspect_tif(aspect_file, algorithm) results['aspect'] = aspect_file # 获取统计信息 aspect_stats = generator.get_aspect_statistics(aspect) print(f"坡向统计: {aspect_stats['min']:.2f}度 - {aspect_stats['max']:.2f}度") elif output_type == 'separate': # 分别生成坡度和坡向 if output_path and os.path.isdir(output_path): # 如果输出路径是目录 slope_file = os.path.join(output_path, f"{os.path.splitext(os.path.basename(dem_path))[0]}_slope.tif") aspect_file = os.path.join(output_path, f"{os.path.splitext(os.path.basename(dem_path))[0]}_aspect.tif") slope, aspect = generator.generate_slope_aspect_separate( slope_file, aspect_file, z_factor, algorithm) else: # 如果输出路径是文件前缀 base = output_path or os.path.splitext(dem_path)[0] slope_file = f"{base}_slope.tif" aspect_file = f"{base}_aspect.tif" slope, aspect = generator.generate_slope_aspect_separate( slope_file, aspect_file, z_factor, algorithm) results['slope'] = slope_file results['aspect'] = aspect_file # 获取统计信息 slope_stats = generator.get_slope_statistics(slope) aspect_stats = generator.get_aspect_statistics(aspect) print(f"坡度范围: {slope_stats['min']:.2f}度 - {slope_stats['max']:.2f}度") print(f"坡向范围: {aspect_stats['min']:.2f}度 - {aspect_stats['max']:.2f}度") elif output_type == 'combined': # 生成合并文件 combined_file = output_path or f"{os.path.splitext(dem_path)[0]}_slope_aspect.tif" slope, aspect = generator.generate_combined_tif(combined_file, z_factor, algorithm) results['combined'] = combined_file # 获取统计信息 slope_stats = generator.get_slope_statistics(slope) aspect_stats = generator.get_aspect_statistics(aspect) print(f"坡度范围: {slope_stats['min']:.2f}度 - {slope_stats['max']:.2f}度") print(f"坡向范围: {aspect_stats['min']:.2f}度 - {aspect_stats['max']:.2f}度") elif output_type == 'all': # 生成所有文件 output_dir = output_path if (output_path and os.path.isdir(output_path)) else None results = generator.generate_all(output_dir, z_factor, algorithm) else: raise ValueError(f"不支持的输出类型: {output_type}") return results def batch_process_slope_aspect(dem_files: list, output_dir: str, output_type: str = 'all', z_factor: float = 1.0, algorithm: str = 'Horn') -> Dict[str, list]: """ 批量处理多个DEM文件 Args: dem_files: DEM文件路径列表 output_dir: 输出目录 output_type: 输出类型 z_factor: 垂直单位转换因子 algorithm: 算法选项 Returns: 包含处理结果的字典 """ if not os.path.exists(output_dir): os.makedirs(output_dir) results = { 'success': [], 'failed': [], 'outputs': [] } total_files = len(dem_files) print(f"开始批量处理 {total_files} 个DEM文件...") for i, dem_file in enumerate(dem_files, 1): if not os.path.exists(dem_file): print(f"[{i}/{total_files}] 文件不存在: {dem_file}") results['failed'].append(dem_file) continue print(f"\n[{i}/{total_files}] 处理: {os.path.basename(dem_file)}") try: # 处理文件 file_results = create_slope_aspect( dem_file, output_type, output_dir, z_factor, algorithm ) results['success'].append(dem_file) results['outputs'].append(file_results) print(f"[{i}/{total_files}] 完成") except Exception as e: print(f"[{i}/{total_files}] 失败: {e}") results['failed'].append(dem_file) # 打印总结 print(f"\n" + "="*50) print(f"批量处理完成!") print(f"成功: {len(results['success'])} 个") print(f"失败: {len(results['failed'])} 个") if results['failed']: print("\n失败的文件:") for failed_file in results['failed']: print(f" {os.path.basename(failed_file)}") return results # 主示例代码 def main_example(): """使用示例""" # 示例1: 基本使用 print("示例1: 基本使用") print("=" * 50) # 配置参数 SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) dem_file = os.path.join(SCRIPT_DIR, "o_dem_23d110d0.tif") if os.path.exists(dem_file): # 方法1: 使用便捷函数 print("\n方法1: 使用便捷函数") slope_file = os.path.join(SCRIPT_DIR, "o_combined.tif") results = create_slope_aspect(dem_file, 'combined', slope_file) # 方法2: 直接使用类 # print("\n方法2: 直接使用类") # with SlopeAspectGenerator(dem_file) as gen: # # 获取DEM信息 # dem_info = gen.get_dem_info() # print(f"DEM信息: {dem_info}") # # 生成合并文件 # slope, aspect = gen.generate_combined_tif("o_combined.tif") # # 获取统计信息 # slope_stats = gen.get_slope_statistics(slope) # aspect_stats = gen.get_aspect_statistics(aspect) # print(f"坡度统计: {slope_stats}") # print(f"坡向统计: {aspect_stats}") else: print(f"DEM文件不存在: {dem_file}") print("请创建一个示例DEM文件或使用您自己的DEM文件") # 创建一个简单的示例DEM print("\n创建一个示例DEM文件用于测试...") create_sample_dem("sample_dem.tif") # 使用示例DEM print("\n使用示例DEM进行测试...") results = create_slope_aspect("sample_dem.tif", 'all') # 示例2: 批量处理 # print("\n\n示例2: 批量处理") # print("=" * 50) # # 假设有多个DEM文件 # dem_files = ["dem1.tif", "dem2.tif", "dem3.tif"] # 替换为实际文件 # # 过滤出存在的文件 # existing_files = [f for f in dem_files if os.path.exists(f)] # if len(existing_files) > 0: # batch_results = batch_process_slope_aspect( # existing_files, # "batch_output", # 'all' # ) # else: # print("未找到DEM文件进行批量处理") # print("创建示例文件进行批量处理演示...") # # 创建几个示例DEM文件 # for i in range(1, 4): # create_sample_dem(f"sample_dem_{i}.tif", offset=i*100) # sample_files = [f"sample_dem_{i}.tif" for i in range(1, 4)] # batch_results = batch_process_slope_aspect( # sample_files, # "batch_output", # 'all' # ) def create_sample_dem(output_path: str, rows: int = 100, cols: int = 100, cell_size: float = 30.0, offset: float = 0.0): """ 创建一个示例DEM文件用于测试 Args: output_path: 输出文件路径 rows: 行数 cols: 列数 cell_size: 像元大小 offset: 高程偏移量 """ from osgeo import osr # 创建一个简单的地形：一个山峰 x = np.linspace(-5, 5, cols) y = np.linspace(-5, 5, rows) X, Y = np.meshgrid(x, y) # 创建高程数据（高斯山峰） elevation = 1000 * np.exp(-(X**2 + Y**2) / 10) + offset # 添加一些噪声 elevation += np.random.normal(0, 10, (rows, cols)) # 创建栅格文件 driver = gdal.GetDriverByName('GTiff') dataset = driver.Create( output_path, cols, rows, 1, gdal.GDT_Float32 ) # 设置地理变换（左下角坐标为(0,0)） dataset.SetGeoTransform((0, cell_size, 0, 0, 0, -cell_size)) # 设置投影（WGS84） srs = osr.SpatialReference() srs.ImportFromEPSG(4326) dataset.SetProjection(srs.ExportToWkt()) # 写入数据 band = dataset.GetRasterBand(1) band.WriteArray(elevation) band.SetNoDataValue(-9999.0) band.SetDescription('Elevation') dataset.FlushCache() dataset = None print(f"示例DEM已创建: {output_path}") if __name__ == "__main__": # 运行示例 main_example()