遥感影像处理全流程
一．预处理
二．几何纠正
三．图像裁剪
四．图像镶嵌和匀色
五．图像增强
六．遥感信息提取
参考资料

一．预处理
1．降噪处理

由于传感器的因素，一些获取的遥感图像中，会出现周期性的噪声，我们必须对其进行消除或减弱方可使用。

（1）除周期性噪声和尖锐性噪声

周期性噪声一般重叠在原图像上，成为周期性的干涉图形，具有不同的幅度、频率、和相位。它形成一系列的尖峰或者亮斑，代表在某些空间频率位置最为突出。一般可以用带通或者槽形滤波的方法来消除。

消除尖峰噪声，特别是与扫描方向不平行的，一般用傅立叶变换进行滤波处理的方法比较方便。
图1 消除噪声前
图2 消除噪声后
（2）除坏线和条带

去除遥感图像中的坏线。遥感图像中通常会出现与扫描方向平行的条带，还有一些与辐射信号无关的条带噪声，一般称为坏线。一般采用傅里叶变换和低通滤波进行消除或减弱。
图3 去条纹前
图4 去条纹后
图5 去条带前
图6 去条带后
2．薄云处理

由于天气原因，对于有些遥感图形中出现的薄云可以进行减弱处理。

3．阴影处理

由于太阳高度角的原因，有些图像会出现山体阴影，可以采用比值法对其进行消除。
二．几何纠正
通常我们获取的遥感影像一般都是Level2级产品，为使其定位准确，我们在使用遥感图像前，必须对其进行几何精纠正，在地形起伏较大地区，还必须对其进行正射纠正。特殊情况下还须对遥感图像进行大气纠正，此处不做阐述。

1．图像配准

为同一地区的两种数据源能在同一个地理坐标系中进行叠加显示和数学运算，必须先将其中一种数据源的地理坐标配准到另一种数据源的地理坐标上,这个过程叫做配准。

（1）影像对栅格图像的配准

将一幅遥感影像配准到相同地区另一幅影像或栅格地图中，使其在空间位置能重合叠加显示。
图7 图像配准前
图8 图像配准后
（2）影像对矢量图形的配准

将一幅遥感影像配准到相同地区一幅矢量图形中，使其在空间位置上能进行重合叠加显示。

2．几何粗纠正

这种校正是针对引起几何畸变的原因进行的，地面接收站在提供给用户资料前，已按常规处理方案与图像同时接收到的有关运行姿态、传感器性能指标、大气状态、太阳高度角对该幅图像几何畸变进行了校正.

3．几何精纠正

为准确对遥感数据进行地理定位，需要将遥感数据准确定位到特定的地理坐标系的，这个过程称为几何精纠正。

（1）图像对图像的纠正

利用已有准确地理坐标和投影信息的遥感影像，对原始遥感影像进行纠正，使其具有准确的地理坐标和投影信息。

（2）图像对地图（栅格或矢量）

利用已有准确地理坐标和投影信息的扫描地形图或矢量地形图，对原始遥感影像进行纠正，使其具有准确的地理坐标和投影信息。
图9 参考地形图
图10 待纠正影像
图11 纠正后影像和地形图套和效果
（3）图像对已知坐标点（地面控制点）

利用已有准确地理坐标和投影信息的已知坐标点或地面控制点，对原始遥感影像进行纠正，使其具有准确的地理坐标和投影信息。

4．正射纠正

利用已有地理参考数据（影像、地形图和控制点等）和数字高程模型数据（DEM、GDEM），对原始遥感影像进行纠正，可消除或减弱地形起伏带来的影像变形，使得遥感影像具有准确的地面坐标和投影信息。
图12 数字正射影像图
三．图像裁剪
在日常遥感应用中，常常只对遥感影像中的一个特定的范围内的信息感兴趣，这就需要将遥感影像裁减成研究范围的大小。
图39 原始影像
1．按ROI裁剪

根据ROI（感兴趣区域）范围大小对被裁减影像进行裁剪。
图40 按ROI（行政区）域裁剪
2．按文件裁剪

按照指定影像文件的范围大小对被裁减影像进行裁剪。

3．按地图裁剪

根据地图的地理坐标或经纬度的范围对被裁减影像进行裁剪。
图41 按地图坐标范围裁剪
四．图像镶嵌和匀色
1．图像镶嵌

也叫图像拼接，是将两幅或多幅数字影像(它们有可能是在不同的摄影条件下获取的)拼在一起，构成一幅整体图像的技术过程。

