xhong/NASA_EarthData_Script
1
0
Fork 0
Code Issues Pull Requests Projects Releases 1 Wiki Activity

fix(DEM_SuPER): 优化ALOSDEM文件检索匹配模式并优化裁剪逻辑.

Browse Source
This commit is contained in:
谢泓 2025-10-12 16:28:28 +08:00
parent 418923f7df
commit fcbfa1bf08
1 changed files with 18 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
DATA_SuPER/DEM_SuPER.py
View File

@ -231,7 +231,7 @@ def process_granule(
dem_file_list = glob.glob(os.path.join(unzip_dir, f"*{name}.hgt")) 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" out_tif_name = f"DEM.{mode_str}.{tile_id}.2000.{name}.tif"
elif mode_str == "ALOSDEM": 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" out_tif_name = f"DEM.{mode_str}.{tile_id}.2011.{name}.tif"
output_file = os.path.join(output_dir, out_tif_name) output_file = os.path.join(output_dir, out_tif_name)
if not os.path.isfile(output_file): if not os.path.isfile(output_file):
@ -250,12 +250,17 @@ def process_granule(
raster.attrs["scale_factor"] = 1 raster.attrs["scale_factor"] = 1
dem_raster_list.append(raster) dem_raster_list.append(raster)
if len(dem_raster_list) >= 1: 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": if name == "slope" or name == "aspect":
dem_mosaiced = mosaic_images(dem_raster_list, nodata=-9999) dem_mosaiced = mosaic_images(dem_raster_list, nodata=-9999)
else: else:
dem_mosaiced = mosaic_images(dem_raster_list, nodata=-32768) 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( dem_mosaiced.rio.to_raster(
output_file, driver="COG", compress="DEFLATE" output_file, driver="COG", compress="DEFLATE"
) )
@ -270,12 +275,13 @@ def process_granule(
def main( def main(
output_root_dir: str, output_root_dir: str,
region: list, region_file: str,
asset_name: list, asset_name: list,
tile_id: str, tile_id: str,
mode_str: str = "NASADEM", mode_str: str = "NASADEM",
): ):
bbox = tuple(list(region.total_bounds)) region = gpd.read_file(region_file)
bbox = tuple(list(region.total_bounds)) # 最大外接矩形
results_urls = [] results_urls = []
mode_str = mode_str.upper() mode_str = mode_str.upper()
output_root_dir = os.path.join(output_root_dir, mode_str) 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) 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: with open(results_urls_file, "w") as f:
json.dump(results_urls, 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: if len(results_urls) == 0:
return return
@ -348,13 +357,13 @@ def main(
if __name__ == "__main__": if __name__ == "__main__":
earthaccess.login(strategy="netrc", persist=True) 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 = "49REL"
tile_id = "Wuhan" 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"] # asset_name = ["NASADEM_HGT", "NASADEM_SC"]
# mode_str = "NASADEM" # mode_str = "NASADEM"
asset_name = ["ALOS_PSR_RTC_HIGH"] asset_name = ["ALOS_PSR_RTC_HIGH"]
mode_str = "ALOSDEM" mode_str = "ALOSDEM"
output_root_dir = ".\\data\\DEM" 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)