目录
现有某区域一幅有投影信息的单波段遥感图像(tm00.img)、DEM数据、两幅多光谱遥感图像(tm01.img\tm02.img)及该区域的矢量数据图层文件(bound.shp),其中有一幅多光谱图像无投影信息(tm02.img),现需要对无投影信息的多光谱图像进行几何校正、镶嵌、裁切等处理,并利用矢量数据裁剪镶嵌后的多光谱数据,利用裁剪后的多光谱数据及DEM数据,提取区域植被覆盖信息。
竞赛要求:
朝北缓坡植被:NDVI大于0.25,坡度小于20°,朝北
非朝北缓坡植被:NDVI大于0.25,坡度小于20°,非朝北
陡坡植被:NDVI大于0.25,坡度大于等于20°
水体:NDVI小于等于0.25,波段4的DN值大于0且小于20
裸地:NDVI小于等于0.25,波段4的DN值大等于20
背景区:NDVI小于等于0.25,波段4的DN值等于0
要求规则表达正确,统计各类植被信息百分比,并对分类结果进行专题制图,专题要求有经纬网、图名、指北针、比例尺、图例,制图符合规范。
现有某地区2000年8月15日与2011年8月22日的两时相TM遥感图像(tm20000815.img\tm20110822.img),要求采用适当的分类方法提取该区域两时相的土地利用信息,提取两时相土地利用变化信息,要求分类合理,列出土地利用变化变化转移矩阵,并分析可能导致土地利用类型变化的主要原因,完成土地利用变更报告,制作两时相土地利用类型专题图,并对整个信息提取过程进行分析与评价。(要求用栅格数据格式进行处理)
从上述两时相土地利用类型数据中提取建设用地,进行栅矢转换,利用GIS空间叠加分析功能提取两时相建设用地空间变化信息,并计算建设用地年均变化面积,对属性数据表进行修改,增加建设用地年均变化字段,制作区域建设用地年均变化空间分布专题图。
第二届和第一届的难度还是有提升的,主要是题型的丰富程度不一样了。
题目:现有某区域一幅有投影信息的单波段遥感图像(tm00.img)、DEM数据、两幅多光谱遥感图像(tm01.img\tm02.img)及该区域的矢量数据图层文件(bound.shp),其中有一幅多光谱图像无投影信息(tm02.img),现需要对无投影信息的多光谱图像进行几何校正、镶嵌、裁切等处理,并利用矢量数据裁剪镶嵌后的多光谱数据,利用裁剪后的多光谱数据及DEM数据,提取区域植被覆盖信息。
梳理一下:
绘制流程图太影响效率了,以后不画了.
(鸽了好久啊, 我这个博客是6月份写的,现在10月份,一直拖了没有写,正巧最近校赛进前四,继续写吧,虽然现在B站还有一些博客已经写过了,但是我写这个博客的初衷就是为了提高自己的水平,所以也不算重复劳动了。当然我也会参考其他博客以及视频内容)
因为问题1数据实际上缺失(实在找不到,刚想做发现数据残缺),这里仅提供简单的操作过程和思路.
这里是以tm00.img为基准,对tm02.img(无投影信息)影像进行几何校正,那么我们知道,几何校正包括两个部分,一个是针对栅格数据的几何校正,他在ArcGIS中称为地理配准;另外一个是针对矢量数据的几何校正,它在ArcGIS中成为空间校正。
在汤国安老师的书中这么描述:
搞错了,我发现不是数据残缺,而是因为ArcGIS无法识别,这是在ArcGIS中的显示:
这是在文件资源管理器中的显示:
我给没有后缀的tm和dem加上了后缀img,这个数据可能是在ENVI中进行处理的,所以默认没有后缀但是有头文件.hdr。操作后可以正常加载了:
所以这里我们应该使用地理配准工具对tm02.img进行几何校正:
为了方便查看两幅影像的对应点,这里可以创建一个查看器窗口:
此处只演示一次如何进行添加控制点:
如此重复操作,注意控制点的位置要均匀和各个位置都需要有.
操作完成后进行更新(因为我们把自动校正关闭了需要手动校正):
另外通过下面操作可以查看RMS等信息知道当前你的控制点选择的质量如何:
这里在镶嵌之前存在一个问题,就是图像中存在大块的黑色部分(我猜测可能在ENVI中这一一部分是显示为NoData,但是在ArcGIS中显示为正常值):
为什么要处理黑色部分?因为现在我们查看其像素值发现是正常值,所以再进行镶嵌的时候,黑色部分会作为正常部分与其他真正正常的部分进行镶嵌最后导致镶嵌的结果错误。所以这里需要将黑色部分全部处理为NoData.当然,如果你从一开始就使用ENVI处理可能不需要这么处理。
我们可以大胆猜测,最初的.hdr文件应该是将0作为忽略值也就是视为NoData了。
因此我们需要将所有的0值视为NoData或者赋值为NoData,这里用栅格计算器一般不可行,因为栅格计算器一般是用于单波段的影像的算术运算,对于这种虽然可以做,但是非常繁琐。
这里我们直接使用复制栅格工具,将0作为NoData值进行另存为:
现在就是正常的了:
现在在tm02_nodata.tif中左侧存在显眼的的蓝边,与周围格格不入:
最后镶嵌完会有比较明显的边界:
所以考虑将这一部分裁剪出去(这里不再演示)。
然后现在貌似有两种镶嵌方法(其实可能是同一种,仅仅是实现方式不一样)
直接使用镶嵌至新栅格工具,操作如下:
由于该方法比较简陋,效果一般,这里不展示结果
在地理数据库中创建镶嵌数据集,需要执行什么操作都在镶嵌数据集中进行(实际上也是调用各种工具,但是贵在它把各种需要用的功能都给你整合在镶嵌数据集这里)
弄完之后如下:
还需要平衡一下色彩:
必须放大了看(不知道为什么不放大看还是会有边界感,不过最后导出没有问题,可能只是显示存在一些问题)
最后导出:
最后镶嵌的效果:
这里的裁剪应该是指的掩膜。
此处的DEM应该也是存在黑色部分(实际为无效值):
通过复制栅格工具解决该问题(不演示),结果如下。
这里确实使用ArcGIS进行分类会比较繁琐,如果使用ENVI的专家决策树会快很多。不过后续的GIS“比赛都是规定ArcGIS这里一类产品,没有ENVI相关的说明,所以这里我们都使用ArcGIS实现吧。
首先,计算出NDVI、坡度、坡向:
(可以简单看出,3、2、1对应红绿蓝三个波段,所以第4个波段为近红外)
Float("tm_mask.tif - Band_4"-"tm_mask.tif - Band_3") / Float("tm_mask.tif - Band_4"+"tm_mask.tif - Band_3")
这里我们似乎通过栅格计算器是最好的选择.
朝北缓坡植被
("ndvi.tif" > 0.25) & ("slope.tif" < 20) & ("aspect.tif" == -1)
非朝北缓坡植被:
("ndvi.tif" > 0.25) & ("slope.tif" < 20) & ("aspect.tif" != -1)
陡坡植被:
("ndvi.tif" > 0.25) & ("slope.tif" >= 20)
水体:
("ndvi.tif" <= 0.25) & (("tm_mask.tif - Band_4" > 0) & ("tm_mask.tif - Band_4" < 20))
裸地:
("ndvi.tif" <= 0.25) & ("tm_mask.tif - Band_4" >= 20)
背景区:
("ndvi.tif" <= 0.25) & ("tm_mask.tif - Band_4" == 0)
我们上面其实有纰漏,不应这么写的,其实可以这么:
SetNull(~(表达式),1)
这样子后面就不需要重新将各个地类的0值重新赋值为NoData(使用栅格计算器或者重分类,这说明我们后续进行操作的时候,如果有思路就应该先把思路大部分明确,而不是做一步感觉可以做出来然后继续下一步,往往这样会出现纰漏甚至掉进编题组挖的坑里去)
得到各种地类之后,就应该统计各类植被信息百分比,然后对分类结果进行专题图的制图。
将各个地类合并,可以使用镶嵌也可以使用栅格计算器(比较麻烦但是很准确当然原理上镶嵌也可以,因为按照分类规则各个区域是不会有重叠的),这里我就使用镶嵌了,过程略。
需要统计植被信息百分比,我猜测这里应该是要统计占比面积,也许你会根据属性表中的各个类别的像元个数直接作为面积占比然后归一化,但是我认为目前影像是投影坐标系,可能在不同地理位置上的像元所代表的面积不一样,所以这里我想要通过计算面积得到更为准确的面积。
或许你会想到属性表中的计算几何,但是计算几何是矢量数据才可以使用的,不过你要是栅格转矢量然后进行可能会存在精度损失。
这里可以使用<以表格显示分区几何统计>进行面积的计算,但是该工具计算的时间稍微长了一些:
这是分区几何统计的结果:
这是单纯的像元个数:
至于后面的专题图时间有限就不进行了。
首先是需要对tm20000815.img\tm20110822.img两幅影像进行监督分类,但是ArcGIS对于监督分类的方法比较少,不过我们还是可以做的,通过影像分类工具条进行。
另外一个就是土地利用变化的变化转移矩阵如何制作?
需要说明的是土地利用变换矩阵是什么?
(参考:【GIS教程】土地利用转移矩阵、土地利用面积变化 - 知乎 (zhihu.com))
土地利用转移矩阵以矩阵的形式将俩个不同时期的土地覆盖类型之间互相转换的数量关系展现出来,可以全面的体现出一个区域土地覆盖类型的数值和转移方向。一般形式如下(下方解释说明可忽略):
其中,TL表示上一时相,TN表示下一时相。C1至Cn表示n种不同的土地覆盖类型。假设行Cn=Ci(上一时相土地类型),列Cn=Cj(下一时相土地类型)。Sij表示上一时相Ci转变为下一时相Cj的面积量,Si*表示上一时相Ci的土地覆盖类型面积的总和,S*j表示下一时相Cj的土地覆盖类型面积的总和。Si*-Sii为Ci土地类型的流出量,即上一时相Ci类型土地中转移为下一时相其它土地类型的面积总和。S*j-Sjj表示Cj土地类型的流入量,即下一时相Cj类型土地中由上一时相其它类型土地转变而来的面积总和。
上面的Sij也可以统一替换为Pij,用来展示转移量占整体的比例,能够体现整体与部分的关系,而避免因原某一土地类型总面积较小而得到转移量也较小,从而面积不能很好地表现此土地类型变化程度的情况。
那么这个 如果想到了就比较简单,应该使用面积制表是最简单的方法。
面积制表的意思稍微难理解一些,所以我简单举一些例子。
假如我们有一个四川省各个县的面要素数据,另外还有一个四川省范围的土地利用利用数据(都是Shp文件),然后我想要统计每一个县下不同土地利用类型的面积,那么就可以使用面积制表。
确实,面积制表多用于这样的场景,但是对于今天我们的两个时像的影像的处理你也可以使用面积制表,他们本质上是一样的。
这个不详细讲述了,虽然ArcGIS里面进行监督分类用的很少,但是其实和ENVI基本上差不多。
详情查看Arcmap之监督分类_arcgis监督分类-CSDN博客
这里直接展示结果:
接着进行面积制表 。
首先需要将两幅栅格图像转化为矢量文件。
另一幅略
接着,由于相同栅格值的栅格块由于位置不同会转化为多个要素(但是GridCode字段值都是一样的)。所以我们为了后面处理时间减短一些,这里我们进行融合。
另一个文件略。
现在开始面积制表:
最后就是表转excel:
然后再稍微整理一下。
至于专题图的制作这里略。
不好意思,看到第三问的要求我就知道我二问做的有点不可理喻,主要还是因为我对面积制表不太理解或者不熟练导致的:
从这里”栅矢转换“我就知道自己做错了,实际上第二问进行变换矩阵的求取不需要转化为矢量就可以我不小心忘记了,大家可以看面积制表的工具:
好了我们正式开始第三问。
第三问我感觉就是在第二问的基础上进行更深入的探索,我们原先不是得到变换矩阵就好了吗?现在不仅仅如此,你还得通过专题图的角度从可视化这个方面彻彻底底的告诉我到底哪些地方变成了什么。所以这就要求我们深入的了解面积制表背后的原理(实际上我们只是通过矢量分析的角度求解它摆了)。
那么到底如何做?
通过矢量分析得到第二问的答案并进一步得到可视化结果。
首先,通过栅格转面将两幅监督分类好的影像转化为矢量面数据,接着使用融合工具基于GridCode字段进行融合并添加字段用于描述要素的类别(方便后面可视化);
接着,使用相交工具进行进行两个面的相交(主要是相交后属性表的变化),接着基于相交结果进行属性表的编辑,最后出土即可。
栅格转面和融合在前面已经操作过了,这里就不重复进行。
这里添加字段:
接着进行相交:
[old]+"-->"+ [new]
然后稍微调整一下符号系统:
好了,其他就不过多赘述了,至于专题图这里略。