通常是先对每幅图像进行几何校正，将它们规划到统一的坐标系中，然后对它们进行裁剪，去掉重叠的部分，再将裁剪后的多幅影像装配起来形成一幅大幅面的影像。
图42镶嵌左影像
图43 镶嵌右影像
图44 镶嵌结果影像
2．影像匀色

在实际应用中，我们用来进行图像镶嵌的遥感影像，经常来源于不同传感器、不同时相的遥感数据，在做图象镶嵌时经常会出现色调不一致，这时就需要结合实际情况和整体协调性对参与镶嵌的影像进行匀色。
图45 匀色前影像
图46 匀色后影像
五．图像增强
为使遥感图像所包含的地物信息可读性更强，感兴趣目标更突出，需要对遥感图像进行增强处理。

1．彩色合成

为了充分利用色彩在遥感图像判读和信息提取中的优势，常常利用彩色合成的方法对多光谱图像进行处理，以得到彩色图像。

彩色图像可以分为真彩色图像和假彩色图像。
图13真彩色合成( TM321)
图14 假彩色合成(TM432)
2．直方图变换

统计每幅图像的各亮度的像元数而得到的随机分布图，即为该幅图像的直方图。

一般来说，包含大量像元的图像，像元的亮度随机分布应是正态分布。直方图为非正态分布，说明图像的亮度分布偏亮、偏暗或亮度过于集中，图像的对比度小，需要调整该直方图到正态分布，以改善图像的质量。
图15 直方图拉伸前（原图偏暗）
图16 直方图拉伸后
图17 直方图拉伸前（原图对比度不强）
图18 直方图拉伸后（线性拉伸）
3．密度分割

将灰度图像按照像元的灰度值进行分级，再分级赋以不同的颜色，使原有灰度图像变成伪彩色图像，达到图像增强的目的。
图19 原始图像
图20 密度分割图像
4．灰度颠倒

灰度颠倒是将图像的灰度范围先拉伸到显示设备的动态范围(如0～255)到饱和状态，然后再进行颠倒，使正像和负像互换。
图21 灰度颠倒前
图22 灰度颠倒后
5．图像间运算

两幅或多幅单波段图像，空间配准后可进行算术运算，实现图像的增强。常见的有加法运算、减法运算、比值运算和综合运算。例如：

减法运算：可突现出两波段差值大的地物，如红外-红，可突现植被信息。

比值运算：常用于计算植被指数、消除地形阴影等。

植被指数：NDVI=（IR-R）/（IR+R）
图23 原始图像
图24 NDVI植被指数图像
6．邻域增强

又叫滤波处理，是在被处理像元周围的像元参与下进行的运算处理，邻域的范围取决于滤波器的大小，如3×3或5×5等。

邻区法处理用于去噪声、图像平滑、锐化和相关运算 。
图25 原始图像
图26 拉普拉斯滤波图像（5×5）
7．主成分分析

也叫PCA变换，可以用来消除特征向量中各特征之间的相关性，并进行特征选择。

主成分分析算法还可以用来进行高光谱图像数据的压缩和信息融合。例如：对LandsatTM的6个波段的多光谱图像(热红外波段除外)进行主成分分析，然后把得到的第1，2，3主分量图像进行彩色合成，可以获得信息量非常丰富的彩色图像。
图27 第一主成分
图28 第二主成分
图29 第三主成分
图30 第四主成分
图31第五主成分
图32 第六主成分
8．K-T变换

即Kauth-Thomas变换，又称为“缨帽变换”。这种变换着眼点在于农作物生长过程而区别于其他植被覆盖，力争抓住地面景物在多光谱空间中的特征。

目前对这个变换的研究主要集中在MSS与TM两种遥感数据的应用分析方面。
图33 第一主分量（亮度）
图34 第二主分量（绿度）
图35第三主分量
9．图像融合

遥感图像信息融合是将多源遥感数据在统一的地理坐标系中，采用一定的算法生成一组新的信息或合成图像的过程。

不同的遥感数据具有不同的空间分辨率、波谱分辨率和时相分辨率，如果能将它们各自的优势综合起来，可以弥补单一图像上信息的不足，这样不仅扩大了各自信息的应用范围，而且大大提高了遥感影像分析的精度。
图36 多光谱影像
图 37高分辨率影像
图38 融合影像（HSV融合）
六．遥感信息提取
遥感图像中目标地物的特征是地物电磁波的辐射差异在遥感影像上的反映。依据遥感图像上的地物特征，识别地物类型、性质、空间位置、形状、大小等属性的过程即为遥感信息提取。

目前信息提取的方法有：目视判读法和计算机分类法。其中目视判读是最常用的方法。

