遥感图像处理—几何校正

遥感图像处理—几何校正本节将从原理和代码两个方面讲解遥感图像的几何校正。原理首先介绍几何校正的概念:在遥感成像过程中,传感器生成的图像像元相对于地面目标物的实际位置发生了挤压、扭曲、拉伸和偏移等问题,这一现象叫做几何畸变。几何畸变会给遥感图像的定量分析、变化检测、图像融合、地图测量或更新等处理带来的很大误差,所以需要

大家好,欢迎来到IT知识分享网。

  本节将从原理和代码两个方面讲解遥感图像的几何校正。

原理

  首先介绍几何校正的概念:在遥感成像过程中,传感器生成的图像像元相对于地面目标物的实际位置发生了挤压、扭曲、拉伸和偏移等问题,这一现象叫做几何畸变。几何畸变会给遥感图像的定量分析、变化检测、图像融合、地图测量或更新等处理带来的很大误差,所以需要针对图像的几何畸变进行校正,即几何校正。

  几何校正分为几何粗校正和几何精校正。粗校正是利用空间位置变化关系,采用计算公式和辅助参数进行的校正,叫做系统几何校正;精校正是在此基础上,使图像的几何位置符合某种地理坐标系统,与地图配准,调整亮度值,即利用地面控制点(GCP)做的几何精校正。

  几何校正步骤:1.空间位置的变换(像元坐标)2.像元灰度值的重新计算,即重采样。

  1. 坐标变换

  坐标变换分为直接法和间接法。

  1)直接法:从原始图像阵列出发,依次计算每个像元在输出图像中的坐标。直接法输出的像元值大小不会发生变化,但输出图像中的像元分布不均匀。

  2)间接法:从输出图像阵列出发,依次计算每个像元在原始图像中的位置,然后计算原始图像在该位置的像元值,再将计算的像元值赋予输出图像像元。此方法保证校正后的图像的像元在空间上均匀分布,但需要进行灰度重采样。该方法是最常用的几何校正方法。

遥感图像处理—几何校正

 

  由上图可见,直接法直接以原始图像的坐标为基准点,坐标偏移到校正后的图像,坐标的位置有很多出现在了像元的中间位置,所以直接输出像元值大小导致像元分布不均匀。

  而对于间接法。以输出图像的坐标为基准点,已经定义在了格点的位置上,此时反算出该点在原始图像上对应的图像坐标,坐标多数落在像元的中间位置。这里采用最邻近法、双线性内插和三次卷积法来计算该点的灰度值,达成重采样的目的。

遥感图像处理—几何校正

 

  2. 重采样

   图像数据经过坐标变换之后,像元中心的位置发生改变,其在原始图像的位置不一定是整数行\列,需要根据输出图像各像元在原始图像中对应的位置,对原始图像重采样,建立新的栅格矩阵。

   主要有最近邻法、双线性内插法、三次卷积法。原理这里不再叙述。其中最近邻法没有改变图像光谱信息,能保证定量分析的准确性,但是图像不够平滑。后两者虽然数据连续平滑了,但是光谱信息遭到了破坏,比较适合于制图,不适用于定量分析。

遥感图像处理—几何校正

 

 

 代码实现

     在熟悉了几何校正的原理后,可以动手实现这一环节。这里主要介绍三种方法,具体如下。

    1.首先介绍最基础也是最底层的方法。(这里是分步切片得到偏移后的坐标,然后重采样)

from osgeo import gdal

def get_indices(source_ds, target_width, target_height):
    source_geotransform = source_ds.GetGeoTransform()
    source_width = source_ds.GetGeoTransform[1]
    source_height = source_ds.GetGeoTransform[5]
    dx = target_width/source_width
    dy = target_height/source_height
    target_x = np.arange(dx/2, source_ds.RasterXSize, dx)
    target_y = np.arange(dy/2, source_ds.RasterYSize, dy)
    return np.meshgrid(target_x, target_y)

def bilinear(in_data, x, y):
    x-=0.5
    y-=0.5
    x0=np.floor(x).astype(int)
    x1=x0+1
    y0=np.floor(y).astype(int)
    y1=y0+1
    
    ul=in_data[y0,x0]*(y1-y)*(x1-x)
    ur=in_data[y0,x1]*(y1-y)*(x-x0)
    ll=in_data[y1,x0]*(y-y0)*(x1-x)
    lr=in_data[y1,x1]*(y-y0)*(x-x0)

    return ul+ur+ll+lr

