栅格数据的批量镶嵌(附Python脚本)

博客小序:在数据处理的过程中,会遇到需要大量镶嵌的情况,当数据较多时手动镶嵌较为麻烦,自己最近对分省的DEM数据进行镶嵌,由于利用python进行镶嵌较为方便,特撰此博文以记之。
参考博客:
https://blog.csdn.net/qq_15642411/article/details/79187787
https://blog.csdn.net/XBR_2014/article/details/85255412

1.脚本处理情况说明

本实例中,处理的数据是分省的DEM数据,每个省由若干数量不同的tif文件,本脚本主要的功能只有一个即: 对分省的DEM数据进行分省镶嵌

2.脚本代码

#添加一个计时器
import time
start = time.time()

import os, shutil, glob
from osgeo import gdal

# 定义一个镶嵌的函数,传入的参数是需要镶嵌的数据的列表,以及输出路径
def mosaic(data_list, out_path):

    # 读取其中一个栅格数据来确定镶嵌图像的一些属性
    o_ds = gdal.Open(data_list[0])
    # 投影
    Projection = o_ds.GetProjection()
    # 波段数据类型
    o_ds_array = o_ds.ReadAsArray()

    if 'int8' in o_ds_array.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in o_ds_array.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    # 像元大小
    transform = o_ds.GetGeoTransform()
    pixelWidth = transform[1]
    pixelHeight = transform[5]

    del o_ds

    minx_list = []
    maxX_list = []
    minY_list = []
    maxY_list = []

    # 对于每一个需要镶嵌的数据,读取它的角点坐标
    for data in data_list:

        # 读取数据
        ds = gdal.Open(data)
        rows = ds.RasterYSize
        cols = ds.RasterXSize

        # 获取角点坐标
        transform = ds.GetGeoTransform()
        minX = transform[0]
        maxY = transform[3]
        pixelWidth = transform[1]
        pixelHeight = transform[5]  # 注意pixelHeight是负值
        maxX = minX + (cols * pixelWidth)
        minY = maxY + (rows * pixelHeight)

        minx_list.append(minX)
        maxX_list.append(maxX)
        minY_list.append(minY)
        maxY_list.append(maxY)

        del ds

    # 获取输出图像坐标
    minX = min(minx_list)
    maxX = max(maxX_list)
    minY = min(minY_list)
    maxY = max(maxY_list)

    # 获取输出图像的行与列
    cols = int((maxX - minX) / pixelWidth)
    rows = int((maxY - minY) / abs(pixelHeight))# 注意pixelHeight是负值

    # 计算每个图像的偏移值
    xOffset_list = []
    yOffset_list = []
    i = 0

    for data in data_list:
        xOffset = int((minx_list[i] - minX) / pixelWidth)
        yOffset = int((maxY_list[i] - maxY) / pixelHeight)
        xOffset_list.append(xOffset)
        yOffset_list.append(yOffset)
        i += 1

    # 创建一个输出图像
    driver = gdal.GetDriverByName("GTiff")
    dsOut = driver.Create(out_path + ".tif", cols, rows, 1, datatype) 
    bandOut = dsOut.GetRasterBand(1)

    i = 0
    #将原始图像写入新创建的图像
    for data in data_list:
        # 读取数据
        ds = gdal.Open(data)
        data_band = ds.GetRasterBand(1)
        data_rows = ds.RasterYSize
        data_cols = ds.RasterXSize

        data = data_band.ReadAsArray(0, 0, data_cols, data_rows)
        bandOut.WriteArray(data, xOffset_list[i], yOffset_list[i])

        del ds
        i += 1

    # 设置输出图像的几何信息和投影信息
    geotransform = [minX, pixelWidth, 0, maxY, 0, pixelHeight]
    dsOut.SetGeoTransform(geotransform)
    dsOut.SetProjection(Projection)

    del dsOut

