Python实现栅格数据拼接、裁剪等操作

从地理空间数据云(http://www.gscloud.cn/)上下载了GDEMV2—30米分辨率的数据,检索按照行政区来的。

下载完的数据解压缩以后,如下:

1.首先把所有的文件夹下的dem数据路径读入到TXT文件中,代码:

此部分参考:https://blog.csdn.net/weixin_42236288/article/details/81810761

import os

# 查询符合条件的文件
f = open(r'D:\ProfessionalProfile\DEMdata\StatisticsDEMfile.txt','w')
j,k=0,0
def file(root, ext):

    for i in os.listdir(root):  # os.listdir(root) 读取root目录下的文件和目录

        # os.path.splitext()将文件名和扩展名分开
        ext = os.path.splitext(root + '\\' + i)  # os.path.splitext()读取文件扩展名

        if ext[1]:  # 有扩展名时 为文件
            global j
            if i.count(ex) > 0:  # 包含特定字符
                pathDEM = root +'\\'+ i
                f.write(pathDEM + '\n')
                f.flush()
                # print(pathDEM)  # 输出
                j += 1
        else:
            global k
            if i.count(ex) > 0:  # 包含特定字符
                pathDEM2 = root + r'\\' + i
                # print( pathDEM2)  # 输出
                k += 1
            root1 = root + '\\' + i  # 拼接目录

            file(root1, ext)  # 递归调用

    print('文件数量为:', j)
    # print('文件夹数量为:', k)
# f.close()

# root = os.getcwd()  # 生成目录 root
# 包含的特定字符
ex = 'dem'
# 调用函数
file('D:\ProfessionalProfile\DEMdata\GDEMV2_30m_shandong','.tif')

运行后得到文件:

2.对该TXT文件里的影像进行拼接

参考:https://blog.csdn.net/qq_43177210/article/details/108402353?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.control&dist_request_id=2a4eb29a-3fae-44ad-8384-ce28a6a41f0a&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.control

# 如果图片大于4G,需要在输入图片路径前加入 -co BIGTIFF=YES
# subprocess 模块首先推荐使用的是它的 run 方法,更高级的用法可以直接使用 Popen 接口。
import subprocess

# 路径
create_slope = r'''E:\r_python\python.exe E:\r_python\Scripts/gdal_merge.py -of GTiff -o '''
list = []
# 打开TXT文件路径
f = open(r'D:\ProfessionalProfile\DEMdata\StatisticsDEMfile.txt', 'r')
# 获取TXT文件中所有DEMTIF文件路径
# readlines()函数:1、一次性读取整个文件。2、自动将文件内容分析成一个行的列表。
lines = f.readlines()
# strip()函数:用于移除字符串头尾指定的字符(默认为空格或换行符)或字符序列。# 注意:该方法只能删除开头或是结尾的字符,不能删除中间部分的字符。
fileDir = [x.strip() for x in lines]
# 遍历所有TIF文件
list = fileDir
# print(list)

# 生成文件路径
tarDir = r'D:\ProfessionalProfile\DEMdata'
#这句没懂
filename = ' '.join(list)
# l3 = len(filename)
# print(filename)
# subprocess.call执行指定的命令,返回命令执行状态,其功能类似于os.system(cmd)。详见:https://www.cnblogs.com/zhou2019/p/10582716.html
# subprocess.call(args, *, stdin=None, stdout=None, stderr=None, shell=False, timeout=None)
subprocess.call(create_slope + r'D:\ProfessionalProfile\DEMdata\demMosaic1224.tif' + ' ' + '-co COMPRESS=LZW ' + filename)

运行结果如下(ENVI中打开):

此外,还学习到了很多,如:

subprocess模块的信息可参考:https://www.cnblogs.com/zhou2019/p/10582716.html

https://blog.csdn.net/m0_38051293/article/details/103218755?utm_medium=distribute.pc_relevant.none-task-blog-baidujs_title-7&spm=1001.2101.3001.4242

https://blog.csdn.net/weixin_40450867/article/details/104191939

https://blog.csdn.net/qq_36178899/article/details/83989031?utm_medium=distribute.pc_relevant_download.none-task-blog-baidujs-1.nonecase&depth_1-utm_source=distribute.pc_relevant_download.none-task-blog-baidujs-1.nonecase

3.裁剪

https://blog.csdn.net/weixin_40501429/article/details/114140680

  • 4
    点赞
  • 29
    收藏
    觉得还不错? 一键收藏
  • 5
    评论
实现按照矢量数据的属性裁剪栅格数据,可以使用Python中的GDAL库和OGR库。具体步骤如下: 1.读取矢量数据和栅格数据,可以使用OGR库和GDAL库中的相关函数进行读取和处理。 2.获取矢量数据的属性信息,可以使用OGR库中的GetField函数获取属性值。 3.遍历栅格数据中的每一个像素,判断该像素是否在矢量数据范围内。 4.如果该像素在矢量数据范围内,则根据矢量数据的属性值进行裁剪操作。 5.将裁剪后的栅格数据保存为新的栅格数据。 下面是一个简单的示例代码: ```python from osgeo import gdal, ogr # 读取矢量数据 vector = ogr.Open('vector.shp') layer = vector.GetLayer() feature = layer.GetFeature(0) field_value = feature.GetField('field_name') # 读取栅格数据 raster = gdal.Open('raster.tif') band = raster.GetRasterBand(1) cols = raster.RasterXSize rows = raster.RasterYSize # 创建输出栅格数据 driver = gdal.GetDriverByName('GTiff') out_raster = driver.Create('out_raster.tif', cols, rows, 1, band.DataType) out_band = out_raster.GetRasterBand(1) # 遍历栅格数据 for i in range(rows): for j in range(cols): # 获取像素值 pixel_value = band.ReadAsArray(j, i, 1, 1)[0][0] # 判断像素是否在矢量数据范围内 if is_in_vector(layer, j, i): # 根据属性值进行裁剪操作 if pixel_value == field_value: out_band.WriteArray(pixel_value, j, i) # 保存输出栅格数据 out_band.FlushCache() out_band.ComputeStatistics(False) out_raster.SetProjection(raster.GetProjection()) out_raster.SetGeoTransform(raster.GetGeoTransform()) out_raster = None ``` 其中,is_in_vector函数可以根据像素的坐标判断该像素是否在矢量数据范围内。具体实现可以参考OGR库中的SpatialReference类和Geometry类。

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值