if __name__ =="__main__":
    in_fn=r''
    out_fn=r''
    cell_size=(0.02,-0.02)

    in_ds = gdal.Open(in_fn)
    x, y=get_indices(in_ds,*cell_size)
    outdata=bilinear(in_ds.ReadAsArray(), x, y)
    
    driver=gdal.GetDriverByName('GTiff')
    rows,cols=outdata.shape
    out_ds=driver.Create(out_fn,cols,rows,1,gdal.GDT_ Int32)

    gt=list(in_ds.GetGeoTransform())
    gt[1]=cell_size[0]
    gt[5]=cell_size[1]
    out_ds.SetGeoTransform(gt)
    
    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(outdata)
    out_band.FlushCache()
    out_band.ComputeStatistics(False)

  2.第二种方法是ReadAsArray()函数,设置buf_xsize, buf_ysize这个缓冲区大小,就可以按这个大小重采样(最邻近法)。同时需要设置新的gt[1]\gt[5]以及行列数,这里代码省略,有兴趣查阅之前的gdal笔记。

  3.第三种方法是gdalwarp或gdal.ReprojectImage。这里有重采样的部分,具体如下。

from osgeo import gdal, gdalconst

# 1.设置分辨率投影
# inputfile = r'E:\桌面文件保存路径\GIS\test\quickbird_test.dat'
# outputfile = r'E:\桌面文件保存路径\GIS\test\resampling_test.tif'
#
# inputds = gdal.Open(inputfile, gdal.GA_ReadOnly)  # 待重采样图像
# band1 = inputds.GetRasterBand(1)
# xsize = inputds.RasterXSize
# ysize = inputds.RasterYSize
# nbands = inputds.RasterCount
# inputProj = inputds.GetProjection()
# inputTrans = inputds.GetGeoTransform()
#
#
# new_x = int(xsize * (inputTrans[1] / 10))  # 计算新行列数
# new_y = int(ysize * (inputTrans[5] / -10)) # 注意抵消负号
#
# print(inputTrans[1])
# print(inputTrans[5])
#
# inputTrans = list(inputTrans)
# inputTrans[1] = 10  # 分辨率
# inputTrans[5] = -10 # 注意设置为-,否则无法SetGeoTransform
#
# print(inputTrans)
#
# driver = gdal.GetDriverByName('GTiff')
# out_ds = driver.Create(outputfile, new_x, new_y, nbands, band1.DataType)
# out_ds.SetGeoTransform(inputTrans)
# out_ds.SetProjection(inputProj)
#
# gdal.ReprojectImage(inputds, out_ds, inputProj, out_ds.GetProjection(), gdalconst.GRA_Bilinear)
#
# del out_ds

# 2.从图1重采样到图2

inputfile = r'E:\桌面文件保存路径\IDL\imgs\TMSPOT\TM-30m.img' #高光谱
referencefile = r'E:\桌面文件保存路径\IDL\imgs\TMSPOT\bldr_sp.img' #高分辨率
outputfile = r'E:\test.tif'  #重采样输出图像

inputras = gdal.Open(inputfile, gdal.GA_ReadOnly)  #待重采样图像,高光谱低分辨率
inputProj = inputras.GetProjection()
inputTrans = inputras.GetGeoTransform()
nbands = inputras.RasterCount

reference = gdal.Open(referencefile, gdal.GA_ReadOnly)  #目标图像,高分辨率全色波段
referenceProj = reference.GetProjection()
referenceTrans = reference.GetGeoTransform()
bandreference = reference.GetRasterBand(1)
x = reference.RasterXSize
y = reference.RasterYSize


driver = gdal.GetDriverByName('GTiff')
output = driver.Create(outputfile, x, y, nbands, bandreference.DataType)
output.SetGeoTransform(referenceTrans)
output.SetProjection(referenceProj)

gdal.ReprojectImage(inputras, output, inputProj, referenceProj, gdalconst.GRA_Bilinear)

del output

  4.此外完整的几何校正和重采样还可以用opencv的一些api。

from matplotlib import pyplot as plt
import cv2
import numpy as np

img = cv2.imread('aier.jpg')
rows,cols = img.shape[:2]
pts1 = np.float32([[50,50],[200,50],[50,200]])  #对应校正点
pts2 = np.float32([[10,100],[200,20],[100,250]])
M = cv2.getAffineTransform(pts1,pts2)  //计算仿射变换矩阵
#第三个参数:变换后的图像大小
res = cv2.warpAffine(img,M,(rows,cols)) //根据矩阵和图像重采样
plt.subplot(121)
plt.imshow(img[:,:,::-1])

plt.subplot(122)
plt.imshow(res[:,:,::-1])  #画出图片,但不显示

plt.show() #显示

  需要注明,本文中原理部分引用了朱文泉老师著书的内容,部分代码摘自csdn博客。

    

 

免责声明:本站所有文章内容,图片,视频等均是来源于用户投稿和互联网及文摘转载整编而成,不代表本站观点,不承担相关法律责任。其著作权各归其原作者或其出版社所有。如发现本站有涉嫌抄袭侵权/违法违规的内容,侵犯到您的权益,请在线联系站长,一经查实,本站将立刻删除。 本文来自网络,若有侵权,请联系删除,如若转载,请注明出处:https://yundeesoft.com/34463.html

(0)

相关推荐

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

关注微信