Google Earth Engine谷歌地球引擎直方图与时间序列图绘制

  本文主要对GEE中的依据栅格图像绘制直方图与时间序列图并调整图像可视化参数操作加以介绍。本文是谷歌地球引擎（Google Earth Engine，GEE）系列教学文章的第八篇，更多GEE文章请参考专栏：GEE学习与应用（https://blog.csdn.net/zhebushibiaoshifu/category_11081040.html）。

  首先，依据Google Earth Engine谷歌地球引擎矢量数据裁剪栅格数据中内容，我们将美国普查局（United States Census Bureau）发布的2018年主要合法边界划分数据产品导入GEE；随后，依据Google Earth Engine谷歌地球引擎数据导入与筛选显示中内容，我们将GMTED2010数据产品导入GEE；并对上述二者做好重命名。同时，依据Google Earth Engine谷歌地球引擎矢量数据裁剪栅格数据中内容，将美国中西部堪萨斯州的矢量图层提取出来。

var kansas=counties.filterMetadata("GEOID","starts_with","20");
Map.addLayer(kansas,{},"Kansas");

  首先，我们进行直方图的绘制。在这里，需要基于ui.Chart.image.histogram()函数进行绘制。

var hist=ui.Chart.image.histogram(DEM,kansas,200,20);
print(hist);

  其中，DEM即为绘制直方图所依据的栅格图像，也就是依据其数值进行绘图；kansas为需要绘图的范围，也就是依据这个范围内DEM图像的数值进行直方图绘制；200为缩放系数，亦即绘图前重采样时新空间分辨率的数值，单位为米——之所以GEE在绘图前会自动对待绘图图层做一次重采样，是因为若没有绘图前的重采样，对于空间分辨率较高的遥感影像（例如Landsat 8的30 m空间分辨率），绘图时所需要的计算时间与空间都太大了，导致绘图效率不高；最后一个20是绘制直方图的柱子数量。

  运行代码，可以看到直方图会出现在右侧“Console”一栏中。

  在这里，为了更清晰地看到ui.Chart.image.histogram()函数几个参数的作用，我们再来两组对比。

  首先，将前述缩放系数由200修改为20，也就是重采样时空间分辨率由200 m改为20 m，就可以看到出现了像元数量较多导致无法绘图的错误。

  其次，将前述柱子数量由20修改为200，可以看到新的直方图柱子数量明显增多，密密麻麻。

  同时，我们可以基于.setOptions()函数对直方图的可视化选项进行修改。

hist=hist.setOptions({
  title:"Histogram of Elevation in Kansas",
  colors:["#f6c7b6"]
});
print(hist);

  其中， title为直方图的标题，colors则为直方图颜色。重新执行代码，可以看到直方图已经具有标题且颜色发生改变。

  同时，依据Google Earth Engine谷歌地球引擎栅格数据可视化代码嵌入中内容，还可以将.setOptions()函数中的若干可视化参数首先保存在一个独立的参数hist_option中，再将hist_option放入.setOptions()函数。

var hist_option={
  title:"Histogram of Elevation in Kansas",
  colors:["#9932CC"],
  is3D:true
}
hist=hist.setOptions(hist_option);
print(hist);

  可以看到，将可视化参数放入一个单独的参数中，效果与前述代码一致。另外，这里还加入了一个is3D参数，是配置绘图结果是否为3D效果的参数；但这个参数对于饼图等比较有效，对于本文直方图而言并没有立体效果。

  接下来，我们绘制时间序列图。

  依据Google Earth Engine谷歌地球引擎数据导入与筛选显示中内容，我们将Landsat 8 Collection 1 Tier 1的大气表观反射率TOA Reflectance产品导入GEE，并对其做好重命名。同时，基于Google Earth Engine谷歌地球引擎栅格代数与NDVI计算中介绍的.select()函数，将Landsat 8 Collection 1 Tier 1大气表观反射率TOA Reflectance产品的前九个波段提取出来；这是因为接下来我们不需要对全部波段进行绘图，选取一部分波段即可

var landsat_band=landsat.select("B[1-9]");

  随后，利用地图左上角的“Add amarker”按钮，在地图任意位置选择一个点要素。

  选择完毕后可以看到点要素已经出现在“Imports”区域中。

  我们对其加以重命名。

  随后，利用.buffer()函数以选择的点为基准，绘制一个缓冲区，1000表示缓冲区长度为1000 m（相当于就是画一个以我们刚刚选择的点为圆心，1000 m为半径的圆形区域）。除了点要素，.buffer()函数还可以以线、面等要素为中心，进行缓冲区域的确立。

var area=point.buffer(1000);
Map.addLayer(area);

  绘制的圆形区域需要将地图缩放比较大时才可以看到。

  接下来，就可以利用ui.Chart.image.series()函数进行时间序列图的绘制。

var series=ui.Chart.image.series(landsat_band,area,ee.Reducer.mean(),50);
print(series);

  其中，landsat_band为需要进行绘制时间序列图的图层，area就是需要绘图的区域；这两个参数的意义分别和前述ui.Chart.image.histogram()函数的前两个参数意义比较接近。随后，ee.Reducer.mean()表示在绘制各波段的折线图时，取我们整个圆形缓冲区域各像元数值的平均；这是由于，因为我们所选的绘图区域并不单单仅有一个像元（如果只选择一个像元的话，就很容易出现部分时间该像元没有遥感影像覆盖，导致时间序列折现图出现“断线”的情况），而是一个区域；而一个区域中自然是有很多个像元了，那么这么多像元的数值取哪一个作为最终出现在时间序列图中的数值呢？一般的，我们就取区域内的平均值、极值等具有代表意义的数值进行绘图。最后，50同样是缩放系数，和前述ui.Chart.image.histogram()函数的缩放系数意义一致。

  执行代码，稍等片刻即可看到时间序列图绘制完毕。

  前面我们用了Landsat 8的9个波段进行绘图，且缓冲区域的半径为1000 m，重采样的空间分辨率（缩放系数）也是比较高的50 m，计算量比较大，导致绘图时间比较长；我们还可以对参数加以适当修改，从而提升绘图效率。

var landsat_band=landsat.select("B[1-3]");
var area=point.buffer(300);
Map.addLayer(area);
var series=ui.Chart.image.series(landsat_band,area,ee.Reducer.mean(),200);
print(series);

  例如，以上代码减少了绘图的波段数量与缓冲区域半径，同时将缩放系数的数值增大（空间分辨率降低），得到新的绘图结果如下。

  前面提到了，对于一个包含多个像元的区域，我们往往采取平均值、极值等方式进行绘图；那么我们就将ee.Reducer.mean()改为ee.Reducer.max()，绘制一个区域像元最大值对应的时间序列图。

var series=ui.Chart.image.series(landsat_band,area,ee.Reducer.max(),200);
print(series);

  可以看到，用像元最大值得到的时间序列图与前述平均值得到的时间序列图比起来，相对要高一些；但是高得也并不明显，毕竟这两幅时间序列图对应的绘图区域半径只有300 m，且重采样后的空间分辨率为200 m，即单个像元的面积在圆形区域内也显得比较大。

  此外，将鼠标放到右侧的图中，还可以交互式显示图中的具体数值。

欢迎关注公众号：疯狂学习GIS