earthaccess，快速下载DEM和拼接

官方网址

安装

conda install -c conda-forge earthacess

import earthaccess
from osgeo import ogr, gdal, osr, gdalconst
import glob


#earthaccess.login()
# 2811273047@qq.com
# Llq@13486871307
def download_dem(lower_left_lon, lower_left_lat , upper_right_lon, upper_right_lat, save_path):
    earthaccess.login()

    # Search Digital Elevation/Digital Terrain Model
    results = earthaccess.search_data(
        # ASTGTM
        short_name='ASTGTM',#'ASTER Global Digital Elevation Model V003',
        doi='10.5067/ASTER/ASTGTM.003',

        # NASADEM
        #short_name='NASADEM',#'ASTER Global Digital Elevation Model V003',
        #doi='10.5067/MEASURES/NASADEM/NASADEM_HGT.001',

        bounding_box=(lower_left_lon, lower_left_lat , upper_right_lon, upper_right_lat),
        #temporal=("1999-02","2019-03"),
        count=-1
    )

    earthaccess.download(results, save_path)


def RasterMosaic(out_path, dem_list):
    """
    out_path: 输出路径
    dem_lis: dem列表
    """
    #print(dem_list[0])
    inputrasfile1 = gdal.Open(dem_list[0], gdal.GA_ReadOnly) # 第一幅影像
    inputProj1 = inputrasfile1.GetProjection()
    # # 构建 gdal.WarpOptions 对象

    output_proj = osr.SpatialReference()
    output_proj.ImportFromEPSG(4326) #
    warp_options = gdal.WarpOptions(srcSRS=inputProj1, dstSRS=output_proj.ExportToWkt(), format='GTiff',resampleAlg=gdalconst.GRA_Bilinear)
    #warp_options = gdal.WarpOptions(srcSRS=inputProj1, dstSRS=output_proj.ExportToWkt(), format='GTiff',
    # cutlineDSName=shp_path, cropToCutline=True, resampleAlg=gdalconst.GRA_Bilinear)

    # # 批量拼接 DEM 文件
    gdal.Warp(out_path, dem_list, options=warp_options)
    del inputrasfile1

shapefile = "../shp/hgc.shp"
ds = ogr.Open(shapefile)
layer = ds.GetLayer()

extent = layer.GetExtent()
xmin, xmax, ymin, ymax = extent
print(xmin, xmax, ymin, ymax)

# xmin=115
# xmax=120
# ymin=38
# ymax=40
# # Download
download_dem(xmin, ymin, xmax, ymax, 'dem')


# 输入DEM文件列表
import os

# 指定路径
path = "dem"

# 获取路径下的所有文件名
#dem_files = [f for f in os.listdir(path) if f.endswith('_dem.tif')]
dem_files = glob.glob(os.path.join(path,'*_dem.tif'))

print(len(dem_files))
if len(dem_files) > 1:
    # 合并DEM文件
    output_file = "merged_dem.tif"
    RasterMosaic(output_file, dem_files)

    # 转换为int16格式
    output_int16_file = "merged_dem_int16.tif"
    ds = gdal.Open(output_file)
    gdal.Translate(output_int16_file, ds, outputType=gdal.GDT_Int16)

else:
    # 转换为int16格式
    output_int16_file = "merged_dem_int16.tif"
    print(os.path.join(path,dem_files[0]))
    ds = gdal.Open(os.path.join(path,dem_files[0]))
    gdal.Translate(output_int16_file, ds, outputType=gdal.GDT_Int16)

# 关闭数据集
ds = None