利用GDAL实现对两幅大小不同的栅格影像相交部分作差值计算

项目需求

项目中需要使用对两幅行列数不同的栅格影像中相交部分做差值计算。

项目构想

• 1、换算出小的影像的坐标原点对应的位置在大影像的行列号
• 2、对两幅二维矩阵进行切片，求出相交区域
• 3、判断对应相交区域的数值，进行差值计算

项目实现

如图分析所示

代码实现

• ReadTheRaster.py

import gdal, gdalconst, numpy
import cv2
import matplotlib.pyplot as plt
class ReadRaster:
    def __init__(self, path, ):
        self.dataset = gdal.Open(path, gdal.GA_ReadOnly)
        self.rows = self.dataset.RasterYSize  # todo  图像宽度
        self.cols = self.dataset.RasterXSize  # todo  图像长度
        self.bands = self.dataset.RasterCount  # TODO 图像波段数量
        self.proj = self.dataset.GetProjection()  # todo 地图投影信息
        self.geotrans = self.dataset.GetGeoTransform()  # todo 仿射矩阵
    def getRasterInformation(self, nband):
        band = self.dataset.GetRasterBand(nband)  # 获取波段对象
        # data = band.ReadAsArray(0, 0, self.cols, self.rows).astype(numpy.float)  #获取波段信息
        data = band.ReadAsArray(0, 0, self.cols, self.rows)  # 获取波段信息
        return data
    def computedoffset(self):
        data = self.geotrans
        return data
    def computeRows(self):
        return self.rows
    def computeCols(self):
        return self.cols
    def writeRasterInformation(self, data, Savepath, nband):
        driver = self.dataset.GetDriver()
        writeable = driver.Create(Savepath, self.cols, self.rows, self.bands, gdal.GDT_Byte)  # TODO  新建数据集
        writeable.SetGeoTransform(self.geotrans)  # 写入仿射变换参数
        writeable.SetProjection(self.proj)  # 写入投影
        for i in range(nband):
            writeable.GetRasterBand(i + 1).WriteArray(data[i], 0, 0)
            writeable.GetRasterBand(i + 1).SetNoDataValue(0)  # todo 给各波段设置nodata值
            writeable.GetRasterBand(i + 1).FlushCache()  # todo 波段统计量
            print(writeable.GetRasterBand(i + 1).GetStatistics(0, 1))  # todo 计算波段统计量  输出为min\max \Mean\stddev
        return 'success'

• 调用

#TODO  RASTERFUNCTION
def computoffset(ds,db): # TODO 计算偏移量从而得出行列号
    data1=ds.computedoffset()
    data2=db.computedoffset()
    xoffset = int((data1[0] - data2[0]) / data1[1])
    yoffset = int((data1[3] - data2[3]) / data1[5])
    data = [xoffset, yoffset]
    return data
def getcludedata(n,ds,db):   # TODO 重新切片切成大小相同的
    # data = db.getRasterInformation(n)[38:25697, 97:10678]  #TODO 调整为多少行多少列  先化为相同行相同列 即二维数组切片
    data = db.getRasterInformation(n)[computoffset(ds,db)[1]:ds.computeRows()+computoffset(ds,db)[1], computoffset(ds,db)[0]:db.computeCols()]  # TODO 调整为多少行多少列  先化为相同行相同列 即二维数组切片
    return data
def getcludedata2(n,ds,db):
    data = ds.getRasterInformation(n)[0:ds.computeRows(),0:db.computeCols()-computoffset(ds,db)[0]]
    return data
def getIntersection(n,ds,db):
    data1=getcludedata2(n,ds,db)
    data2=getcludedata(n,ds,db)
    data2[data1==0]=0 # todo 将data1中为0位置的元素对应把data2中的元素也置为0
    data1[data2==0]=0 # todo 将data2中为0位置的元素对应把data1中的元素也置为0
    result=abs(data2-data1)
    return  result
    
ds = ReadTheRaster.ReadRaster(TheFirstfilepath) #TODO 文件路径
db = ReadTheRaster.ReadRaster(TheSecondefilepath)
status = ds.writeRasterInformation([getIntersection(1, ds, db), getIntersection(2, ds, db), getIntersection(3, ds, db)],outputfilepath, 3) 

效果展示