1．目视判读

也叫人工解译，即用人工的方法判读遥感影像，对遥感影像上目标地物的范围进行手工勾绘，达到信息提取的目的。
图47 人工解译水系
2．图像分类

是依据是地物的光谱特征，确定判别函数和相应的判别准则，将图像所有的像元按性质分为若干类别的过程。

（1）监督分类

在研究区域选有代表性的训练场地作为样本，通过选择特征参数（如亮度的均值、方差等），建立判别函数，对样本进行分类，依据样本的分类特征来识别样本像元的归属类别的方法。
图48 原图像
图49 监督分类图像
（2）非监督分类

没有先验的样本类别，根据像元间的相似度大小进行归类，将相似度大的归为一类的方法。

（3）其他分类方

包括神经网络分类、分形分类、模糊分类等分类方法，以及他数据挖掘方法如模式识别、人工智能等，在这里不做进一步阐述。
参考资料
https://blog.csdn.net/u010666165/article/details/81449081

0.前言
因为没有喝上“秋天的第一份奶茶”，准备来更新一篇推送。
在上一篇推文中，我展示了如何使用Python结合Landsat制作遥感影像图(Python干货 | 制作遥感影像图)。
对于Landsat数据来说，对某个区域的重访周期为16天，每个位置使用全球参考系(WRS)进行索引，即每一个位置都会对应一个Path和Row，相邻的影像之间会有部分区域是重叠的。
Fig.1 World Reference System
在某些遥感影像的应用场景中，如果我们关注的区域正好处于两景影像的交界处，如下图中的象山港，那我们就需要将影像拼接起来才可以使用。
单张影像是这样。
本文合并后是这样。
1.准备工作
相较于上一篇推送，我们这次为了实现遥感影像的镶嵌拼接，我们使用到了两个库，rasterio和gdal。
import rasterio as rioimport gdal
先介绍一下我们实现两组遥感影像拼接的思路，首先选取两景相邻的影像，分别得到他们的空间范围，再得到两景组合到一起之后的空间范围，使用gdal新建一个tif文件(数据中转用)，分别得到原来两景影像在新建的tif文件中的起始位置，将对应的数据写入新的tif文件中，即实现镶嵌拼接。
上面说的是两景影像的拼接，如果是更多影像拼接同样适用，但是现阶段的方法如果拼接多的影像的话，需要的内存空间很大，容易导致内存溢出，感兴趣的朋友可以思考一下如何高效实现多景影像的拼接。
其中还有两处关键处理，一是如何去除重叠区域的无效信息，二是重叠区域的数据如何选择。希望各位看官能从代码里面找到答案。
2.动起手来
得到输入影像的四个角点。
def tiffileList2filename(tiffileList):filename = []prefix = []for ifile in tiffileList:file0 = ifile.split("\\")[-1]prefix.append(os.path.join(ifile, file0))filename.append(os.path.join(ifile, file0) + "_B1.TIF")return filename, prefixdef get_extent(tiffileList):filename, prefix = tiffileList2filename(tiffileList)rioData = rio.open(filename[0])left = rioData.bounds[0]bottom = rioData.bounds[1]right = rioData.bounds[2]top = rioData.bounds[3]for ifile in filename[1:]:rioData = rio.open(ifile)left = min(left, rioData.bounds[0])bottom = min(bottom, rioData.bounds[1])right = max(right, rioData.bounds[2])top = max(top, rioData.bounds[3])return left, bottom, right, top, filename, prefix
得到新建tif文件的size，这里已知Landsat空间分辨率为30m，如果是其他遥感数据，需对应进行修改。
def getRowCol(left, bottom, right, top):cols = int((right - left) / 30.0)rows = int((top - bottom) / 30.0)return cols, rows
主程序，其中plot_rgb为上一篇推送中用到的函数。
if __name__ == '__main__':tiffileList = [r'PathofLandsat8\LC08_L1TP_118039_20160126_20170330_01_T1',r'PathofLandsat8\LC08_L1TP_118040_20160126_20170330_01_T1']left, bottom, right, top, filename, prefix = get_extent(tiffileList)cols, rows= getRowCol(left, bottom, right, top)bands = ['B7', 'B5', 'B3']n_bands = len(bands)arr = np.zeros((n_bands, rows, cols), dtype=np.float)# 打开一个tif文件in_ds = gdal.Open(filename[0])for i in range(len(bands)):ibands = bands[i]# 新建一个tif文件driver = gdal.GetDriverByName('gtiff')out_ds = driver.Create(ibands + 'mosaic.tif', cols, rows)# 设置tif文件的投影out_ds.SetProjection(in_ds.GetProjection())out_band = out_ds.GetRasterBand(1)# 设置新tif文件的地理变换gt = list(in_ds.GetGeoTransform())gt[0], gt[3] = left, topout_ds.SetGeoTransform(gt)# 对要拼接的影像进行循环读取for ifile in prefix:in_ds = gdal.Open(ifile + '_' + ibands + '.tif')# 计算新建的tif文件及本次打开的tif文件之间的坐标漂移trans = gdal.Transformer(in_ds, out_ds, [])# 得到偏移起始点success, xyz = trans.TransformPoint(False, 0, 0)x, y, z = map(int, xyz)# 读取波段信息fnBand = in_ds.GetRasterBand(1)data = fnBand.ReadAsArray()# 写入tif文件之前，最大值设置为255，这一步很关键data = data / 65535 * 255data[np.where(data == 255)] = 0# 影像重合部分处理，重合部分取最大值xSize = fnBand.XSizeySize = fnBand.YSizeoutData = out_band.ReadAsArray(x, y, xSize, ySize)data = np.maximum(data, outData)out_band.WriteArray(data, x, y)del out_band, out_dsfile2read = ibands + 'mosaic.tif'arr[i, :, :] = tiff.imread(file2read)os.remove(file2read)plot_rgb(arr, rgb=(0, 1, 2))
3.小结
大功告成！
因遥感数据较大，可后台回复 Landsat拼接获取数据及源码。
数据处理、分析、可视化
行业资讯、学习资料
长按关注不迷路

