学习笔记——GEE\USGS\地理空间数据云\ENVI反复横跳的心酸过程

在用GEE做无监督分类时，发现影像拼接色彩差距较大、出现明显拼接缝的问题，在尝试了网上已有的直方图匹配算法之后效果不佳且经常溢出，所以回归本心用ENVI去校正色彩。本文主要记录了在多个网站和教程之间反复跳坑的过程。。

时而支楞时而摆烂的小刘

2333人浏览 · 2023-02-23 18:16:27

时而支楞时而摆烂的小刘 · 2023-02-23 18:16:27 发布

整个问题实际上是在用GEE做无监督分类时，发现影像拼接色彩差距较大、出现明显拼接缝的问题，在尝试了网上已有的直方图匹配算法之后效果不佳且经常溢出，所以回归本心用ENVI去校正，试图获得色彩统一的整个影像。现在半失败，还在尝试当中。。。【更：晚上回去辗转反侧，突然想起C2L2本身就已经大气校正过了，which means 一顿操作猛如虎，一看结果二百五。所以以下纯当记录方法了。。结果没有实际意义。。】

在GEE上根据NDVI查找色调相近影像

不知道为啥一上来就确定8、9月份。。最后色调不一致弄得贼痛苦。。之后一定先搞7、8月份。。单纯记录一下处理过程。。

var region = ee.FeatureCollection('projects/ee-lzying1999/assets/Boundary-wgs');
// Display the sample region.
Map.centerObject(region,9)
/*
//缩放 
function applyScaleFactors(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true);
}*/

//NDVI 
function NDVI(image) { 
  return image.addBands( 
    image.normalizedDifference(["SR_B4", "SR_B3"]) 
         .rename("NDVI")); 
} 
//remove cloud by select bitwise
function bitwiseExtract(value, fromBit, toBit) {
  if (toBit === undefined) toBit = fromBit;
  var maskSize = 1 + toBit - fromBit; //位数
  var mask = ee.Number(1 << maskSize).subtract(1); 
  return value.rightShift(fromBit).bitwiseAnd(mask);
}

function cloudfree_landsat(image){
  var qa = image.select('QA_PIXEL')
  var cloudState = bitwiseExtract(qa, 3)    //5
  var cloudShadowState = bitwiseExtract(qa, 4)  //3
  var mask = cloudState.eq(0) // Clear
  .and(cloudShadowState.eq(0)) // No cloud shadow
  return image.updateMask(mask)  
}
//load landsat
var collection = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2").filterDate('2000-08-01','2000-09-30').filterBounds(region);
var preinput = collection.filterBounds(region).sort('CLOUD_COVER', false).filterMetadata('CLOUD_COVER','less_than', 3.0)
.map(cloudfree_landsat).map(NDVI)
print(preinput)
Map.addLayer(preinput.max(),{min:7000,max:16000,bands:['SR_B3','SR_B2','SR_B1']},'max')
Map.addLayer(preinput.median(),{min:7000,max:16000,bands:['SR_B3','SR_B2','SR_B1']},'median')
Map.addLayer(preinput.mean(),{min:7000,max:16000,bands:['SR_B3','SR_B2','SR_B1']},'mean')
Map.addLayer(preinput.mosaic(),{min:7000,max:16000,bands:['SR_B3','SR_B2','SR_B1']},'mosaic')
Map.addLayer(preinput.qualityMosaic('NDVI'),{min:7000,max:16000,bands:['SR_B3','SR_B2','SR_B1']},'qualitymosaic')
Map.addLayer(ee.Image().paint(region, 0, 2), {}, 'region');

先比较了不同合成方法和镶嵌的效果，最后发现mean()和mosaic()效果最好。。在去云方法上也考虑了很多种：1. 网上找的去云函数+mosaic；2. sort云量排序，筛选云量小于3%的影像之后mosaic；3. 直接qualityMosaic('NDVI')，按照NDVI值最大进行镶嵌（云的NDVI小）。还做了几种方法的组合尝试，最后发现这就是玄学。。

最后还是确定了2000.8.31、2000.9.7两天的图像能完全覆盖研究区。

USGS——上天入地苦求无门

网站：USGS网址

这个网站倒是功能很全，找了一堆教程终于搞差不多明白之后发现下载才是大坑。

最上面shapefile upload这里要求贼多：Your shapefile must be compressed into a single file using the ZIP archive file format and must contain a shape format file (.shp), shape index format file (.shx), and attribute format file (.dbf). A shapefile may also include an optional projection format file (.prj). 就是要把shp\shx\dbf\prj文件打包成zip上传，同时shp文件中不能有sub_shp，也就是要一整个完成的轮廓，且这个轮廓的点不能超过500个。。于是直接下面polygon，将右边地图放到我的研究区，直接use map搞定。

下面date range开始选了一整年的时间，后面就用了GEE发现的时间（实际上不如直接USGS一步到位，在这个上面慢慢选，还能看缩略图）。cloudcover选了0-5%。Result Options也是坑，results per page最好选大点（长时间的话，50或100），不然后面翻页了选中的足迹就会消失。

然后选择数据集。数据集还是很全的，选择了Landsat collection 2 level-1，这个collection1或2，level1或2，TOA或SR数据也是查了好久的资料，概括一下就是 Landsat collection 2 level-1数据是处理最好的。因为要2000年的数据，选择Landsat4/5 TM C2 L1。

additional criteria这里可选可不选，选了之后更准确。比如选择行列号、卫星、数据格式、collection这些。

