python gdal遥感影像基础操作（读写、增加波段、添加坐标系）

1.读写

def read_img(filename):

    dataset=gdal.Open(filename)

    im_width = dataset.RasterXSize

    im_height = dataset.RasterYSize

    im_geotrans = dataset.GetGeoTransform()

    im_proj = dataset.GetProjection()

    im_data = dataset.ReadAsArray(0,0,im_width,im_height)

    del dataset 

    return im_proj,im_geotrans,im_width, im_height,im_data

 

def write_img(filename, im_proj, im_geotrans, im_data):

    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

    if len(im_data.shape) == 3:

        im_bands, im_height, im_width = im_data.shape

    else:

        im_bands, (im_height, im_width) = 1,im_data.shape 

    driver = gdal.GetDriverByName("GTiff")

    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)

    dataset.SetProjection(im_proj)

    if im_bands == 1:

        dataset.GetRasterBand(1).WriteArray(im_data)

    else:

        for i in range(im_bands):

            dataset.GetRasterBand(i+1).WriteArray(im_data[i])

2.增加波段

# -*- coding: utf-8 -*-

import os

import cv2

from osgeo import gdal

import numpy as np

import time

 

#多波段

def add_band(imgPath1, imgPath2, outpath):

    (filepath,fullname) = os.path.split(imgPath1)

    filename = fullname.split('.')[0]

    driver = gdal.GetDriverByName('GTiff')

    driver.Register()

    dataset = gdal.Open(imgPath1)

    bands = dataset.RasterCount

    geotrans = dataset.GetGeoTransform()

    proj = dataset.GetProjection()

    im_width=dataset.RasterXSize  #col

    im_height=dataset.RasterYSize  #row

    datas = dataset.ReadAsArray(0,0,1,1)

    if 'int8' in datas.dtype.name:

        datatype = gdal.GDT_Byte

    elif 'int16' in datas.dtype.name:

        datatype = gdal.GDT_UInt16

    else:

        datatype = gdal.GDT_Float32

    dataset = gdal.Open(imgPath2)

    bands2 = dataset.RasterCount

    new_bands = bands + bands2

    slice_img = driver.Create(outpath + filename + '.tif',im_width,im_height,new_bands,datatype)

    for band in range(new_bands-1):

        slice_img.GetRasterBand(band + 1).WriteArray(dataset.GetRasterBand(band + 1).ReadAsArray(0,0,im_width,im_height))

    for band in range(bands,new_bands):

        slice_img.GetRasterBand(band + 1).WriteArray(dataset.GetRasterBand(band + 1).ReadAsArray(0,0,im_width,im_height))

    slice_img.SetGeoTransform(geotrans)

    slice_img.SetProjection(proj)

 

#单波段

def add_band2(imgPath, addPath, outPath):

    driver = gdal.GetDriverByName('GTiff')

    driver.Register()

    dataset = gdal.Open(imgPath)

    bands = dataset.RasterCount

    im_width=dataset.RasterXSize  #col

    im_height=dataset.RasterYSize  #row

    geotrans = dataset.GetGeoTransform()

    proj = dataset.GetProjection()

    data = dataset.ReadAsArray(0,0,im_width,im_height)

    dataset2 = gdal.Open(addPath)

    data2 = dataset.ReadAsArray(0,0,im_width,im_height)

    slice_img = driver.Create(outPath,im_width,im_height,4,gdal.GDT_Byte)

    for band in range(3):

        slice_img.GetRasterBand(band + 1).WriteArray(dataset.GetRasterBand(band + 1).ReadAsArray(0,0,im_width,im_height))

    for band in range(3,4):

        slice_img.GetRasterBand(band+1).WriteArray(dataset2.GetRasterBand(1).ReadAsArray(0,0,im_width,im_height))

    slice_img.SetGeoTransform(geotrans)

    slice_img.SetProjection(proj)

if __name__=='__main__':

    add_band2('D:/wcs/test/1/0717/image_3/fuse3.tif', 'D:/wcs/test/1/0717/image_3/img_log.tif', 'D:/wcs/test/1/0717/image_3/fuse4.tif')

3.添加坐标系

# -*- coding: utf-8 -*-

import os

import cv2

from osgeo import gdal

import numpy as np

import time

 

def addcoor(img1,img2):

    dataset1 = gdal.Open(img1)

    geotrans = dataset1.GetGeoTransform()

    proj = dataset1.GetProjection()

    dataset2 = gdal.Open(img2,1)

    dataset2.SetGeoTransform(geotrans)

    dataset2.SetProjection(proj)

    del dataset1

    del dataset2

if __name__=='__main__':

    addcoor('D:/wcs/test/outhillshd.tif','D:/wcs/test/outhillshd_edge.tif')

    # path = 'D:/wcs/test/dem_test/sample_16bit/'

    # imgs = os.listdir(path)

    # for img in imgs:

    #   addcoor(os.path.join(path,img),'D:/wcs/test/dem_test/img/')