def main():

    input_folder = "D:\\cnblogs\\data\\china"
    file_list = glob.glob(input_folder + "\\*")

    out_file = "D:\\cnblogs\\data\\china_moasic"
    if os.path.exists(out_file):
        shutil.rmtree(out_file)
        os.mkdir(out_file)
    else:
        os.mkdir(out_file)

    for file in file_list:
        basename = os.path.basename(file)
        out_path = out_file + "\\" + basename

        data_list = glob.glob(file + "\\" + "*.tif")
        print(data_list)

        try:
            mosaic(data_list, out_path)
            print(file + "镶嵌结束")
        except:
            bad_list.append(file)
            print(file + "数据超过4G或其他原因导致无法镶嵌")

bad_list = []
main()
print("无法镶嵌的文件包括如下")
print (bad_list)

end = time.time()
print ("程序运行时间{:.2f}分钟".format((end-start)/60.0))

3.问题总结

1.多幅栅格数据镶嵌时,由于镶嵌后的栅格数据大小超过4G,python的会报错。
2.此脚本只考虑了栅格数据只有一个波段的情形,多波段栅格数据的镶嵌未来有机会遇到再补充。
3.由于新生成的镶嵌数据是规则的矩形,但是用于镶嵌的数据不一定能刚好完美的覆盖新的图层,导致没有数据用来镶嵌的部分的值0,此问题较为重要,需要注意。


本文作者:DQTDQT
限于作者水平有限,如文中存在任何错误,欢迎不吝指正、交流。

联系方式:
QQ:1426097423
e-mail:duanquntaoyx@163.com

本文版权归作者和博客园共有,欢迎转载、交流,但未经作者同意必须保留此段声明,且在文章页面明显位置给出原文链接,如果觉得本文对您有益,欢迎点赞、探讨。

