python如何使用gdal扩充影像

8 篇文章 1 订阅
import os
import gdal
import numpy as np
import pandas as pd
import datetime
import glob

def readTif(fileName):
    dataset = gdal.Open(fileName)
    if dataset == None:
        print(fileName + "文件无法打开")
        
    width = dataset.RasterXSize 
    height = dataset.RasterYSize 
    bands = dataset.RasterCount 

    data = dataset.ReadAsArray()
    geotrans = dataset.GetGeoTransform()
    proj = dataset.GetProjection()
    del dataset
    return width, height, bands, data, geotrans, proj

def writeTiff(im_data,input_data,im_geotrans,xOffset,yOffset,im_proj, path):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32
        
    im_bands, im_height, im_width = im_data.shape
        
    if len(input_data.shape) == 3:
        input_bands, input_height, input_width = input_data.shape
    elif len(input_data.shape) == 2:
        input_data = np.array([input_data])
        input_bands, input_height, input_width = input_data.shape
    else:
        print("image is error")
            
    #创建文件
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(path, int(im_width), int(im_height), int(input_bands), datatype)
    if(dataset!= None):
        dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数
        dataset.SetProjection(im_proj) #写入投影
        
    for i in range(input_bands):   
        
        dataset.GetRasterBand(i + 1).WriteArray(input_data[i],xOffset,yOffset)

    del dataset

if __name__ == '__main__':
    start = datetime.datetime.now()

    ref_tif = r"G:\xdn\class_result\sy_index_result_pro.tif"
    samplefile = r'G:\xdn\sample\sy_sample_pro211.tif'
    path = r'G:\xdn\sample\extecd_sample5.tif'

    width_ref, height_ref, bands_ref, data_ref, geotrans_ref, proj_ref = readTif(ref_tif)
    width_s, height_s, bands_s, data_s, geotrans_s, proj_s = readTif(samplefile)

    print("开始扩充影像。。。。。。。。")
    #计算待扩充影像在扩充影像左上角的偏移量
    xOffset=int((geotrans_s[0]-geotrans_ref[0])/geotrans_ref[1])
    yOffset=int((geotrans_s[3]-geotrans_ref[3])/geotrans_ref[5])

    im_data = data_ref
    input_data = data_s
    im_geotrans = geotrans_ref
    im_proj = proj_ref
    print(xOffset,yOffset)
    print(width_ref,height_ref)
    print(width_s,height_s)

    writeTiff(im_data,input_data,im_geotrans,xOffset,yOffset,im_proj, path)
    end = datetime.datetime.now()
    print(end-start)

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值