How to make high resolution remote sensing image dataset？
github地址
0.引言
根据自身需要制作遥感影像数据还是很具有现实意义：  
第一，高分辨率遥感影像数据集目前整体上是缺乏的，主流的有UC-Merced dataset,WHU-RS dataset,RSSCN7 dataset以及2017年由武汉大学制作的目前最大包括10000张30类的高分影像的AID数据集。  
第二，主流数据集缺乏地区特点，比如针对合肥地区的地物做分类，我希望的我数据集尽可能接近我可能使用的测试数据集。人工制作数据集的工作量虽然比较大，但是结果会相对准确。
第三，扩充地物类型，AID数据集目前也只有30类，地物类别缺乏。
第四，设计高分遥感数据集，更有利于模型训练，遥感影像数据集必须要面临的一个很重要的问题就是标签样本数据缺乏的问题，相对来说，使用迁移学习的预训练模型效果会好一些，我使用过的Inception-v3做高分分类效果就很好。
1.制作流程
1.1 目标地区的遥感影像获取
使用LocalSpaceViewer来获取谷歌影像数据：  

下载类型选择谷歌影像，这个影像是多个数据源的叠加影像，更准确。
下载级别建议采用16，这样分块后的图片大小是256*256，是非常适合做场景数据集的。
截取区域的选取和要分类的数据源以及要分类的类别有关。
1.2 遥感影像的选取原则
1.Higher intra-class variations  
2.Smaller inter-class dissimilarity  
3.Relative large-scale dataset  
以上三个原则是我们选取分类场景影像最重要的指标，也就是类间的差距要小，类内的差距要大，另外就是相对比较大的标签数据量，这三个同时考虑了模型的泛化能力。

Buildings（类内差距要大）
                                              Greenbelt                           Water                             Farmland
类间差距要小，类间相似性要高

1.3  标签分类与图片编号
尽量多分类，我这次做的合肥地区的遥感影像数据集并没有做很大，只有5类地物，每类100张左右。
图片编号可使用FasStone Image View进行统一名称替换就可以。


至此就完成了数据集的制作，自制的高分影像数据集对于我们做特定地区的图像分类问题很有意义。
我制作的合肥地区的高分辨率遥感影像数据数据集在我的github上，欢迎关注！！！

