利用envi与arcmap/arcgisPro制作一张植被覆盖指数专题地图（地图学作业）
前言
本篇文章是作者的第一篇文章，对于软件的使用和对专业知识的了解都很浅薄，作为遥感专业的入门级玩家，我认为发布一篇博客来记录自己的学习过程与成果是很有必要的，一是可以重新梳理完成学习任务的基本思路，二是总结完成这次任务中自己的收获并巩固提高，三是在写文章的过程中，会有更多被略过的知识点重新展现在面前，会发现更多新的知识点。文章漏洞百出是一定的，但是只要有总结提高，就是值得的。最后，希望大家能包容文章中的漏洞并且向我提出，这对我很有帮助，也希望我的这篇文章能够帮助大家一起进步，学习简单的遥感数据处理与GIS制图。预祝大家学习愉快。
备注：博主来自山东科技大学，因此，这次的地图数据区选取的是青岛市黄岛区。
数据准备
我们需要准备的实验数据很简单，有以下这两个就够了：
- Landsat8遥感影像数据
- shape file文件
相信大家都知道该去哪儿下载遥感数据，但是，遥感影像数据的下载面临一个问题，就是下载速度太慢，我自己有的遥感数据也不是很多，如果有同学仅将遥感数据用于学习的话，可以在QQ上找到我，我可以发给你一份简单的练习数据。QQ：1337250362
shp文件大家可以去规划云自主选择并下载对应的json文件，下载好后可以来到这个网址进行文件格式的转换（json转shp）。同时，给大家避个坑，我在第一次学习打开shp文件的时候，仅下载了shp文件，而其他的文件没有下载（完整版参考下图），仅有shp文件，是无法在GIS中打开的。
利用envi进行数据处理
Lc8数据的打开
首先打开envi，相信大家已经忘记要怎么打开遥感数据了，这里直接手把手教给大家！
依次进入如下板块：文件-打开为-光学传感器-Landsat8-带有GeoTIFF的元数据
进入后打开Landsat8数据文件夹里面的MTL.txt文件
辐射定标
数据正常打开后，我们来到工具栏
这里选择辐射校正中的辐射定标工具，表现为上图蓝色高亮位置，双击打开就行，打开后来到一个GUI界面，如下图所示：
这里我们选择有七个波段的图层就可以。点击OK后，又出现下一个GUI：
图像裁剪
下面来到图像裁剪，这里，我们就要用到shapefile文件了。
在envi中导入shape file：文件-打开-xxx.shp，打开后呈现如下图所示效果：
图像不拟合很正常，黄岛区一直在填海造陆，遥感数据比较老，shape file比较新，所以边缘不是重合的，不过不影响，后面计算ndvi的时候这些都不会对我们的数据造成影响。
我们需要一个evf文件，在最后的大气校正环节，没有这个文件会报错
首先打开ENVI Classic，File – Open Vector File
双击进入后，打开shp文件所在的文件夹，你会惊喜的发现，里面是空的。其实不是没有shp文件，是因为它默认只显示evf文件，因此，我们需要点击右下角，将显示文件更改为shp或者All files
打开后将投影方式选择为UTM，将投影坐标更改为WGS-84，Zone修改为51如图所示
看到轮廓的效果图后，ENVI Classic所有的界面就可以关掉了。
回到envi
看到我们的工具栏，进入感兴趣去中的Subset Data From ROIs工具
选择我们定标后的.dat文件（这里的dat文件和shp文件一定都要在envi中打开，不然envi找不到）
选中shape file文件，Mask pixels outside of ROIs选择Yes，下面一行的参数0表示为白色北京，1表示为黑色背景，然后选择一个文件夹保存文件就好了。
在envi打开我们刚刚保存的文件后就是这个效果（如果没看到这个效果图可以试着把其他图层关掉再看）
图像的裁剪就到此为止。
大气校正
终于来到了大气校正，这一步的参数有很多，大家要细心修改参数
首先还是来到工具栏
双击蓝色高亮部分打开大气校正工具
选择输入图像，为刚刚裁剪过的图像文件，我的叫pic_cut，选好后会弹出下图的对话框，选择第二个，然后ok。
将Output Reflectance File的路径设置为与Input相同，同时这里的输出作为我们的大气校正结果，应给他起个名字，Output Directory for FLAASH Files的路径也要与第一列一致。
下面修改我们的传感器类型
下面来到了最复杂的查看所属地区的高程信息。
envi中打开envi安装路径下的envi53-data，找到GMTED2010.jp2，双击打开。
找到工具箱中的compute statistics工具
选中我们的jp2文件，然后打开掩膜选项，选择创建掩膜
选项中
选项中找到引用感兴趣区（灰色高亮位置)
选中我们的老朋友shp文件后点ok
输出选择为输出到内存（memory）
接下来弹出的对话框依旧选择Memory
然后ok
然后勾选histograms，然后ok
然后会得到一系列的DN值，这里的DN值表示高程，下面是主要数据表，包括最大值最小值与均值，我们可以看到，均值为53。
现在可以先关一关没有用的图层，比如说刚刚的JP2，黑不溜秋的而且真滴麻烦，看着就闹心。
现在你是不是已经忘记我们要干什么了？没错，反正不是做大气校正。
好了废话少说，首先回到桌面，重新点击我们后台的envi图标，这样刚刚没有进行完的GUI界面又会重新回来
希望你在关掉刚刚的图层后没有因为只顾着瞎探索没有用的功能而忘记刚刚的最重要的平均高程数据（没错我给忘了，只能再来一次）
现在，你终于可以理直气壮地把0.053填进Ground Elevation。
此时，打开遥感数据中的MTL.txt
看到file_date这一行，填入遥感数据的日期和时间
大气校正的最后一步，进入Multispectral Settings…点击Kaufman-Tanre…进入下面的界面
点击defaults选择第一个，然后ok
出来后点击Apply!!!
大功告成
我感觉这是最麻烦的一部分
Band Math获取NDVI
NDVI=(NIR-R)/(NIR+R)
因此，我们的公式可以为：
(float(b5)-b4)/(b5+b4)
找到工具栏中的Band Math
双击打开，输入我们的公式然后ADD TO LIST，选中公式选OK
依次点击选中Red和Nera Infrared分别对应B4B5
我给他取名为NDVI，选中Compress，然后OK
这样就得到了NDVI的DN值图
大家可以打开十字丝，按住左键滑动一下，浏览一下DN值（十字丝是左上角工具栏里面红色那个）
我们可以发现，DN值集中在（-1，1），DN值小于0表现为水体或裸露岩石、土壤表面，大于0表现为绿色植被，值越大覆盖程度越高。
计算植被覆盖指数
接着打开我们的compute statistics工具
选择NDVI，然后OK
勾选Histograms然后ok
得到这样一张图表
我们应该设定一个置信区间，来去除异常值。我们设置5%~95%为置信区间，小于置信区间最小值的值被定义为无植被覆盖，值为0，大于置信区间最大值的被定义为全部被植被覆盖，值为1。
公式体现为：
(b1 It min)*0+(b1 It max)1+(b1 ge min and b1 le max)((b1-min)/(max-min))
It表示小于，ge表示大于，le表示小于等于
在我的表中，按照置信区间，得到：
min = -0.703938 max = 0.528594
因此，公式为：
(b1 lt -0.703938)*0+(b1 gt 0.528594)1+(b1 ge -0.703938 and b1 le 0.528594)((b1+0.703938)/(0.528594+0.703938))
在Band Math中导入这条函数语句，对NDVI文件进行处理，保存处理后的文件，取个名字并命名文件后缀为.tif
打开我们最后的tif文件
调整十字丝观察DN值，发现现在DN值全部集中在（0， 1）之间。
数据处理阶段到此结束，恭喜大家！
ArcMap or ArcGIS Pro成图
导入tif文件并将其可视化
最后啦兄弟们
打开ArcMAP，选择添加数据
添加我们的tif文件
得到如图所示的效果
单机左侧色带，选择一个你喜欢的色带，就可以得到一个没有其他地图元素的地图了，接下来，你只需要对他个性化，便能得到自己喜欢的专题图。
ArcGIS Pro同理，不过要先创建一张地图，不然tif文件没法打开。
效果自我感觉要比ArcMAP好很多，因为有更多的个性化选择。
如果你认真看完了这篇文章，并且觉得有帮助，可以积极评论留言，每一句鼓励，每一点错误的指出，都是你们对我最大的支持，谢谢大家。
文章出处登录后可见！
