【python】+【gdal】+【numpy】批量遥感影像数据处理

介绍了怎么用gdal来处理数据(主要是分块读取、存储的思想；根据目标值对影像像素进行拉伸；去除设置了dataIgnore后的“白噪点”；还有部分对影像的裁剪)，并且结合numpy对处理效率进行提高
gdal只是拿来处理数据，对于缩减时间、提高效率还得是numpy这种数据处理库
专门的东西拿来解决专门的问题！
而且对于数据量很大的文件，还得是分块读取！

分块读取数据

def read_block(file, start_index, end_index, block, function):
    """
    :param file: 要分段读取的函数
    :param function: 要进行的操作，类似回调函数，因为我这里一定要进行一些操作，所以就没有进行判断
    :param start_index: 整个影像中开始读的位置 [行号,列号]，一般是[0,0]
    :param end_index: 结束读的位置
    :param block: 分成多少行读取
    :return:
    """
    l_row = end_index[0] - start_index[0]  # 要读的行数
    l_col = end_index[1] - start_index[1]
    index = 0
    img = gdal.Open(file)

    groups = l_row // block
    excess = l_row % block
    
    while index < groups:
        print("开始读第%d片-----------" % index, datetime.datetime.now())
        start_row = start_index[0] + block * index
        block_array = img.ReadAsArray(start_index[1], start_row, l_col, block)  # 分行读取
        function(block_array, block, l_col, index, start_row)
        index += 1
        # 对这一块操作完之后就删除，避免一直留在内存里面，导致内存不够
        block_array = None
        del block_array
    
    if excess > 0:
        start_row = start_index[0] + block * index
        block_array = img.ReadAsArray(start_index[1], start_row, l_col, excess)  # 读剩下的行和列
        function(block_array, excess, l_col, index, start_row)
        block_array = None
        del block_array

开始处理

先明确自己要干什么：
因为我们处理的是高分数据，存储的像素深度是16位，我们是将影像作为地图，能看就行，不需要很多像素值，所以第一个问题：将像素深度降下来，把值变成0~255之间。
但是因为处理的是多张影像，得保证颜色统一啊，所以第二个问题：对像素值进行拉伸。
这里有个技巧，如果要自动计算拉伸范围的话，怎么找到每张影像合适的范围比较麻烦，所以我们可以通过envi来，在envi上拉伸后，查看范围值，直接把这个值作为RGB的范围就行（省了好多事~）。
就是这个值，我这里已经处理过了，没处理之前的找不到了。。。

处理好之后，发现高分影像的阴影值有些是[0,0,0]，做遥感影像应该知道，为了只看影像范围内的值，通常要设置个dataIgnore，设置0之后发现，啧，怎么那么多“白点”。所以第三个问题：去除“白点”。
差不多就这些，其他的裁剪什么的可以去看别人的文章（hhhhhh~）。

拉伸

拉伸看了很多网上的处理方式，之前一直想自动计算拉伸后的值，可以但是很麻烦，不知道这个值怎么才是最合适的，所以就用了上面说的小技巧。
思路： 原始的波段中最大的值变成你新拉伸的值的最大值，最小的值变成新拉伸的最小的值，然后因为要把值变成0-255之间（8位），所以就整体缩放到255，就是这两句（在变了区间后再进行缩放，这样得到的结果才会和envi上看上去差不多，我记得似乎是这样的，可以自己尝试一下）：

compress_scale = (new_max[band] - new_min[band]) / 255
new_array[band, :, :] = (new_array[band, :, :] - new_min[band]) / compress_scale

def linear_stretch(array, rows, cols, new_min, new_max):
    """
    拉伸
    :param array:
    :param rows:
    :param cols:
    :param new_min: 要拉伸到的范围的波段对应最小值[band1,band2,band3]
    :param new_max: 要拉伸到的范围的波段对应最大值[band1,band2,band3]
    :return:
    """
    print('线性拉伸------------', datetime.datetime.now())
    bands = 3
    new_array = np.zeros((3, rows, cols))

    band_num = [2, 1, 0]	# 因为我的波段对应R、G、B的顺序是3，2，1 如果你对应是1，2，3就换成0，1，2

    for band in range(bands):
        compress_scale = (new_max[band] - new_min[band]) / 255
        # 所以我说专门的事交给专业的做，用np把要处理的挑出来，不处理的就不管，省时省力啊！
        # 小于最小值的变成最小值
        new_array[band, :, :] = np.where(array[band_num[band], :, :] < new_min[band],
                                         new_min[band], array[band_num[band], :, :])
        # 大于最大值的变成最大值
        new_array[band, :, :] = np.where(new_array[band, :, :] > new_max[band], new_max[band], new_array[band, :, :])

        new_array[band, :, :] = (new_array[band, :, :] - new_min[band]) / compress_scale
    print('线性拉伸结束------------', datetime.datetime.now())

    del array
    return new_array

裁剪、去白点