本文的文字及图片来源于网络,仅供学习、交流使用,不具有任何商业用途,版权归原作者所有,如有问题请及时联系我们以作处理
以下文章来源于腾讯云，作者：bugsuse。
0.前言
因为没有喝上“秋天的第一份奶茶”，准备来更新一篇推送。
在上一篇推文中，我展示了如何使用Python结合Landsat制作遥感影像图(Python干货 | 制作遥感影像图)。
对于Landsat数据来说，对某个区域的重访周期为16天，每个位置使用全球参考系(WRS)进行索引，即每一个位置都会对应一个Path和Row，相邻的影像之间会有部分区域是重叠的。
Fig.1 World Reference System
在某些遥感影像的应用场景中，如果我们关注的区域正好处于两景影像的交界处，如下图中的象山港，那我们就需要将影像拼接起来才可以使用。
单张影像是这样。
本文合并后是这样。
1.准备工作
相较于上一篇推送，我们这次为了实现遥感影像的镶嵌拼接，我们使用到了两个库， rasterio和gdal。
importrasterio as rioimport gdal
先介绍一下我们实现两组遥感影像拼接的思路，首先选取两景相邻的影像，分别得到他们的空间范围，再得到两景组合到一起之后的空间范围，使用gdal新建一个tif文件(数据中转用)，分别得到原来两景影像在新建的tif文件中的起始位置，将对应的数据写入新的tif文件中，即实现镶嵌拼接。
上面说的是两景影像的拼接，如果是更多影像拼接同样适用，但是现阶段的方法如果拼接多的影像的话，需要的内存空间很大，容易导致内存溢出，感兴趣的朋友可以思考一下如何高效实现多景影像的拼接。
其中还有两处关键处理，一是如何去除重叠区域的无效信息，二是重叠区域的数据如何选择。希望各位看官能从代码里面找到答案。
2.动起手来
得到输入影像的四个角点。
deftiffileList2filename(tiffileList):
filename=[]
prefix=[]for ifile intiffileList:
file0= ifile.split("\\")[-1]
prefix.append(os.path.join(ifile, file0))
filename.append(os.path.join(ifile, file0)+ "_B1.TIF")returnfilename, prefixdefget_extent(tiffileList):
filename, prefix=tiffileList2filename(tiffileList)
rioData=rio.open(filename[0])
left=rioData.bounds[0]
bottom= rioData.bounds[1]
right= rioData.bounds[2]
top= rioData.bounds[3]for ifile in filename[1:]:
rioData=rio.open(ifile)
left=min(left, rioData.bounds[0])
bottom= min(bottom, rioData.bounds[1])
right= max(right, rioData.bounds[2])
top= max(top, rioData.bounds[3])return left, bottom, right, top, filename, prefix
得到新建tif文件的size，这里已知Landsat空间分辨率为30m，如果是其他遥感数据，需对应进行修改。
```pythondefgetRowCol(left, bottom, right, top):
cols= int((right - left) / 30.0)
rows= int((top - bottom) / 30.0)return cols, rows
主程序，其中plot_rgb为上一篇推送中用到的函数。
if __name__ == '__main__':
tiffileList= [r'PathofLandsat8\LC08_L1TP_118039_20160126_20170330_01_T1',
r'PathofLandsat8\LC08_L1TP_118040_20160126_20170330_01_T1']
left, bottom, right, top, filename, prefix=get_extent(tiffileList)
cols, rows=getRowCol(left, bottom, right, top)
bands= ['B7', 'B5', 'B3']
n_bands=len(bands)
arr= np.zeros((n_bands, rows, cols), dtype=np.float)#打开一个tif文件
in_ds =gdal.Open(filename[0])for i inrange(len(bands)):
ibands=bands[i]#新建一个tif文件
driver = gdal.GetDriverByName('gtiff')
out_ds= driver.Create(ibands + 'mosaic.tif', cols, rows)#设置tif文件的投影
out_ds.SetProjection(in_ds.GetProjection())
out_band= out_ds.GetRasterBand(1)#设置新tif文件的地理变换
gt =list(in_ds.GetGeoTransform())
gt[0], gt[3] =left, top
out_ds.SetGeoTransform(gt)#对要拼接的影像进行循环读取
for ifile inprefix:
in_ds= gdal.Open(ifile + '_' + ibands + '.tif')#计算新建的tif文件及本次打开的tif文件之间的坐标漂移
trans =gdal.Transformer(in_ds, out_ds, [])#得到偏移起始点
success, xyz =trans.TransformPoint(False, 0, 0)
x, y, z=map(int, xyz)#读取波段信息
fnBand = in_ds.GetRasterBand(1)
data=fnBand.ReadAsArray()#写入tif文件之前，最大值设置为255，这一步很关键
data = data / 65535 * 255data[np.where(data== 255)] =0#影像重合部分处理，重合部分取最大值
xSize =fnBand.XSize
ySize=fnBand.YSize
outData=out_band.ReadAsArray(x, y, xSize, ySize)
data=np.maximum(data, outData)
out_band.WriteArray(data, x, y)delout_band, out_ds
file2read= ibands + 'mosaic.tif'arr[i, :, :]=tiff.imread(file2read)
os.remove(file2read)
plot_rgb(arr, rgb=(0, 1, 2))
3.小结
大功告成