diff --git a/DATA_SuPER/DEM_SuPER.py b/DATA_SuPER/DEM_SuPER.py index 32a2391..0fe6abf 100644 --- a/DATA_SuPER/DEM_SuPER.py +++ b/DATA_SuPER/DEM_SuPER.py @@ -139,7 +139,7 @@ def process_granule( name: str, roi: list, clip=True, - tile_id: str = None, + tile_id: str = "", ) -> bool: """ 读取解压并重命名处理后的指定类型 NASADEM 数据并进行预处理, 包括读取, 裁剪, 镶嵌, 并对坡度坡向进行缩放 @@ -215,12 +215,10 @@ def main(region: list, asset_name: list, tile_id: str): os.makedirs(unzip_dir, exist_ok=True) os.makedirs(output_dir, exist_ok=True) results_urls_file = f"{output_root_dir}\\NASADEM_{tile_id}_results_urls.json" - if not os.path.isfile(results_urls_file): - results_urls = earthdata_search(asset_name, roi=bbox) - with open(results_urls_file, "w") as f: - json.dump(results_urls, f) - else: - results_urls = json.load(open(results_urls_file)) + # 默认覆盖上一次检索记录 + results_urls = earthdata_search(asset_name, roi=bbox) + with open(results_urls_file, "w") as f: + json.dump(results_urls, f) # 构造待解压的文件列表 zip_file_list = [os.path.join(download_dir, os.path.basename(result[0])) for result in results_urls] diff --git a/README.md b/README.md index 4fbf0b7..599b3dc 100644 --- a/README.md +++ b/README.md @@ -59,7 +59,7 @@ envs_dirs: - D:\program\miniforge3\envs - 其他路径地址(可选，创建虚拟环境时将会按照顺序查找) channels: - - defaults + - conda-forge show_channel_urls: true channel_alias: https://mirrors.tuna.tsinghua.edu.cn/anaconda default_channels: @@ -117,6 +117,13 @@ mamba env create -f setup/environment.yml mamba activate lpdaac ``` +- 导出当前虚拟环境 + +```sh +mamba env export > setup/lpdaac.yml +``` + + ## 3 设计思路 ### 3.1 数据组织 diff --git a/utils/common_utils.py b/utils/common_utils.py index 803dd18..c5ad324 100644 --- a/utils/common_utils.py +++ b/utils/common_utils.py @@ -5,7 +5,7 @@ ------------------------------------------------------------------------------- Authors: Hong Xie -Last Updated: 2025-07-07 +Last Updated: 2025-08-13 =============================================================================== """ @@ -506,9 +506,9 @@ def reproject_image( target_image.shape[1] == image.shape[1] and target_image.shape[2] == image.shape[2] ): - # 若判断为降尺度/等尺度, 则直接使用 cubic_spline 重采样投影到目标影像 + # 若判断为降尺度/等尺度, 则直接使用 cubic 重采样投影到目标影像 image_reprojed = image.rio.reproject_match( - target_image, resampling=Resampling.cubic_spline + target_image, resampling=Resampling.cubic ) else: # print("target_image shape is not match with image shape", image.shape, "to", target_image.shape) @@ -520,7 +520,7 @@ def reproject_image( # 使用 target_crs 进行重投影 reproject_kwargs = { "dst_crs": target_crs, - "resampling": Resampling.cubic_spline, + "resampling": Resampling.cubic, } if target_shape is not None: reproject_kwargs["shape"] = target_shape @@ -563,7 +563,7 @@ def mosaic_images( tif_list, nodata=nodata, method=method, - resampling=Resampling.cubic_spline, + resampling=Resampling.cubic, ) # 将结果重新构建为 xarray 数据集 # 单张SAR影像直接读取 transform: 233400.0 30.0 0.0 3463020.0 0.0 -30.0 diff --git a/utils/sr2rgb.py b/utils/sr2rgb.py new file mode 100644 index 0000000..6934fa2 --- /dev/null +++ b/utils/sr2rgb.py @@ -0,0 +1,97 @@ +""" +将 COG 格式的 Red, Green, Blue 单波段地表反射率图像合成 RGB 图像 + +1. 对比度线性拉伸至 0-255; +2. 合并波段; +3. 保存为 uint8 格式 RGB 图像; +""" + +import os +import rasterio as rio +import numpy as np + +def sr2rgb(red_path, green_path, blue_path, output_path): + """ + 将红、绿、蓝三个单波段地表反射率图像合成为RGB图像 + + 参数: + red_path (str): 红色波段文件路径 + green_path (str): 绿色波段文件路径 + blue_path (str): 蓝色波段文件路径 + output_path (str): 输出RGB图像路径 + """ + # 检查文件是否存在 + for path in [red_path, green_path, blue_path]: + if not os.path.exists(path): + raise FileNotFoundError(f"文件不存在: {path}") + + print(f"正在处理文件:") + print(f" 红波段: {red_path}") + print(f" 绿波段: {green_path}") + print(f" 蓝波段: {blue_path}") + + # 读取三个波段数据 + with rio.open(red_path) as red_src: + red_band = red_src.read(1) + profile = red_src.profile + print(f"红波段形状: {red_band.shape}, 数据类型: {red_band.dtype}") + + with rio.open(green_path) as green_src: + green_band = green_src.read(1) + print(f"绿波段形状: {green_band.shape}, 数据类型: {green_band.dtype}") + + with rio.open(blue_path) as blue_src: + blue_band = blue_src.read(1) + print(f"蓝波段形状: {blue_band.shape}, 数据类型: {blue_band.dtype}") + + # 线性拉伸到0-255范围 + def stretch_band(band): + # 处理NaN值 + band_no_nan = np.where(np.isnan(band), 0, band) + band_min = np.min(band_no_nan) + band_max = np.max(band_no_nan) + + # 避免除零错误 + if band_max == band_min: + stretched = np.zeros_like(band_no_nan, dtype=np.uint8) + else: + stretched = ((band_no_nan - band_min) / (band_max - band_min) * 255).astype(np.uint8) + return stretched + + red_stretched = stretch_band(red_band) + green_stretched = stretch_band(green_band) + blue_stretched = stretch_band(blue_band) + + # 合并三个波段为RGB图像 + rgb_array = np.dstack((red_stretched, green_stretched, blue_stretched)) + + # 更新元数据 + profile.update( + dtype=rio.uint8, + count=3, + photometric='RGB', + nodata=0 # 设置nodata为0，因为uint8的范围是0-255 + ) + + # 保存RGB图像 + with rio.open(output_path, 'w', **profile) as dst: + for i in range(3): + dst.write(rgb_array[:, :, i], i + 1) + + print(f"RGB图像已保存到: {output_path}") + + +if __name__ == "__main__": + # tif_dir = "D:\\NASA_EarthData_Script\\data\\HLS\\2024\\2024012" + # red_path = os.path.join(tif_dir, "HLS.S30.T49RGP.2024012T031101.v2.0.RED.subset.tif") + # green_path = os.path.join(tif_dir, "HLS.S30.T49RGP.2024012T031101.v2.0.GREEN.subset.tif") + # blue_path = os.path.join(tif_dir, "HLS.S30.T49RGP.2024012T031101.v2.0.BLUE.subset.tif") + # output_path = os.path.join(tif_dir, "HLS.S30.T49RGP.2024012T031101.v2.0.RGB.tif") + + tif_dir = "D:\\NASA_EarthData_Script\\data\\HLS\\2025\\2025011" + red_path = os.path.join(tif_dir, "HLS.S30.T49RGP.2025011T031009.v2.0.RED.subset.tif") + green_path = os.path.join(tif_dir, "HLS.S30.T49RGP.2025011T031009.v2.0.GREEN.subset.tif") + blue_path = os.path.join(tif_dir, "HLS.S30.T49RGP.2025011T031009.v2.0.BLUE.subset.tif") + output_path = os.path.join(tif_dir, "HLS.S30.T49RGP.2025011T031009.v2.0.RGB.tif") + + sr2rgb(red_path, green_path, blue_path, output_path) \ No newline at end of file