#! /usr/bin/env python import os import sys import numpy as np import argparse import subprocess import glob from skimage import io # 尝试导入 NASADEM 库 try: from NASADEM import NASADEM from rasters import RasterGrid HAS_NASADEM = True except ImportError: HAS_NASADEM = False print("Warning: NASADEM library not installed. NASADEM download will not be available.") print("Install with: pip install nasadem rasters") # 尝试导入 srtm 库 try: import srtm HAS_SRTM = True except ImportError: HAS_SRTM = False print("Warning: srtm library not installed. SRTM download will not be available.") print("Install with: pip install srtm (requires Python >= 3.12)") print(f"Current Python version: {sys.version_info.major}.{sys.version_info.minor}") resolutions = 30 # 90 def write_demrsc_file(FILE, Corner_LON, Corner_LAT, X_STEP, Y_STEP, WIDTH, LENGTH): """Write ROI_PAC DEM resource file (.dem.rsc) Args: FILE (str): Output resource file path Corner_LON (str): Corner longitude Corner_LAT (str): Corner latitude X_STEP (str): Longitude step size Y_STEP (str): Latitude step size WIDTH (str): Number of columns LENGTH (str): Number of lines """ f = open(FILE, 'w') f.write('DATE12 111111-222222\n') f.write('FILE_LENGTH ' + str(int(LENGTH)) + '\n') f.write('FILE_TYPE .dem\n') f.write('PROCESSOR roipac\n') f.write('PROJECTION LATLON\n') f.write('RLOOKS 1\n') f.write('WIDTH ' + str(int(WIDTH)) + '\n') f.write('XMAX ' + str(int(int(WIDTH)-1)) + '\n') f.write('XMIN 0\n') f.write('X_FIRST ' + str(float(Corner_LON)) + '\n') f.write('X_STEP ' + str(float(X_STEP)) + '\n') f.write('X_UNIT degrees\n') f.write('YMAX ' + str(int(int(LENGTH)-1)) + '\n') f.write('YMIN 0\n') f.write('Y_FIRST ' + str(float(Corner_LAT)) + '\n') f.write('Y_STEP ' + str(float(Y_STEP)) + '\n') f.write('Y_UNIT degrees\n') f.write('Z_OFFSET 0\n') f.write('Z_SCALE 1\n') f.close() def write_dempar_file(FILE, Corner_LON, Corner_LAT, X_STEP, Y_STEP, WIDTH, LENGTH, DATA_FORMAT): """Write Gamma DEM parameter file (.dem.par) Args: FILE (str): Output parameter file path Corner_LON (str): Corner longitude Corner_LAT (str): Corner latitude X_STEP (str): Longitude step size Y_STEP (str): Latitude step size WIDTH (str): Number of columns LENGTH (str): Number of lines DATA_FORMAT (str): Data format (INTEGER*2 or REAL*4) """ DEM_TYPE = 'Copernicus30' Proj = 'EQA' f = open(FILE, 'w') f.write("Gamma DIFF&GEO DEM/MAP parameter file\n") f.write("title:\tIMPORTED DEM FROM %s\n" % DEM_TYPE) f.write("DEM_projection: %s\n" % Proj) f.write("data_format: %s\n" % DATA_FORMAT) f.write("DEM_hgt_offset: 0.00000\n") f.write("DEM_scale: 1.00000\n") f.write("width: %s\n" % WIDTH) f.write("nlines: %s\n" % LENGTH) f.write("corner_lat: %s decimal degrees\n" % Corner_LAT) f.write("corner_lon: %s decimal degrees\n" % Corner_LON) f.write("post_lat: %s decimal degrees\n" % Y_STEP) f.write("post_lon: %s decimal degrees\n" % X_STEP) f.write("\n") f.write("ellipsoid_name: WGS 84\n") f.write("ellipsoid_ra: 6378137.000 m\n") f.write("ellipsoid_reciprocal_flattening: 298.2572236\n") f.write("\n") f.write("datum_name: WGS 1984\n") f.write("datum_shift_dx: 0.000 m\n") f.write("datum_shift_dy: 0.000 m\n") f.write("datum_shift_dz: 0.000 m\n") f.write("datum_scale_m: 0.00000e+00\n") f.write("datum_rotation_alpha: 0.00000e+00 arc-sec\n") f.write("datum_rotation_beta: 0.00000e+00 arc-sec\n") f.write("datum_rotation_gamma: 0.00000e+00 arc-sec\n") f.write("datum_country_list Global Definition, WGS84, World\n") f.write("\n") f.close() def convert_to_gamma(input_file, output_name, byteorder='big', processor='gamma'): """Convert TIF file to Gamma or ROI_PAC format (.dem and .dem.par/.dem.rsc) Args: input_file (str): Input TIF file path output_name (str): Output name (without extension) byteorder (str): Byte order ('big' or 'little') processor (str): Processor type ('gamma' or 'roi_pac') """ processor_name = 'Gamma' if processor == 'gamma' else 'ROI_PAC' print(f"\n开始转换到 {processor_name} 格式: {input_file}") # 读取 DEM 数据 - 优先使用 GDAL（更稳定），然后是 rasterio dem_data = None # 方法1: 使用 GDAL（最稳定） try: from osgeo import gdal print(" 尝试使用 GDAL 读取数据...") ds = gdal.Open(input_file) if ds is not None: band = ds.GetRasterBand(1) dem_data = band.ReadAsArray() print(f" ✓ 使用 GDAL 读取成功") ds = None # 关闭文件 except Exception as e: print(f" ✗ GDAL 读取失败: {e}") # 方法2: 如果GDAL失败，尝试rasterio if dem_data is None: try: import rasterio print(" 尝试使用 rasterio 读取数据...") with rasterio.open(input_file) as src: dem_data = src.read(1) print(f" ✓ 使用 rasterio 读取成功") except Exception as e: print(f" ✗ rasterio 读取失败: {e}") # 方法3: 最后尝试skimage if dem_data is None: try: print(" 尝试使用 skimage 读取数据...") dem_data = io.imread(input_file) print(f" ✓ 使用 skimage 读取成功") except Exception as e: print(f" ✗ skimage 读取失败: {e}") # 如果所有方法都失败 if dem_data is None: raise ValueError(f"无法读取DEM文件: {input_file}，所有读取方法都失败了") # 确定数据格式 if dem_data.dtype == 'float32': DATA_FORMAT = 'REAL*4' else: DATA_FORMAT = 'INTEGER*2' # 字节序转换 if not sys.byteorder == byteorder: dem_data.byteswap(True) # 输出文件路径 dem_file = output_name + '.dem' if processor == 'gamma': dem_par_file = output_name + '.dem.par' else: dem_par_file = output_name + '.dem.rsc' # 写入二进制 DEM 数据 dem_data.tofile(dem_file) # 使用 gdalinfo 获取地理信息 info_file = 'temp_gdalinfo.txt' cmd = f"gdalinfo {input_file} > {info_file}" os.system(cmd) # 解析地理信息 Corner_LON = None Corner_LAT = None Post_LON = None Post_LAT = None WIDTH = None FILE_LENGTH = None with open(info_file, 'r') as f: for line in f: if 'Origin =' in line: parts = line.split('Origin =')[1].strip().split('(')[1].split(')')[0].split(',') Corner_LON = parts[0] Corner_LAT = parts[1] elif 'Pixel Size ' in line: parts = line.split('Pixel Size =')[1].strip().split('(')[1].split(')')[0].split(',') Post_LON = parts[0] Post_LAT = parts[1] elif 'Size is' in line: parts = line.split('Size is')[1].strip().split(',') WIDTH = parts[0] FILE_LENGTH = parts[1] # 删除临时文件 if os.path.exists(info_file): os.remove(info_file) # 写入参数文件 if processor == 'gamma': write_dempar_file(dem_par_file, Corner_LON, Corner_LAT, Post_LON, Post_LAT, WIDTH, FILE_LENGTH, DATA_FORMAT) else: write_demrsc_file(dem_par_file, Corner_LON, Corner_LAT, Post_LON, Post_LAT, WIDTH, FILE_LENGTH) print(f"{byteorder} endian {dem_file} and {dem_par_file} are generated.") print(f"{processor_name} 格式转换完成!") def get_sufix(STR): """Get file extension""" n = len(STR.split('.')) SUFIX = STR.split('.')[n-1] return SUFIX def read_region(STR): """Parse region string 'west/east/south/north'""" WEST = STR.split('/')[0] EAST = STR.split('/')[1].split('/')[0] SOUTH = STR.split(EAST+'/')[1].split('/')[0] NORTH = STR.split(EAST+'/')[1].split('/')[1] WEST = float(WEST) SOUTH = float(SOUTH) EAST = float(EAST) NORTH = float(NORTH) return WEST, SOUTH, EAST, NORTH def cmd_init(lon, lat, save_path): s_lon = str(abs(lon)) s_lat = str(abs(lat)) if abs(lon) < 10: s_lon = "00" + str(abs(lon)) elif abs(lon) < 100: s_lon = "0" + str(abs(lon)) if abs(lat) < 10: s_lat = "0" + str(abs(lat)) if lon < 0: c_lon = "W" + str(s_lon) else: c_lon = "E" + str(s_lon) if lat < 0: c_lat = "S" + str(s_lat) else: c_lat = "N" + str(s_lat) cmd = "aws s3 cp --no-sign-request" + " s3://copernicus-dem-{0}m/Copernicus_DSM_COG_{1}_{2}_00_{3}_00_DEM/ {4} --recursive".format( str(resolutions), str(int(resolutions / 3)), str(c_lat), str(c_lon), save_path) return cmd def get_remote_file(lon, lat, save_path, max_retries=3): """Get Copernicus Dem by lon and lat with retry mechanism Args: lon (number): lontitude lat (number): latitude save_path (str): directory to save data max_retries (int): maximum retry attempts """ lon = int(lon) lat = int(lat) cmd = cmd_init(lon, lat, save_path) for attempt in range(max_retries): try: result = subprocess.run(cmd, shell=True, capture_output=True, text=True, timeout=300) if result.returncode == 0: print(f"✓ Successfully downloaded: lon={lon}, lat={lat}") return True else: if attempt < max_retries - 1: print(f"✗ Failed (attempt {attempt+1}/{max_retries}): lon={lon}, lat={lat}, retrying...") import time time.sleep(2) else: print(f"✗ Failed after {max_retries} attempts: lon={lon}, lat={lat}") return False except subprocess.TimeoutExpired: if attempt < max_retries - 1: print(f"⏱ Timeout (attempt {attempt+1}/{max_retries}): lon={lon}, lat={lat}, retrying...") import time time.sleep(2) else: print(f"✗ Timeout after {max_retries} attempts: lon={lon}, lat={lat}") return False except Exception as e: print(f"✗ Exception occurred for lon={lon}, lat={lat}: {e}") return False return False def get_remote_file_batch(west, east, south, north, save_path, num_workers=4): """Get Copernicus Dems by WESN extent with parallel download Args: west (number): western longitude boundary east (number): eastern longitude boundary south (number): southern latitude boundary north (number): northern latitude boundary save_path (str): directory to save data num_workers (int): number of parallel download workers (default: 4) """ from concurrent.futures import ThreadPoolExecutor, as_completed from tqdm import tqdm west = int(west) east = int(east) south = int(south) north = int(north) # 生成所有需要下载的坐标对 coords = [(lon, lat) for lon in range(west, east + 1) for lat in range(south, north + 1)] total_tiles = len(coords) print(f"\n{'='*80}") print(f"并行下载 Copernicus DEM 数据") print(f"{'='*80}") print(f"区域范围: {west}°E - {east}°E, {south}°N - {north}°N") print(f"总瓦片数: {total_tiles}") print(f"并行线程: {num_workers}") print(f"{'='*80}\n") success_count = 0 failed_count = 0 # 使用线程池并行下载 with ThreadPoolExecutor(max_workers=num_workers) as executor: # 提交所有下载任务 future_to_coord = {executor.submit(get_remote_file, lon, lat, save_path): (lon, lat) for lon, lat in coords} # 使用 tqdm 显示进度条 with tqdm(total=total_tiles, desc="下载进度", unit="瓦片") as pbar: for future in as_completed(future_to_coord): coord = future_to_coord[future] try: result = future.result() if result: success_count += 1 else: failed_count += 1 except Exception as e: print(f"\n✗ 下载失败 {coord}: {e}") failed_count += 1 pbar.update(1) # 打印统计信息 print(f"\n{'='*80}") print(f"下载完成统计") print(f"{'='*80}") print(f"成功: {success_count}/{total_tiles} 瓦片") print(f"失败: {failed_count}/{total_tiles} 瓦片") if failed_count > 0: print(f"\n⚠️ 有 {failed_count} 个瓦片下载失败,可能会影响 DEM 完整性") print(f"{'='*80}\n") def download_nasadem(west, south, east, north, save_path, cell_size_deg=0.000277778): """Download NASADEM data using NASADEM library Args: west (float): western longitude boundary south (float): southern latitude boundary east (float): eastern longitude boundary north (float): northern latitude boundary save_path (str): directory to save data cell_size_deg (float): cell size in degrees (default: 0.000277778 ≈ 30m at equator) Returns: str: path to the downloaded DEM file, or None if failed """ if not HAS_NASADEM: print("✗ NASADEM library not available. Cannot download NASADEM.") print(" Install with: pip install nasadem rasters") return None print(f"\n{'='*80}") print(f"下载 NASADEM 数据") print(f"{'='*80}") print(f"区域范围: {west}°E - {east}°E, {south}°N - {north}°N") print(f"分辨率: {cell_size_deg}° (约 {cell_size_deg * 111000:.0f} 米)") print(f"{'='*80}\n") try: from NASADEM.NASADEM import NASADEMConnection from rasters import RasterGrid import time start_time = time.time() # 定义目标区域 geometry = RasterGrid.from_bbox( xmin=west, ymin=south, xmax=east, ymax=north, cell_size=cell_size_deg, crs="EPSG:4326" ) print(f"RasterGrid 维度: {geometry.shape}") print(f"正在从NASA服务器下载高程数据...") print(f"提示: 首次下载需要从NASA服务器获取数据,可能需要几分钟时间...") # 创建NASA DEM连接并获取高程数据 conn = NASADEMConnection() elevation = conn.elevation_m(geometry) elapsed_time = time.time() - start_time print(f"\n✓ 下载完成!") print(f" 高程范围: {float(elevation.min()):.1f} 到 {float(elevation.max()):.1f} 米") print(f" 数据维度: {elevation.shape}") print(f" 用时: {elapsed_time:.1f} 秒") # 保存为 GeoTIFF output_file = os.path.join(save_path, f"NASADEM_{west}_{east}_{south}_{north}.tif") print(f"\n保存到: {output_file}") # 检查elevation对象是否有save方法 if hasattr(elevation, 'save'): elevation.save(output_file) print(f"✓ 文件保存成功: {output_file}") if os.path.exists(output_file): file_size = os.path.getsize(output_file) / (1024 * 1024) print(f" 文件大小: {file_size:.2f} MB") return output_file # 否则使用 rasterio 保存 try: import rasterio from rasterio.crs import CRS from rasterio.transform import from_bounds # 获取数据数组 if hasattr(elevation, 'read'): elev_array = elevation.read() else: elev_array = np.array(elevation) # 计算变换矩阵 height, width = elev_array.shape transform = from_bounds(west, south, east, north, width, height) # 写入文件 with rasterio.open( output_file, 'w', driver='GTiff', height=height, width=width, count=1, dtype=elev_array.dtype, crs=CRS.from_epsg(4326), transform=transform, nodata=-32767, compress='LZW' ) as dst: dst.write(elev_array, 1) print(f"✓ 文件保存成功: {output_file}") if os.path.exists(output_file): file_size = os.path.getsize(output_file) / (1024 * 1024) print(f" 文件大小: {file_size:.2f} MB") return output_file except ImportError: print("⚠️ rasterio 未安装,尝试使用 GDAL...") # 使用 GDAL 创建 GeoTIFF width, height = geometry.width, geometry.height gdal_cmd = f"gdal_create -outsize {width} {height} -of GTiff -co COMPRESS=LZW -a_srs EPSG:4326 -a_ullr {west} {north} {east} {south} {output_file}" result = os.system(gdal_cmd) if result == 0 and os.path.exists(output_file): print(f"✓ 文件创建成功: {output_file}") return output_file else: print(f"✗ 文件创建失败") return None except Exception as e: print(f"\n✗ 下载失败: {e}") import traceback traceback.print_exc() return None def download_srtm(west, south, east, north, save_path, srtm_data_dir=None, cell_size=0.000277778): """Download SRTM data using srtm library (requires pre-downloaded .hgt files) Args: west (float): western longitude boundary south (float): southern latitude boundary east (float): eastern longitude boundary north (float): northern latitude boundary save_path (str): directory to save output DEM srtm_data_dir (str): directory containing SRTM .hgt files (required) cell_size (float): cell size in degrees (default: 0.000277778 ≈ 30m) Returns: str: path to the generated DEM file, or None if failed Note: This function requires SRTM .hgt files to be pre-downloaded. You can download SRTM data from: - https://srtm.csi.cgiar.org/ - https://earthexplorer.usgs.gov/ Python version requirement: - srtm library requires Python >= 3.12 - For Python < 3.12, use process_srtm_hgt_files() instead """ if not HAS_SRTM: print("✗ srtm library not available. Cannot use SRTM.") print(" Install with: pip install srtm (requires Python >= 3.12)") print(f" Current Python version: {sys.version_info.major}.{sys.version_info.minor}") print(" Alternative: Use process_srtm_hgt_files() to manually process .hgt files") return None if srtm_data_dir is None or not os.path.isdir(srtm_data_dir): print("✗ SRTM data directory not provided or does not exist.") print(" Please download SRTM .hgt files and specify the directory.") print(" Download sources:") print(" - https://srtm.csi.cgiar.org/") print(" - https://earthexplorer.usgs.gov/") return None print(f"\n{'='*80}") print(f"使用 SRTM 数据生成 DEM") print(f"{'='*80}") print(f"区域范围: {west}°E - {east}°E, {south}°N - {north}°N") print(f"SRTM数据目录: {srtm_data_dir}") print(f"分辨率: {cell_size}° (约 {cell_size * 111000:.0f} 米)") print(f"{'='*80}\n") try: import time from tqdm import tqdm start_time = time.time() # 创建 SrtmService service = srtm.SrtmService(srtm_data_dir, cache_size=100) print(f"✓ SrtmService 创建成功") # 生成坐标网格 lats = np.arange(south, north + cell_size, cell_size) lons = np.arange(west, east + cell_size, cell_size) print(f"生成网格: {len(lats)} × {len(lons)} = {len(lats) * len(lons)} 个点") # 批量查询高程 coords = [(lat, lon) for lat in lats for lon in lons] print(f"正在查询高程数据...") elevations = service.get_elevations_batch(coords, default=0) # 重塑为数组 elevation_array = np.array(elevations, dtype=np.float32).reshape(len(lats), len(lons)) # 替换0值为NaN(可能是缺失数据) elevation_array[elevation_array == 0] = np.nan elapsed_time = time.time() - start_time print(f"\n✓ 查询完成!") print(f" 有效数据: {np.sum(~np.isnan(elevation_array))} / {elevation_array.size} 点") if np.sum(~np.isnan(elevation_array)) > 0: print(f" 高程范围: {np.nanmin(elevation_array):.1f} 到 {np.nanmax(elevation_array):.1f} 米") print(f" 用时: {elapsed_time:.1f} 秒") # 保存为 GeoTIFF output_file = os.path.join(save_path, f"SRTM_{west}_{east}_{south}_{north}.tif") print(f"\n保存到: {output_file}") try: import rasterio from rasterio.crs import CRS from rasterio.transform import from_bounds height, width = elevation_array.shape transform = from_bounds(west, south, east, north, width, height) # 替换NaN为nodata值 elevation_array[np.isnan(elevation_array)] = -32767 with rasterio.open( output_file, 'w', driver='GTiff', height=height, width=width, count=1, dtype=np.float32, crs=CRS.from_epsg(4326), transform=transform, nodata=-32767, compress='LZW' ) as dst: dst.write(elevation_array.astype(np.float32), 1) print(f"✓ 文件保存成功: {output_file}") if os.path.exists(output_file): file_size = os.path.getsize(output_file) / (1024 * 1024) print(f" 文件大小: {file_size:.2f} MB") return output_file except ImportError: print("⚠️ rasterio 未安装,无法保存文件") return None except Exception as e: print(f"\n✗ 生成失败: {e}") import traceback traceback.print_exc() return None def split_large_region(west, south, east, north, tile_size=1.0): """将大区域分割成小块 Args: west, south, east, north: 区域边界 tile_size: 分块大小（度），默认1°×1° Returns: list: 分块列表 [(west, south, east, north), ...] """ tiles = [] w = west while w < east: e = min(w + tile_size, east) s = south while s < north: n = min(s + tile_size, north) tiles.append((w, s, e, n)) s += tile_size w += tile_size return tiles def download_srtm_cgiar(west, south, east, north, save_path, api_key=None, num_workers=4, dem_type='SRTMGL1', auto_tile=True): """Download DEM data from OpenTopography automatically Args: west (float): western longitude boundary south (float): southern latitude boundary east (float): eastern longitude boundary north (float): northern latitude boundary save_path (str): directory to save output DEM api_key (str): OpenTopography API key (optional, will prompt if not provided) num_workers (int): number of parallel downloads (not used for API, kept for compatibility) dem_type (str): DEM type (default: 'SRTMGL1') - SRTMGL1: SRTM GL1 30m (recommended for SRTM) - SRTMGL3: SRTM GL3 90m - NASADEM: NASADEM 30m - COP30: Copernicus 30m - COP90: Copernicus 90m auto_tile (bool): 自动分块下载大区域 Returns: str: path to the generated DEM file, or None if failed Note: OpenTopography API provides: - SRTMGL1: SRTM GL1 30m (recommended) - SRTMGL3: SRTM GL3 90m - NASADEM: NASADEM Global DEM 30m - COP30: Copernicus 30m - COP90: Copernicus 90m API Documentation: https://portal.opentopography.org/apidocs/ Get API key: https://opentopography.org/myOpenTopo """ # DEM类型信息 dem_info = { 'SRTMGL1': {'resolution': '30m', 'name': 'SRTM GL1'}, 'SRTMGL3': {'resolution': '90m', 'name': 'SRTM GL3'}, 'NASADEM': {'resolution': '30m', 'name': 'NASADEM'}, 'COP30': {'resolution': '30m', 'name': 'Copernicus 30m'}, 'COP90': {'resolution': '90m', 'name': 'Copernicus 90m'} } info = dem_info.get(dem_type, {'resolution': 'unknown', 'name': dem_type}) print(f"\n{'='*80}") print(f"从 OpenTopography 自动下载 DEM 数据") print(f"{'='*80}") print(f"区域范围: {west}°E - {east}°E, {south}°N - {north}°N") print(f"数据源: {info['name']} ({info['resolution']} resolution)") print(f"{'='*80}\n") # 计算区域面积 region_area = (east - west) * (north - south) # 如果区域较大，自动分块 if auto_tile and region_area > 4: print(f"⚠️ 区域较大 ({region_area:.1f} 平方度)，将自动分块下载以提高稳定性") # 计算分块 tiles = split_large_region(west, south, east, north, tile_size=1.0) print(f" 将分成 {len(tiles)} 个 1°×1° 的区块下载\n") # 创建临时目录存放分块文件 temp_dir = os.path.join(save_path, f"temp_tiles_{dem_type}_{west}_{east}_{south}_{north}") if not os.path.exists(temp_dir): os.makedirs(temp_dir) # 下载每个分块 downloaded_files = [] failed_tiles = [] for i, (w, s, e, n) in enumerate(tiles, 1): print(f"\n{'='*80}") print(f"下载分块 {i}/{len(tiles)}: {w}°E-{e}°E, {s}°N-{n}°N") print(f"{'='*80}") tile_file = download_single_tile(w, s, e, n, temp_dir, api_key, dem_type, info) if tile_file: downloaded_files.append(tile_file) print(f"✓ 分块 {i}/{len(tiles)} 下载成功") else: failed_tiles.append((w, s, e, n)) print(f"✗ 分块 {i}/{len(tiles)} 下载失败") # 如果有失败的分块 if failed_tiles: print(f"\n⚠️ 有 {len(failed_tiles)} 个分块下载失败") print("失败的分块:") for w, s, e, n in failed_tiles: print(f" {w}°E-{e}°E, {s}°N-{n}°N") # 合并所有分块 if downloaded_files: print(f"\n{'='*80}") print(f"合并 {len(downloaded_files)} 个分块...") print(f"{'='*80}\n") merged_file = merge_dem_tiles(downloaded_files, save_path, dem_type, west, east, south, north) # 清理临时文件 import shutil try: shutil.rmtree(temp_dir) print("✓ 临时文件已清理") except: pass if merged_file: print(f"\n✓ 最终输出文件: {merged_file}") return merged_file else: print("✗ 合并失败") return None else: print("✗ 所有分块下载失败") return None # 小区域直接下载 else: return download_single_tile(west, south, east, north, save_path, api_key, dem_type, info) def download_single_tile(west, south, east, north, save_path, api_key, dem_type, dem_info): """下载单个DEM块 Args: west, south, east, north: 区域边界 save_path: 保存路径 api_key: OpenTopography API key dem_type: DEM类型 dem_info: DEM信息字典 Returns: str: 下载的文件路径，失败返回None """ try: import time import urllib.request import urllib.parse start_time = time.time() # OpenTopography API endpoint base_url = "https://portal.opentopography.org/API/globaldem" # 如果没有提供API key，提示用户获取 if api_key is None: print("⚠️ 需要OpenTopography API key才能下载SRTM数据") print("\n获取免费API key的步骤:") print(" 1. 访问: https://opentopography.org/myOpenTopo") print(" 2. 注册免费账户") print(" 3. 在 'My Account' 页面请求API key") print(" 4. 使用 --opentopo-api-key 参数提供API key") print("\n或者使用其他DEM数据源:") print(" - Copernicus (推荐): --dem-source copernicus") print(" - NASADEM: --dem-source nasadem") return None # 构建请求参数 params = { 'demtype': dem_type, 'south': south, 'north': north, 'west': west, 'east': east, 'outputFormat': 'GTiff', 'API_Key': api_key } # 构建完整URL url = base_url + '?' + urllib.parse.urlencode(params) # 输出文件名 output_file = os.path.join(save_path, f"{dem_type}_{west}_{east}_{south}_{north}.tif") # 下载文件 - 使用curl以获得更好的稳定性 max_retries = 3 retry_count = 0 download_success = False while retry_count < max_retries and not download_success: try: retry_count += 1 print(f"正在从OpenTopography下载数据... (尝试 {retry_count}/{max_retries})") # 尝试使用 curl 下载（更稳定） try: import subprocess curl_cmd = ['curl', '-L', '-o', output_file, '-s', '--show-error', url] result = subprocess.run(curl_cmd, capture_output=True, text=True, timeout=1800) # 30分钟超时 if result.returncode != 0: raise Exception(f"curl下载失败: {result.stderr}") print(" ✓ curl下载完成") download_success = True except (subprocess.TimeoutExpired, FileNotFoundError): # 如果curl不可用或超时，使用urllib print(" 使用urllib下载...") urllib.request.urlretrieve(url, output_file) download_success = True except urllib.error.HTTPError as e: print(f"\n✗ HTTP错误: {e.code} {e.reason}") if e.code == 401: print(" API key无效或已过期") print(" 请获取新的API key: https://opentopography.org/myOpenTopo") elif e.code == 403: print(" 访问被拒绝，可能是:") print(" - 请求区域过大") print(" - API key配额已用完") print(" - 需要注册获取免费API key") if retry_count < max_retries: print(f" 将在3秒后重试...") time.sleep(3) continue return None except Exception as e: print(f"\n✗ 下载失败: {e}") if retry_count < max_retries: print(f" 将在3秒后重试...") time.sleep(3) continue return None # 检查文件是否有效 if os.path.exists(output_file): file_size = os.path.getsize(output_file) # 如果文件很小，可能是错误消息 if file_size < 1000: with open(output_file, 'r') as f: error_msg = f.read() print(f"\n✗ 下载失败: {error_msg}") os.remove(output_file) return None elapsed_time = time.time() - start_time print(f"\n✓ 下载完成!") print(f" 输出文件: {output_file}") print(f" 文件大小: {file_size / (1024 * 1024):.2f} MB") print(f" 用时: {elapsed_time:.1f} 秒") # 验证文件完整性 print("\n验证文件完整性...") try: import rasterio with rasterio.open(output_file) as src: # 尝试读取多个位置的数据验证 print(f" 文件大小: {src.width} x {src.height} 像素") # 读取多个测试点 test_points = [ (0, 0, min(100, src.height), min(100, src.width)), # 左上角 (max(0, src.height-100), max(0, src.width-100), src.height, src.width), # 右下角 (src.height//2-50, src.width//2-50, src.height//2+50, src.width//2+50), # 中间 ] for i, (row_start, col_start, row_end, col_end) in enumerate(test_points): try: window = ((row_start, row_end), (col_start, col_end)) test_data = src.read(1, window=window) print(f" ✓ 测试点 {i+1}/3 验证成功") except Exception as e: print(f" ✗ 测试点 {i+1}/3 验证失败: {e}") raise print(f" ✓ 所有验证点通过，文件完整") return output_file except Exception as e: print(f" ✗ 文件验证失败: {e}") print(f" 文件可能已损坏，正在删除...") try: os.remove(output_file) except: pass return None else: print("✗ 文件下载失败") return None except Exception as e: print(f"\n✗ 下载失败: {e}") import traceback traceback.print_exc() return None def merge_dem_tiles(tile_files, save_path, dem_type, west, east, south, north): """合并多个DEM分块文件 Args: tile_files: 分块文件列表 save_path: 保存路径 dem_type: DEM类型 west, east, south, north: 最终区域边界 Returns: str: 合并后的文件路径，失败返回None """ if not tile_files: print("✗ 没有文件需要合并") return None try: import rasterio from rasterio.merge import merge import time print(f"开始合并 {len(tile_files)} 个分块文件...") # 读取所有分块 src_files_to_mosaic = [] for tile_file in tile_files: src = rasterio.open(tile_file) src_files_to_mosaic.append(src) # 合并 start_time = time.time() mosaic, out_trans = merge(src_files_to_mosaic) # 获取输出元数据 out_meta = src_files_to_mosaic[0].meta.copy() out_meta.update({ "driver": "GTiff", "height": mosaic.shape[1], "width": mosaic.shape[2], "transform": out_trans, "compress": "lzw" }) # 输出文件名 output_file = os.path.join(save_path, f"{dem_type}_{west}_{east}_{south}_{north}.tif") # 写入合并后的文件 with rasterio.open(output_file, "w", **out_meta) as dest: dest.write(mosaic) # 关闭所有源文件 for src in src_files_to_mosaic: src.close() elapsed_time = time.time() - start_time file_size = os.path.getsize(output_file) / (1024 * 1024) print(f"✓ 合并完成!") print(f" 输出文件: {output_file}") print(f" 文件大小: {file_size:.2f} MB") print(f" 用时: {elapsed_time:.1f} 秒") return output_file except ImportError: print("✗ rasterio 未安装，无法合并文件") print(" 安装方法: pip install rasterio") return None except Exception as e: print(f"✗ 合并失败: {e}") import traceback traceback.print_exc() return None def download_srtm_cgiar_old(west, south, east, north, save_path, num_workers=4): """Download SRTM data from CSI-CGIAR (DEPRECATED - tiles not accessible) This function is kept as backup but CSI-CGIAR tile downloads are not reliable. Use download_srtm_cgiar() with OpenTopography API instead. """ print("警告: CSI-CGIAR瓦片下载不可靠，建议使用OpenTopography API") return None def process_srtm_hgt_files(west, south, east, north, srtm_data_dir, save_path, output_name=None): """Process SRTM .hgt files manually (alternative for Python < 3.12) Args: west (float): western longitude boundary south (float): southern latitude boundary east (float): eastern longitude boundary north (float): northern latitude boundary srtm_data_dir (str): directory containing SRTM .hgt files save_path (str): directory to save output DEM output_name (str): output filename (without extension) Returns: str: path to the generated DEM file, or None if failed Note: SRTM .hgt files are 1°×1° tiles with 1201×1201 pixels (3 arc-second) or 3601×3601 pixels (1 arc-second). This function reads and merges them. """ print(f"\n{'='*80}") print(f"手动处理 SRTM .hgt 文件") print(f"{'='*80}") print(f"区域范围: {west}°E - {east}°E, {south}°N - {north}°N") print(f"SRTM数据目录: {srtm_data_dir}") print(f"{'='*80}\n") try: import rasterio from rasterio.merge import merge from rasterio.crs import CRS import time start_time = time.time() # 查找需要的SRTM文件 hgt_files = [] for lat in range(int(south), int(north) + 1): for lon in range(int(west), int(east) + 1): # 构造文件名 lat_prefix = 'N' if lat >= 0 else 'S' lon_prefix = 'E' if lon >= 0 else 'W' # SRTM文件名格式: N39E116.hgt filename = f"{lat_prefix}{abs(lat):02d}{lon_prefix}{abs(lon):03d}.hgt" filepath = os.path.join(srtm_data_dir, filename) if os.path.exists(filepath): hgt_files.append(filepath) print(f" ✓ 找到文件: {filename}") else: print(f" ⚠ 文件缺失: {filename}") if not hgt_files: print("\n✗ 没有找到任何SRTM文件!") return None print(f"\n找到 {len(hgt_files)} 个SRTM文件") # 使用GDAL合并.hgt文件 merged_file = os.path.join(save_path, f"SRTM_merged_{west}_{east}_{south}_{north}.tif") # 构建gdal_merge命令 gdal_merge_cmd = ['gdal_merge.py', '-o', merged_file, '-of', 'GTiff', '-co', 'COMPRESS=LZW'] gdal_merge_cmd.extend(hgt_files) print(f"\n使用GDAL合并文件...") result = subprocess.run(gdal_merge_cmd, capture_output=True, text=True) if result.returncode != 0 or not os.path.exists(merged_file): print(f"✗ GDAL合并失败: {result.stderr}") return None print(f"✓ 合并完成: {merged_file}") # 裁剪到目标范围 if output_name is None: output_name = os.path.join(save_path, f"SRTM_{west}_{east}_{south}_{north}") output_file = output_name + '.tif' gdal_warp_cmd = [ 'gdalwarp', '-te', str(west), str(south), str(east), str(north), '-of', 'GTiff', '-co', 'COMPRESS=LZW', merged_file, output_file ] print(f"\n裁剪到目标范围...") result = subprocess.run(gdal_warp_cmd, capture_output=True, text=True) if result.returncode != 0 or not os.path.exists(output_file): print(f"✗ GDAL裁剪失败: {result.stderr}") return None elapsed_time = time.time() - start_time print(f"\n✓ SRTM DEM 生成完成!") print(f" 输出文件: {output_file}") if os.path.exists(output_file): file_size = os.path.getsize(output_file) / (1024 * 1024) print(f" 文件大小: {file_size:.2f} MB") print(f" 用时: {elapsed_time:.1f} 秒") # 清理临时文件 if os.path.exists(merged_file) and merged_file != output_file: os.remove(merged_file) print(f" 清理临时文件: {merged_file}") return output_file except ImportError: print("✗ rasterio 未安装,无法处理文件") return None except Exception as e: print(f"\n✗ 处理失败: {e}") import traceback traceback.print_exc() return None def merge_tif_files(input_pattern, output_file): """合并多个 TIF 文件为一个文件 Args: input_pattern (str): 输入 TIF 文件的匹配模式 output_file (str): 输出合并后的 TIF 文件路径 """ # 查找所有匹配的 TIF 文件 tif_files = sorted(glob.glob(input_pattern)) if not tif_files: print(f"未找到匹配的 TIF 文件: {input_pattern}") return False print(f"\n{'='*80}") print(f"合并 DEM 瓦片") print(f"{'='*80}") print(f"找到 {len(tif_files)} 个 TIF 文件待合并") print(f"输出文件: {output_file}") print(f"{'='*80}\n") # 使用 gdal_merge.py 合并文件 try: start_time = __import__('time').time() cmd = ["gdal_merge.py", "-o", output_file, "-co", "COMPRESS=LZW", "-co", "BIGTIFF=YES"] + tif_files print(f"执行合并命令...") # 使用 tqdm 显示进度 result = subprocess.run(cmd, capture_output=True, text=True) elapsed_time = __import__('time').time() - start_time if result.returncode == 0: # 检查输出文件大小 if os.path.exists(output_file): file_size = os.path.getsize(output_file) / (1024 * 1024) # MB print(f"\n✓ 成功合并文件!") print(f" 输出文件: {output_file}") print(f" 文件大小: {file_size:.2f} MB") print(f" 用时: {elapsed_time:.1f} 秒") print(f"{'='*80}\n") return True else: print(f"\n✗ 合并失败: {result.stderr}") print(f"{'='*80}\n") return False except Exception as e: print(f"\n✗ 合并过程中发生错误: {e}") print(f"{'='*80}\n") return False def convert_format(input_file, output_file, format_type="GTiff", options=None): """转换 TIF 文件格式 Args: input_file (str): 输入文件路径 output_file (str): 输出文件路径 format_type (str): 输出格式类型 (默认: GTiff) options (list): gdal_translate 的额外选项 """ if options is None: options = [] try: cmd = ["gdal_translate", "-of", format_type] + options + [input_file, output_file] print(f"执行转换命令: {' '.join(cmd)}") result = subprocess.run(cmd, capture_output=True, text=True) if result.returncode == 0: print(f"成功转换文件到: {output_file}") return True else: print(f"转换失败: {result.stderr}") return False except Exception as e: print(f"转换过程中发生错误: {e}") return False def process_dem_files(save_path, west, east, south, north, merge=True, convert=True, output_format="GTiff", convert_options=None, gamma=False, gamma_byteorder='big', processor='gamma', output_name=None, num_workers=4, srtm_data_dir=None, opentopo_api_key=None, opentopo_dem_type='SRTMGL1', fabdem_dir=None): """处理下载的 DEM 文件：下载、合并和转换 Args: save_path (str): 数据保存目录 west (float): 西经边界 east (float): 东经边界 south (float): 南纬边界 north (float): 北纬边界 merge (bool): 是否合并文件 convert (bool): 是否转换格式 output_format (str): 输出格式 convert_options (list): 转换选项 gamma (bool): 是否转换为 Gamma/ROI_PAC 格式 gamma_byteorder (str): 字节序 ('big' 或 'little') processor (str): 处理器类型 ('gamma' 或 'roi_pac') output_name (str): 输出文件名（不含扩展名） num_workers (int): 并行下载线程数 (default: 4) dem_source (str): DEM数据源 ('copernicus', 'nasadem' 或 'srtm', default: 'copernicus') srtm_data_dir (str): SRTM数据目录 (仅srtm源需要) """ if convert_options is None: convert_options = ["-co", "COMPRESS=LZW", "-co", "TILED=YES"] # ============ 优先尝试本地 FABDEM 瓦片库 ============ if fabdem_dir and os.path.isdir(fabdem_dir): print(f"\n{'='*80}") print(f"优先使用本地 FABDEM 瓦片库: {fabdem_dir}") print(f"{'='*80}") try: from make_local_dem import (find_needed_tiles, extract_tiles_from_zips, tiles_to_gamma_dem) import tempfile, shutil tiles = find_needed_tiles(west, south, east, north) print(f"需要 {len(tiles)} 个 1°×1° FABDEM 瓦片") temp_dir = tempfile.mkdtemp(prefix='fabdem_tiles_') try: tif_files = extract_tiles_from_zips(tiles, fabdem_dir, temp_dir) if tif_files: # 一步完成: VRT → ENVI 二进制 → .dem（仅一次大文件写入） out_name = output_name if output_name else os.path.join(save_path, 'out') dem_file, dem_par_file = tiles_to_gamma_dem( tif_files, out_name, west, south, east, north, gamma_byteorder) if dem_file and dem_par_file: print(f"\n✓ 本地 FABDEM 生成 GAMMA DEM 成功!") print(f" DEM: {dem_file}") print(f" PAR: {dem_par_file}") shutil.rmtree(temp_dir, ignore_errors=True) return else: print("⚠ 本地 FABDEM 转换失败，回退到网络下载") else: print("⚠ 未从本地 FABDEM 提取到瓦片，回退到网络下载") finally: shutil.rmtree(temp_dir, ignore_errors=True) except ImportError as e: print(f"⚠ 无法导入 make_local_dem 模块 ({e})，回退到网络下载") except Exception as e: print(f"⚠ 本地 FABDEM 处理异常 ({e})，回退到网络下载") # ============ 回退: 网络下载 ============ # 检查是否使用预下载的SRTM数据 if srtm_data_dir and os.path.isdir(srtm_data_dir): print(f"\n{'='*80}") print(f"使用预下载的 SRTM 数据") print(f"{'='*80}") if HAS_SRTM: # Python >= 3.12, 使用srtm库 dem_file = download_srtm(west, south, east, north, save_path, srtm_data_dir) else: # Python < 3.12, 使用手动处理方法 print("注意: srtm库需要Python >= 3.12,使用手动处理方法") dem_file = process_srtm_hgt_files(west, south, east, north, srtm_data_dir, save_path, output_name) else: # 自动从OpenTopography下载 print(f"\n{'='*80}") print(f"开始下载 DEM 数据...") print(f"{'='*80}") print("自动从OpenTopography下载DEM数据") dem_file = download_srtm_cgiar(west, south, east, north, save_path, api_key=opentopo_api_key, num_workers=num_workers, dem_type=opentopo_dem_type) if dem_file is None: print("✗ DEM 生成失败,退出") return current_file = dem_file # 生成输出文件名 if output_name is None: output_name = os.path.join(save_path, f"DEM_W{int(west)}_E{int(east)}_S{int(south)}_N{int(north)}") # 转换到 Gamma/ROI_PAC 格式 if gamma: print("\n开始转换到 {} 格式...".format(processor.upper())) convert_to_gamma(current_file, output_name, gamma_byteorder, processor) # 转换到其他格式（非 Gamma/ROI_PAC） if not gamma and convert: print("\n开始转换文件格式...") converted_file = os.path.join(save_path, f"final_DEM_W{int(west)}_E{int(east)}_S{int(south)}_N{int(north)}.tif") convert_format(current_file, converted_file, output_format, convert_options) if __name__ == "__main__": parser = argparse.ArgumentParser( description='Download DEM data from OpenTopography and convert to Gamma/ROI_PAC format.', formatter_class=argparse.RawTextHelpFormatter ) # 与 makedem.py 相同的参数 parser.add_argument('-r', dest='region', help='Research region, west/east/south/north (e.g., 106/110/36/40)') parser.add_argument('-d', dest='dem', help='Raw dem file that used for further processing') parser.add_argument('-s', dest='par', help='SLC parameter file of SAR image used for determining research region') parser.add_argument('-p', dest='processor', help='Interferometry processor. [ gamma or roi_pac ] [default: gamma]') parser.add_argument('-o', dest='out', help='Output name of the generated DEM') parser.add_argument('--byteorder', dest='byteorder', choices=['big', 'little'], help='Byteorder of the generated DEM: big or little. [default: big for gamma and little for roi_pac]') parser.add_argument('--dir', dest='PATH', help='Processing directory for generating DEM. [default: Current directory]') parser.add_argument('--num-workers', dest='num_workers', type=int, default=4, help='Number of parallel download workers (used for batch processing). [default: 4]') parser.add_argument('--srtm-data-dir', dest='srtm_data_dir', default=None, help='Directory containing SRTM .hgt files (optional).\n' 'If not provided, DEM data will be automatically downloaded from OpenTopography.\n' 'Download SRTM data from: https://srtm.csi.cgiar.org/') parser.add_argument('--opentopo-api-key', dest='opentopo_api_key', default='09ad77d34545607fdf5cb182b64ac64e', help='OpenTopography API key for downloading DEM data.\n' 'Get a free API key at: https://opentopography.org/myOpenTopo\n' 'Default key is provided, but you can use your own key.') parser.add_argument('--opentopo-dem-type', dest='opentopo_dem_type', choices=['SRTMGL1', 'SRTMGL3', 'NASADEM', 'COP30', 'COP90'], default='SRTMGL1', help='DEM type for OpenTopography download. [default: SRTMGL1]\n' ' - SRTMGL1: SRTM GL1 30m (recommended)\n' ' - SRTMGL3: SRTM GL3 90m\n' ' - NASADEM: NASADEM 30m\n' ' - COP30: Copernicus 30m\n' ' - COP90: Copernicus 90m') parser.add_argument('--fabdem-dir', dest='fabdem_dir', default=None, help='Directory containing local FABDEM ZIP tiles (e.g., /mnt/ZYD/全球FABDEM).\n' 'If provided, local FABDEM tiles will be used FIRST before network download.\n' 'Falls back to OpenTopography download if local tiles are unavailable.') args = parser.parse_args() # 确定工作目录 if args.PATH: workdir = args.PATH else: workdir = os.getcwd() os.chdir(workdir) # 确定处理器类型 if args.processor: processor = args.processor else: processor = 'gamma' # 确定字节序 if args.byteorder: Byteorder = args.byteorder else: if processor == 'gamma': Byteorder = 'big' else: Byteorder = 'little' # 确定输出名称 if args.out: Name = args.out else: Name = "out" # 处理已有 DEM 文件的情况 if args.dem: dem = args.dem print('Raw dem file is provided: %s.' % dem) # 转换为 TIF 格式（如果不是 TIF） SUFIX = get_sufix(dem) if SUFIX != 'tif': DTIF = dem.replace('.' + SUFIX, '.tif') call_str = f'gdal_translate {dem} -of GTiff {DTIF}' os.system(call_str) DEM = DTIF else: DEM = dem # 转换到 Gamma/ROI_PAC 格式 convert_to_gamma(DEM, Name, Byteorder, processor) BB = Byteorder + ' endian' print('') print('%s %s and %s are generated.' % (BB, Name + '.dem', Name + ('.dem.par' if processor == 'gamma' else '.dem.rsc'))) print('Congratulations! Done!') sys.exit(0) # 处理从 SLC 参数文件确定区域的情况 if args.par: Par = args.par print("SLC_par file is provided: %s" % Par) print("DEM over research region will be downloaded automatically based on %s" % Par) print("DEM data will be downloaded from OpenTopography") call_str = "SLC_corners " + Par + " > corners.txt" os.system(call_str) File = open("corners.txt", "r") InfoLine = File.readlines()[8:10] File.close() MinLat = float(InfoLine[0].split(':')[1].split(' max. ')[0]) MaxLat = float(InfoLine[0].split(':')[2]) MinLon = float(InfoLine[1].split(':')[1].split(' max. ')[0]) MaxLon = float(InfoLine[1].split(':')[2]) north = int(MaxLat) + 2 south = int(MinLat) east = int(MaxLon) + 2 west = int(MinLon) elif args.region: region = args.region west, south, east, north = read_region(region) else: parser.print_usage() sys.exit(os.path.basename(sys.argv[0]) + ': error: research region, raw_demfile and SLC parameter file, at least one is needed.') print('Research region: %s(west) %s(south) %s(east) %s(north)' % (west, south, east, north)) print('>>> Ready to download DEM over research region.') # 处理数据（下载、合并和转换到 Gamma/ROI_PAC 格式） process_dem_files( save_path=workdir, west=west, east=east, south=south, north=north, merge=True, convert=False, gamma=True, gamma_byteorder=Byteorder, processor=processor, output_name=os.path.join(workdir, Name), num_workers=args.num_workers, srtm_data_dir=args.srtm_data_dir, opentopo_api_key=args.opentopo_api_key, opentopo_dem_type=args.opentopo_dem_type, fabdem_dir=args.fabdem_dir )