因为涉及到的像素范围是差不多的，所以就放在一起了，也可以分开，但是似乎也没必要
思路： 高分影像的阴影里面有些值是的RGB值是[0,0,0]，设置DataIgnore的时候设置的是0，所以就会把RGB中为0的值忽略掉（变成透明的），这样会把影像内部不该忽略的值也忽略掉，所以有白点，这些白点其实就是没有像素值的地方。那把原来RGB是[0,0,0]但是不能被忽略的值变成一个接近黑色的值就不会被忽略了，所以就去除了白点了。 同理，如果是255也差不多。直接将shp文件的array和要去白点的array相加也可以去白点，但是相当于所有像素值都加了值，颜色会改变一点，不介意的话可以直接加。
在别人的帮助下发现的新思路（表示感谢），这样处理会方便很多，不用一行一行的去判断了，所以把之前的删掉了。

def remove_white_point(array, shp_array):
    """
    原理：用原数组array加上对应的shp数组后找到值为1的位置，这些位置对应的就是原数组array中应该不为0但是是0的位置，
    把这些变成和0差不多的值就行了
    :param array: 要去除白点的数组
    :param shp_array: 对应范围内的shp数组
    :return:
    """
    added_array = array * shp_array # 相乘求交，直接裁剪掉了
    index_of_1 = np.where(added_array == 1)
    array[index_of_1] = 1
    return array

更新位置

自己更新一下transform属性，保证影像显示的位置是正确的

# 更新transform
def refresh_transform(transform, row, col):
    # 新的transform只是左上角坐标不一样 行号和列号的变化
    new_transform = [transform[0] + col * transform[1], transform[1], transform[2],
                     transform[3] + row * transform[5],
                     transform[4], transform[5]]
    return new_transform

整合起来

这里保存文件的时候可以顺便把位深降了。
降位深原理： 把原本16位的像素范围（0-65535）转换成8位的像素范围（0-255）（拉伸），然后在写入文件的时候创建一个8位的文件格式，把拉伸后的数据放进去就行。

datatype = gdal.GDT_Byte  # gdal.GDT_UInt16是16位

class Handler(object):
    def __init__(self, tif_file, shp_file, save_file, cutmin, cutmax, offset, block):
        self.tif_file = tif_file
        self.shp_file = shp_file
        if self.shp_file:
            self.shp_array = read_tiff(self.shp_file)
            self.s_projection, self.s_transform, self.s_rows, self.s_cols = get_tiff_message(self.shp_file)

        self.save_file = save_file
        self.projection, self.transform, self.rows, self.cols = get_tiff_message(self.tif_file)

        self.cutmin = cutmin
        self.cutmax = cutmax
        self.dataset = None
        self.block = block
        self.offset = offset
        self.construct()

    def construct(self):
        # 创建文件
        # 用这个可以实现降位深度
        datatype = gdal.GDT_Byte  # gdal.GDT_UInt16
        driver = gdal.GetDriverByName("GTiff")
        transform = refresh_transform(self.transform, self.offset[0], self.offset[1])
        self.dataset = driver.Create(self.save_file, self.s_cols, self.s_rows, 3, datatype)
        self.dataset.SetGeoTransform(transform)  # 写入仿射变换参数
        self.dataset.SetProjection(self.projection)  # 写入投影

    def handler(self, array, block, l_col, index, *args):
        # 拉伸
        stretch_array = linear_stretch(array, block, l_col, self.cutmin, self.cutmax)
        array = None
        del array

        # 去除白点
        shp_index = index * self.block
        #这里就是把对应的shp传进去，因为没有数据了，也没测试这样能不能截取对应部分的shp，可以自己测试一下
        shp_array = self.shp_array[shp_index::, :]
        remove_white_array = remove_white_point(stretch_array, shp_array)

        # 保存
        yoff = index * self.block
        self.write(remove_white_array, yoff)

        remove_white_array = None
        del remove_white_array

    def write(self, array, yoff):
        """注意！！！保存的时候，因为我们是分块读取（读几行的全部列），
        	而创建文件的时候是创建了一个大的可以放所有文件的容器，写入的时候也是分块的，就得注意写入的位置。
        	yoff"""
        for i in range(3):  # 修改这里可以改变保存的波段数 只保存前三个
            self.dataset.GetRasterBand(i + 1).WriteArray(array=array[i], xoff=0, yoff=yoff)

        del array

调用

if __name__ == '__main__':
    process_file = r"E:\data\xxx.dat"  # 要处理的文件
    shp_file = r"D:\data\shp_to_tif\xxx.tif"  # shp转的tif文件
    save_file = r"D:\xxx.tif"  # 输出的文件
    print('处理文件------', save_file)

    # tif文件和shp文件不一定匹配，所以读tif文件的开始行列号和结束的行列号不一定都是最开始到最结尾
    # 但这里还需要更详细的判断，可以自己领悟一下，哈哈哈，特殊条件特殊判断嘛

    projection1, transform1, shp_rows, shp_cols = get_tiff_message(shp_file)
    projection, transform, tif_rows, tif_cols = get_tiff_message(process_file)

    start_index = [int(np.fabs(transform1[3] - transform[3])), int(np.fabs(transform1[0] - transform[0]))]
    end_index = [min(start_index[0] + shp_rows, tif_rows), min(start_index[1] + shp_cols, tif_cols)]

    obj = Handler(tif_file=process_file, shp_file=shp_file,
                  save_file=save_file, cutmin=[83, 156, 214], cutmax=[355, 399, 465], offset=start_index, block=1000)

    read_block(process_file, start_index=start_index, end_index=end_index, block=1000,
               function=obj.handler)

至此，差不多就这些，可能有其他更好的方式，可以一起学习一下~