当前位置: 首页 > news >正文

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()

http://www.dtcms.com/a/16061.html

相关文章:

  • C++ 多态
  • 2.5 使用注解进行单元测试详解
  • 【vue3】实现pdf在线预览的几种方式
  • 论软件评估方法
  • rv1103b编译opencv
  • 家里装修想用投影仪,如何选择?装修中应该注意什么?
  • 浮点数二分
  • 数智化的力量:柏强制药构建医药高质量发展的新引擎
  • 场外个股期权下单后多久成交?场外个股期权对投资组合的影响
  • 《pyqt+open3d》第三章——icp配准点对面——稳健性提升
  • 【前端框架】Vue3 中 `setup` 函数的作用和使用方式
  • 【Elasticsearch】Token Graphs
  • 【LeetCode】739. 每日温度
  • rust学习二、入门之运行单个脚本
  • vue前端可视化大屏页面适配方案
  • 日本90年代经济泡沫初期是什么情况?
  • 面试经典150题——字典树
  • 深入探索xtquant:账户信息查询的全面指南
  • WPS中如何批量上下居中对齐word表格中的所有文字
  • DeepSeek的深度解析:由来、研发过程、公司背景、优势、劣势与总结
  • WSL Ubuntu 安装 CUDA 教程
  • 嵌入式AI革命:DeepSeek开源如何终结GPU霸权,开启单片机智能新时代?
  • 芯谷 D2761:专为扬声器保护设计的音频限幅器
  • 网络将内网服务转换到公网上
  • vue学习10
  • Windows 图形显示驱动开发-WDDM 2.0 -GpuMmu 寻址方式
  • 蓝桥杯单片机组第十三届初赛试题-程序题(第2批)
  • 我用AI做数据分析之数据清洗
  • 【python】连接Jira获取token以及jira对象
  • C++-----------酒店客房管理系统