【Python&RS】基于Python批量镶嵌拼接遥感影像/栅格数据

        我之前分享过【Python&RS】基于GDAL镶嵌拼接遥感影像,但是没有加入批量处理的代码。最近正好有这个需求,所以就对原来的代码进行了优化加入了批量拼接的代码。现在只需输入一个文件夹即可将其中的影像全部镶嵌起来。

 一、导入GDAL库

from osgeo import gdal

二、查看影像信息

        为了凑字数的,可以查看影像的投影、宽度、高度、波段数等信息。不过需要注意的是在ENVI中没有投影坐标系,只有地理坐标系是做不了镶嵌拼接的。

        这个代码我还不太清楚能不能不要投影坐标系进行拼接,你们可以自己试试。但最好还是用包含投影坐标系的影像进行拼接。所以在拼接之前就可以用这段代码先看一看。

def Get_data(filepath):
    ds = gdal.Open(filepath)  # 打开数据集dataset
    ds_width = ds.RasterXSize  # 获取数据宽度
    ds_height = ds.RasterYSize  # 获取数据高度
    ds_bands = ds.RasterCount  # 获取波段数
    ds_geo = ds.GetGeoTransform()  # 获取仿射地理变换参数
    ds_prj = ds.GetProjection()  # 获取投影信息
    print("影像的宽度为:" + str(ds_width))
    print("影像的高度为:" + str(ds_height))
    print("仿射地理变换参数为:" + str(ds_geo))
    print("投影坐标系为:" + str(ds_prj))
    # data = ds.ReadAsArray(0, 0, ds_width, ds_height)  # 以数组的形式读取整个数据集

三、镶嵌模块

        代码中的srcSRS,dstSRS分别是输入投影和输出投影,这里用一样的就行了。因为我们做的是镶嵌操作,肯定是不用动原始坐标系的。其他的参数都在代码中表明了,这里就不介绍了,如果大家有什么问题,可以留言交流。

def Mosaic_GDAL(path_image):
    """
    :param path_image: 需要镶嵌影像的路径
    :return: None
    """
    path = path_image
    path_lists = os.listdir(path)
    for i in range(0, len(path_lists)):
        print("正在处理第%s幅影像......" % i)
        if i == 0:
            continue
        elif i == 1:
            img1 = gdal.Open(path+path_lists[0], gdal.GA_ReadOnly)
            img2 = gdal.Open(path+path_lists[1], gdal.GA_ReadOnly)
            input_proj = img2.GetProjection()
            options = gdal.WarpOptions(srcSRS=input_proj, dstSRS=input_proj, format='GTiff',
                                       resampleAlg=gdal.GRA_NearestNeighbour, callback=Show_Progress)
            # 输入投影,输出投影,输出格式,重采样方法
            gdal.Warp(path+"%s.tif" % (i+1), [img1, img2], options=options)  # 输出路径,需要镶嵌的数据,参数配置
            img1 = None
            img2 = None
            del img1, img2
        else:
            img1 = gdal.Open(path+"%s.tif" % i, gdal.GA_ReadOnly)
            img2 = gdal.Open(path + path_lists[i], gdal.GA_ReadOnly)
            input_proj = img2.GetProjection()
            options = gdal.WarpOptions(srcSRS=input_proj, dstSRS=input_proj, format='GTiff',
                                       resampleAlg=gdal.GRA_NearestNeighbour, callback=Show_Progress)
            # 输入投影,输出投影,输出格式,重采样方法
            gdal.Warp(path + "%s.tif" % (i + 1), [img1, img2], options=options)  # 输出路径,需要镶嵌的数据,参数配置
            img1 = None
            img2 = None
            del img1, img2

        使用时直接修改Mosaic_GDAL函数的入参就行了,选择数据存放的路径会自动拼接,命名也会自己设置无需额外修改。

        如果大家在学习Python或者RS时有什么问题,可以随时留言交流!如果大家对批量处理有兴趣同样可以留言给博主,博主会分享相关代码以供学习!

以下是一个基于GDAL库的python代码,用于批量遥感影像栅格化为矢量文件: ```python import os import glob import gdal from osgeo import ogr path = 'path/to/raster/files' output_path = 'path/to/vector/files' def raster_to_vector(raster_file, vector_file): # Open raster file and get metadata ds = gdal.Open(raster_file) band = ds.GetRasterBand(1) geotransform = ds.GetGeoTransform() proj = ds.GetProjection() cols = ds.RasterXSize rows = ds.RasterYSize # Create output shapefile and layer driver = ogr.GetDriverByName('ESRI Shapefile') out_ds = driver.CreateDataSource(vector_file) out_layer = out_ds.CreateLayer(vector_file, srs=ogr.osr.SpatialReference(proj), geom_type=ogr.wkbPolygon) # Create field for pixel value field_def = ogr.FieldDefn('Value', ogr.OFTInteger) out_layer.CreateField(field_def) # Create polygons from raster cells for y in range(rows): for x in range(cols): # Get pixel value and check if it is nodata value = band.ReadAsArray(x, y, 1, 1)[0, 0] if value == band.GetNoDataValue(): continue # Calculate polygon vertices ulx, xres, xskew, uly, yskew, yres = geotransform xcoord = ulx + (x + 0.5) * xres ycoord = uly + (y + 0.5) * yres ring = ogr.Geometry(ogr.wkbLinearRing) ring.AddPoint(xcoord - 0.5 * xres, ycoord - 0.5 * yres) ring.AddPoint(xcoord + 0.5 * xres, ycoord - 0.5 * yres) ring.AddPoint(xcoord + 0.5 * xres, ycoord + 0.5 * yres) ring.AddPoint(xcoord - 0.5 * xres, ycoord + 0.5 * yres) ring.AddPoint(xcoord - 0.5 * xres, ycoord - 0.5 * yres) poly = ogr.Geometry(ogr.wkbPolygon) poly.AddGeometry(ring) # Create feature and add to layer feature = ogr.Feature(out_layer.GetLayerDefn()) feature.SetGeometry(poly) feature.SetField('Value', value) out_layer.CreateFeature(feature) # Clean up ds = None out_ds = None # Loop through raster files for raster_file in glob.glob(os.path.join(path, '*.tif')): # Create output vector file name vector_file = os.path.join(output_path, os.path.splitext(os.path.basename(raster_file))[0] + '.shp') # Convert raster to vector raster_to_vector(raster_file, vector_file) ``` 这个代码将遥感影像文件夹中的所有.tif文件转换为矢量文件,输出到指定的路径下。可以根据需要对输出的矢量文件进行进一步处理,如合并、剪裁等操作。
评论 18
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

RS迷途小书童

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值