转载于:https://www.cnblogs.com/xsman/p/11439894.html

  • 0
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
### 回答1: 批量裁剪栅格数据是指对多个栅格数据进行统一的裁剪操作,可以使用Python编写代码来实现该功能。 首先,需要导入必要的库,如`osgeo`库用于读取和处理栅格数据,`os`库用于文件操作。 ```python from osgeo import gdal, gdalnumeric import os ``` 接下来,可以定义一个函数来实现裁剪操作。 ```python def crop_raster(input_path, output_path, extent): # 打开输入栅格数据 ds = gdal.Open(input_path) if ds is None: print("打开栅格数据失败!") return # 读取输入栅格数据的信息 geotransform = ds.GetGeoTransform() projection = ds.GetProjection() x_size = ds.RasterXSize y_size = ds.RasterYSize # 计算裁剪后的输出栅格数据大小 ulx, uly = extent[0], extent[3] lrx, lry = extent[2], extent[1] px_start = int((ulx - geotransform[0]) / geotransform[1]) px_end = int((lrx - geotransform[0]) / geotransform[1]) py_start = int((uly - geotransform[3]) / geotransform[5]) py_end = int((lry - geotransform[3]) / geotransform[5]) px_count = px_end - px_start py_count = py_end - py_start # 创建输出栅格数据 driver = gdal.GetDriverByName("GTiff") output_ds = driver.Create(output_path, px_count, py_count, 1, gdal.GDT_Float32) # 设置输出栅格数据的信息 output_ds.SetProjection(projection) output_ds.SetGeoTransform((ulx, geotransform[1], 0, uly, 0, geotransform[5])) # 读取输入栅格数据的像素值,并写入输出栅格数据 input_data = ds.GetRasterBand(1).ReadAsArray(px_start, py_start, px_count, py_count) output_band = output_ds.GetRasterBand(1) output_band.WriteArray(input_data) # 关闭栅格数据集 ds = None output_ds = None ``` 在主程序中,可以通过遍历需要裁剪的栅格数据路径列表,调用`crop_raster`函数进行裁剪。 ```python if __name__ == "__main__": # 定义输入栅格数据路径列表 input_paths = ["input_1.tif", "input_2.tif", "input_3.tif"] # 定义输出栅格数据路径列表 output_paths = ["output_1.tif", "output_2.tif", "output_3.tif"] # 定义裁剪范围,格式为 [xmin, ymin, xmax, ymax] extent = [100, 50, 200, 150] # 遍历所有输入栅格数据 for i in range(len(input_paths)): input_path = input_paths[i] output_path = output_paths[i] # 调用裁剪函数进行裁剪 crop_raster(input_path, output_path, extent) print("裁剪完成!") ``` 以上代码片段可以实现对多个栅格数据进行批量裁剪,并将裁剪结果保存到指定的输出路径中。裁剪范围可以根据实际需求进行调整。 ### 回答2: 批量裁剪栅格数据是指对多个栅格数据进行裁剪操作,并将结果保存为新的栅格数据。以下是一个用 Python 编写的简单批量裁剪栅格数据的代码示例: ```python import arcpy # 设置工作空间 arcpy.env.workspace = "C:/RasterData" # 设置裁剪区域 clip_features = "C:/ClipData/clip.shp" # 获取所有栅格数据的文件名 raster_files = arcpy.ListRasters() # 遍历每个栅格数据进行裁剪 for raster_file in raster_files: # 构建裁剪后的栅格数据的输出路径 out_raster = "C:/OutputData/" + raster_file.split('.')[0] + "_clipped.tif" # 裁剪栅格数据 arcpy.Clip_management(raster_file, "#", out_raster, clip_features, "#", "ClippingGeometry") print("成功裁剪栅格数据:" + raster_file) print("批量裁剪栅格数据完成!") ``` 以上代码中,首先通过`arcpy.env.workspace`设置工作空间为包含待裁剪栅格数据的文件夹路径。然后使用`arcpy.ListRasters()`获取所有栅格数据的文件名。接下来通过遍历每个栅格数据文件,使用`arcpy.Clip_management()`函数进行裁剪操作,并将结果保存到指定文件夹。在裁剪完成后,程序会打印出成功裁剪的栅格数据文件名,并最终显示“批量裁剪栅格数据完成!”的提示。 ### 回答3: 批量裁剪栅格数据是一种对多个栅格数据进行裁剪操作的需求,可以使用Python实现。下面是一个大致的代码示例: ```python import os from osgeo import gdal from osgeo import ogr from osgeo import gdalconst # 设置要裁剪的区域范围 extent = [xmin, xmax, ymin, ymax] # 根据实际需要设置裁剪范围坐标 # 设置输入栅格数据的文件夹路径 input_folder = '/path/to/input/folder/' # 设置输出栅格数据的文件夹路径 output_folder = '/path/to/output/folder/' # 获取输入文件夹中的所有栅格数据文件 file_list = os.listdir(input_folder) tif_list = [file for file in file_list if file.endswith('.tif')] # 遍历每个栅格数据文件 for tif_file in tif_list: # 打开栅格数据 input_path = os.path.join(input_folder, tif_file) input_ds = gdal.Open(input_path, gdalconst.GA_ReadOnly) # 获取栅格数据的投影信息、地理变换等 prj = input_ds.GetProjection() geotransform = input_ds.GetGeoTransform() # 根据裁剪范围创建输出栅格数据的文件名和路径 output_path = os.path.join(output_folder, tif_file) # 创建输出栅格数据 output_ds = gdal.Warp(output_path, input_ds, format='GTiff', cutlineDSName=extent, cropToCutline=True) # 关闭栅格数据文件 input_ds = None output_ds = None ``` 以上代码使用了`gdal`库来进行栅格数据的读取和裁剪操作。首先,设置要裁剪的区域范围(extent)和输入、输出文件夹的路径。然后,遍历输入文件夹中的每个栅格数据文件,在每次循环中,打开栅格数据文件,获取其投影信息和地理变换参数,并根据裁剪范围创建输出栅格数据的文件名和路径。最后,使用`gdal.Warp`函数对栅格数据进行裁剪,并关闭输入和输出栅格数据文件。 注意:以上代码仅为大致示例,并未完全测试,实际使用中可能需要根据具体数据格式和需求进行适当调整和优化。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值