diff --git a/DATA_SuPER/DEM_SuPER.py b/DATA_SuPER/DEM_SuPER.py index 36a6675..53a643c 100644 --- a/DATA_SuPER/DEM_SuPER.py +++ b/DATA_SuPER/DEM_SuPER.py @@ -231,7 +231,7 @@ def process_granule( dem_file_list = glob.glob(os.path.join(unzip_dir, f"*{name}.hgt")) out_tif_name = f"DEM.{mode_str}.{tile_id}.2000.{name}.tif" elif mode_str == "ALOSDEM": - dem_file_list = glob.glob(os.path.join(unzip_dir, f"*{name}.tif")) + dem_file_list = glob.glob(os.path.join(unzip_dir, "*.dem.tif")) out_tif_name = f"DEM.{mode_str}.{tile_id}.2011.{name}.tif" output_file = os.path.join(output_dir, out_tif_name) if not os.path.isfile(output_file): @@ -250,12 +250,17 @@ def process_granule( raster.attrs["scale_factor"] = 1 dem_raster_list.append(raster) if len(dem_raster_list) >= 1: + # 先逐一裁剪再镶嵌合成(仅在提供ROI且要求裁剪时执行) + if roi is not None and clip: + clipped_list = [] + for raster in dem_raster_list: + clipped = clip_image(raster, roi, clip_by_box=True) + clipped_list.append(clipped) + # 镶嵌合成 if name == "slope" or name == "aspect": dem_mosaiced = mosaic_images(dem_raster_list, nodata=-9999) else: dem_mosaiced = mosaic_images(dem_raster_list, nodata=-32768) - if roi is not None and clip: - dem_mosaiced = clip_image(dem_mosaiced, roi, clip_by_box=True) dem_mosaiced.rio.to_raster( output_file, driver="COG", compress="DEFLATE" ) @@ -270,12 +275,13 @@ def process_granule( def main( output_root_dir: str, - region: list, + region_file: str, asset_name: list, tile_id: str, mode_str: str = "NASADEM", ): - bbox = tuple(list(region.total_bounds)) + region = gpd.read_file(region_file) + bbox = tuple(list(region.total_bounds)) # 最大外接矩形 results_urls = [] mode_str = mode_str.upper() output_root_dir = os.path.join(output_root_dir, mode_str) @@ -297,9 +303,12 @@ def main( # 默认覆盖上一次检索记录 results_urls = earthdata_search(asset_name, roi=bbox) + if mode_str == "ALOSDEM": + # ALOSDEM 数据仅保存含 "FBS" 字符串的 URL + results_urls = [url for url in results_urls if "FBS" in url[0]] with open(results_urls_file, "w") as f: json.dump(results_urls, f) - logging.info(f"Found {len(results_urls)} {mode_str} granules.") + logging.info(f"Found {len(results_urls)} {mode_str} granules after filtering.") if len(results_urls) == 0: return @@ -348,13 +357,13 @@ def main( if __name__ == "__main__": earthaccess.login(strategy="netrc", persist=True) - # region = gpd.read_file("./data/vectors/wuling_guanqu_polygon.geojson") + # region_file = "./data/vectors/wuling_guanqu_polygon.geojson" # tile_id = "49REL" tile_id = "Wuhan" - region = gpd.read_file(f"./data/vectors/{tile_id}.geojson") + region_file = f"./data/vectors/{tile_id}.geojson" # asset_name = ["NASADEM_HGT", "NASADEM_SC"] # mode_str = "NASADEM" asset_name = ["ALOS_PSR_RTC_HIGH"] mode_str = "ALOSDEM" output_root_dir = ".\\data\\DEM" - main(output_root_dir, region, asset_name, tile_id, mode_str) + main(output_root_dir, region_file, asset_name, tile_id, mode_str)