最左侧脚印代表足迹，可以看覆盖区域；右边这个图是看卫星影像缩略图（覆盖在右侧地图上）；倒数第三个表示下载，可以选择下载格式。

product options，点开之后发现和很多教程显示的都不一样！！！本来以为是数据源的问题，重新去找了collection1的数据结果发现也不一样，没有教程里面的GEOTIFF格式选项！！后来听说好像是USGS改版了。。硬着头皮下载了完整版的，然后去ENVI打开MTL_txt（头一次知道这种格式打开居然是遥感影像O.O），发现打不开。。。然后又去求助万能的CSDN，知道要修改里面的内容，把第一行和倒数第二行的LANDSAT改成L1。然而还是打不开，又报错了，又去搜，发现要改的更多了，于是毫无求索精神的本人果断换了地理空间数据云。。。谴责自己一分钟。。。。

地理空间数据云——兜兜转转又回到了原点

没想到在全是中文的界面还操作不明白，在行列号那里输入了区间，结果转圈圈转了好久。。求助了之后才知道只能一个一个行列号输入，好在之前在USGS上已经了解的差不多了。。也算因祸得福？

找到对应的影像之后点击下载（这里还发现了问题，就是云量和USGS上的显示好像不太匹配，也不知道是咋回事）。下载之后是一个安装包，本以为是geotiff格式，没想到还是和USGS没区别。。。

好在这次用ENVI终于能打开这个文件了。

ENVI——玄学，都是玄学

将需要的影像都打开，选择linear 2%线性拉伸更好看一点。

然后就去搜怎么辐射定标、大气校正、匀色镶嵌。

辐射定标

在这一步可以subset，按照我们已有的研究区shp文件进行裁剪（范围内的辐射定标）。

选择apply FLAASH settings，以便于进行下一步的大气校正。output文件结尾是dat文件（也是头一回知道dat可以打开影像。。孤陋寡闻了。。）

将三幅影像全部进行辐射定标，之后准备大气校正。

大气校正

在右边工具栏搜索FLAASH，打开工具。input这里输入上一步的dat文件，然后按图选择，由于上一步辐射定标选择了FLAASH settings，这一步scale factor不用改动。

【1-2】output file也是dat文件，下面的directory要和上面output file文件路径一致（只需要到文件夹，路径名最好保证全部英文）。【3】选择卫星，这里选择landsatTM5。【4】放到下面详细记录。【5】选择对应影像的metadata（右键view metadata）查看时间填入。【6】选择与研究区一致的气候。【7】一般都是rural农村，根据研究区来。【8】有些教程说是气溶胶校正，可以不选，但有的又选了。这里干脆按照默认选择了。【9】多波段设置，如下。【10】save可以保存当前设置的参数，也可以不选。全部设置好之后可以点击Apply运行，需要等一段时间才能校正好。

这里展开讲讲【4】ground elevation。要注意这里的单位是km。

首先加载elevation。在右边点击compute statistics，选择elevation，然后subset，选择ROI/EVF，选择我们研究区的shp。然后一路ok。

出现研究区域elevation的统计值，这里选择平均值，单位是m。因此在输入【4】的时候，要变成1.165km。

到这一步，大气校正就做完了，此时的颜色还是相对正常的，但也能看出来有一张已经开始走上不归路了。。。

匀色镶嵌

然后。。最要命的来了。。匀色镶嵌。。这也是我从GEE转战ENVI的原因。。

【1】添加大气校正好的dat文件。【2】将影像边上的黑色——背景值设置为0。【3】设置色彩参考图层。这里要先勾选color correction中的直方图匹配才可以调整。我原本以左侧面积较大的影像作为reference，但后面镶嵌颜色很奇怪。这里就放一下用较暗的影像作为reference的结果。【4】色彩校正。直方图匹配算法，一般只校正overlap area。【5】拼接缝设置。我第一次做的时候没有设置这个，结果出现了一条很明显的拼接缝。。这次学会了，点击auto generate seamlines，自己生成了拼接缝。【6】具体设置，包括羽化等等。这里不是了解的太多，直接默认设置了。【7】先看一眼。。结果直接看破防了。。【8】输出。这次用较暗的影像作为参考的效果比之前那个好多了。看一下对比：

之前用面积较大的影像做参考（突出一个阴间）：

这次用较暗的影像做参考：

这次好像还可以诶。。

裁剪

要在上面镶嵌完的图层上进行。设置ROI——>裁剪。import vector这里最下面有一个new或active，如果换了图层还是设置一下，为new，然后裁剪。

也可以调整右边那个红红的改变颜色等等，这里为了省事，什么都没调。

这里设置background为0，也是将背景去掉。output就自己设置一下，也是dat文件。

裁剪完之后。。宝娟，你的颜色怎么成这样了。。

拉伸5%之后勉强能看。

也不知道问题出在哪。每次ENVI运行好久我就觉得软件静悄悄必定在作妖。。

后续

原本是想以此为底图来做无监督分类的，因为无监督分类对底图色彩（其实是波段）还是很敏感的。结果做一张图发现如此费劲，大大小小坑踩了一圈。这个方法还有一些可以优化的，比如修改时间。选取时间更接近的几景图像。讲道理做粗分类应该用7、8月作物生长期的最合适，也不知道为啥脑子一抽选了个这个时间。后面可能不会再搞这么复杂了，直接用GEE平台做监督分类，不需要这么多校正，模型会校正偏差。记录这个也是对最近学的东西的一个记录吧，也省得之后再去找各种各样的教程。。