介绍了怎么用gdal来处理数据，并且结合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~）。
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]
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
裁剪、去白点
因为涉及到的像素范围是差不多的，所以就放在一起了，也可以分开，但是似乎也没必要
def clip_compress(array, shp_array, rows, cols, index, cutmin, cutmax, block):
"""
裁剪啥的。。。
:param array: 读出来要处理的数据
:param shp_array: 裁剪的矢量数据转成的栅格数据，再用gdal读出来的数据，我这里转换的边界内是1，边界外是0
:param rows: 处理了多少行
:param cols:
:param index: 处理到了第几块
:param cutmin:
:param cutmax:
:param block:
:return:
"""
bands = 3
shp_row = index * block
tif_row = -1
# 把原来的tif文件波段值直接压缩 对tif的值进行改变
# 将值归一化，然后再拉伸
# 将黑点（rgb=0 的地方替换的值）
color = [1, 12, 6]
array_data = linear_stretch(array, rows, cols, cutmin, cutmax)
del array
for i in range(shp_row, rows + shp_row):
tif_row = tif_row + 1
# 一行中shp数组的最开始和最后的位置
index_of_1 = np.where(np.array(shp_array[i, :]) == 1)
if index_of_1[0].size > 0:
shp_start = index_of_1[0][0]
shp_end = index_of_1[0][index_of_1[0].size - 1]
index_of_0 = np.where(array_data[0, tif_row, :] == 0) # tif文件中这一行为0的位置
# 把边界外的值变成0 放到找0之后实现可以减少不必要的位置
for band in range(bands):
# 把边界外的值变成0
array_data[band, tif_row, :] = np.multiply(array_data[band, tif_row, :], shp_array[i, :])
# 边界内为0的点加上合适的值 边界内为0的点可以判断值在不在shp的范围内
new_index1 = np.where(index_of_0[0][:] >= shp_start)
if new_index1[0].size > 0:
new_start = new_index1[0][0]
if index_of_0[0][new_start] <= shp_end:
new_index2 = np.where(index_of_0[0][:] <= shp_end)
if new_index2[0].size > 0:
new_end = new_index2[0][new_index2[0].size - 1]
for index in range(new_start, new_end + 1):
new_col = index_of_0[0][index]
# 只需要判断这个位置是不是1就行
if shp_array[i, new_col] == 1:
for band in range(bands):
array_data[band, tif_row, new_col] = color[band]
else:
for band in range(bands):
array_data[band, tif_row, :] = 0
return array_data
更新位置
自己更新一下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
整合起来
class Handler(object):
def __init__(self, tif_file, shp_file, name, save_file, cutmin, cutmax, offset, block):
"""
:param tif_file 处理的tif文件，只是用来获取信息的
:param shp_file 矢量文件
:param name 保存的文件名
:param save_file 保存的路径地址
:param cutmin RGB的最小值
:param cutmax RGB的最大值
:param offset
:param 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.name = name
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.offset = offset
self.block = block
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 clip_compress(self, array, block, l_col, index, *args):
print("开始处理第%d片-----------" % index, datetime.datetime.now())
new_array = clip_compress(array, self.shp_array, block, l_col, index, self.cutmin, self.cutmax, self.block)
print("处理结束-------------", datetime.datetime.now())
yoff = index * self.block
self.write(new_array, yoff)
new_array = None
del new_array
"""注意！！！保存的时候，因为我们是分块读取（读几行的全部列），
而创建文件的时候是创建了一个大的可以放所有文件的容器，写入的时候也是分块的，就得注意写入的位置。"""
def write(self, array, yoff):
print("开始存储-----------", datetime.datetime.now())
for i in range(3): # 修改这里可以改变保存的波段数 只保存前三个
self.dataset.GetRasterBand(i + 1).WriteArray(array=array[i], xoff=0, yoff=yoff)
del array
print("存储结束-----------", datetime.datetime.now())
调用
if __name__ == '__main__':
process_file = r"E:dataabattt.dat" # 要处理的文件
shp_file = r"D:datashp_to_tifsss.tif" # shp转的tif文件
save_file = r"D:dataabaxxx.tif" # 输出的文件
print('处理文件------', save_file)
projection1, transform1, shp_rows, shp_cols = get_tiff_message(shp_file)
projection, transform, tif_rows, tif_cols = get_tiff_message(process_file)
# tif文件和shp文件不一定匹配，所以读tif文件的开始行列号和结束的行列号不一定都是最开始到最结尾
# 但这里还需要更详细的判断，可以自己领悟一下，哈哈哈，特殊条件特殊判断嘛
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, name="xxx",
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.clip_compress)
至此，差不多就这些，可能有其他更好的方式，可以一起学习一下~
