GF1 遥感影像镶嵌拼接
#拼接影像
# -*- coding: cp936 -*-
import math
import numpy as np
from osgeo import gdal, ogr
from tqdm import tqdm
#获取影像边界范围
def Get_Extent(in_fn):
'''
这个函数主要用于获取一幅遥感影像的经纬度范围
'''
ds = gdal.Open(in_fn)
geotrans = list(ds.GetGeoTransform())
xsize = ds.RasterXSize
ysize = ds.RasterYSize
min_x = geotrans[0]
max_y = geotrans[3]
max_x = geotrans[0] + xsize * geotrans[1]
min_y = geotrans[3] + ysize * geotrans[5]
ds = None
return min_x, max_y, max_x, min_y
tiffiles_path_list = [r"J:\gz\gfvwf1\GF1_WFV2_E94.6_N41.0_20240811_L1A13520258001\\ref_rpcortho46.dat",
r"J:\gz\gfvwf1\GF1_WFV2_E94.6_N41.0_20240811_L1A13520258001\ref_rpcortho.dat"]
tiff_out_path = r'ouputtif.tiff'
def mosiacmulbandtiff( tiffiles_path_list, tiff_out_path):
'''----------------第一步：在内存中创建“空”的镶嵌结果的空栅格:有多大、几个波段、数据类型是是什么、添加坐标意义（投影、地理转换信息）-------------------------'''
global out_ds
# 获取待镶嵌栅格的最大最小的坐标值
in_fn = tiffiles_path_list[0]
min_x, max_y, max_x, min_y = Get_Extent(in_fn)
for in_fn in tiffiles_path_list[1:]:
minx, maxy, maxx, miny = GetExtent(in_fn)
min_x = min(min_x, minx)
min_y = min(min_y, miny)
max_x = max(max_x, maxx)
max_y = max(max_y, maxy)
max_x=max_x+1
max_y=max_y+1
print(f'镶嵌结果图的经纬度范围是：{min_x} {min_y} {max_x} {max_y}')
# 计算镶嵌后影像的行列数
in_ds = gdal.Open(tiffiles_path_list[0])
geotrans = list(in_ds.GetGeoTransform())
width = geotrans[1]
height = geotrans[5]
print(f'像素东西方向分辨率为{width}，南北方向分辨率为{height}')
columns = math.ceil((max_x - min_x) / width)
rows = math.ceil((max_y - min_y) / (-height))
print(f'镶嵌结果图的行列数是：行：{rows} 列：{columns}')
in_band = in_ds.GetRasterBand(1)
#driver = gdal.GetDriverByName('MEM') #创建驱动程序：创建的数据集保存到内存中，而不是磁盘
driver = gdal.GetDriverByName('gtiff')
num_bands = in_ds.RasterCount
print(f'波段数为：{num_bands}')
print(f'第一张影像的投影信息是{in_ds.GetProjection()}')
out_ds = driver.Create('d://gf1mosiac.tiff', columns, rows, num_bands, in_band.DataType) #这里只能设置输出结果的形状大小，波段等，现在是没有坐标意义的
# 给结果添加坐标意义
out_ds.SetProjection(in_ds.GetProjection())
geotrans[0] = min_x
geotrans[3] = max_y
out_ds.SetGeoTransform(geotrans)
'''------------------------------------------第二步：往空栅格里面写数据----------------------------------------------'''
for b in tqdm(range(1,num_bands ), desc='正在镶嵌真彩色近红外波段影像文件'): # 从所有影像的第一个波段开始镶嵌
for tif in tiffiles_path_list:
in_tif = gdal.Open(tif)
# 计算in_tif中左上角像元对应out_ds中的行列号
trans = gdal.Transformer(in_tif, out_ds, []) # in_tif是输入栅格，out_ds是目标栅格,第三个参数设置转换器选项，这里不需要
# success接收是否转换成功，False表示计算从源栅格到目标栅格的偏移
success, xyz = trans.TransformPoint(False, 0, 0)
x, y, z = map(int, xyz)
fnband=in_tif.GetRasterBand(b)
data =fnband .ReadAsArray() #将in_tif第b波段读成数组
# 使用np.where将0替换为np.nan
#data = np.where(data == 0, np.nan, data)
print(f'data:{data}')
# 找到不包含NaN值的行的索引
#valid_rows_index = np.all(~np.isnan(data), axis=1)
# 提取有效的像元
#valid_data = data[valid_rows_index]
#print(f'valid_data:{valid_data}')
#row_difference = len(data) - len(valid_data)
#print(f'行差：{row_difference}')
#if row_difference > 1:
# row_difference = 1
#print(f'子影像的有效值在空栅格中的像素坐标为第{y + row_difference}行，第{x}列')
# 将删除NaN后的有效像元写入到输出数据集中
#out_ds.GetRasterBand(b).WriteArray(valid_data, x, y + row_difference)
#
xSize = fnband.XSize
ySize = fnband.YSize
outData = out_ds.GetRasterBand(b).ReadAsArray(x, y, xSize, ySize)
data = np.maximum(data, outData)
out_ds.GetRasterBand(b).WriteArray(data, x, y) # 将data从空栅格的第{y}行，第{x}列开始写入
del in_tif
# cutline_ds = ogr.Open(geojson)
# cutline_layer = cutline_ds.GetLayer()
#input_raster = out_ds
# gdal.PushErrorHandler('CPLQuietErrorHandler')
mosiacmulbandtiff( tiffiles_path_list, tiff_out_path)
def clipbygeojson(tiff_out_path,input_raster,geojson = r"F:\地调系统\Test_Data(Beijing)\shp\朝阳区.geojson"):
result = gdal.Warp(
tiff_out_path,
input_raster,
format='GTiff',
cutlineDSName=geojson,
cropToCutline=True,
dstNodata=0
)
result.FlushCache()