0%

Google Earth Engine research

目前,超级计算机和高性能计算系统逐渐性能过剩,大量级的云计算以商品的形式变得触手可及。与此同时,来自NASA、 U.S. Geological Survey 和NOAA等多个美国政府机构和European Space Agency的PB级别存档遥感数据已免费对外提供,相应的地理信息数据大量级处理技术(如TerraLib、Hadoop、GeoSpark和GeoMesa)也日臻成熟。但是要充分利用这些海量数据资源仍然有相当高的技术要求。其中一个主要障碍就是基础的IT管理:数据的获取和存储;解析模糊文件格式;管理数据库、机器分布、任务及队列、CPUs、GPU、网络;地理信息数据处理框架等等;这就导致了很多研究者难以触及这些大数据处理工具,从而限制了研究人员对于海量遥感数据的研究。

针对以上问题,Google Earth Engine(GEE)应运而生。由谷歌、卡内基梅隆大学和美国地质调查局联合开发的GEE是目前世界上先进的PB级地理数据科学分析及可视化平台。GEE面向用户提供海量卫星影像数据集与地理数据集,包括40多年历史卫星影像数据与欧空局的卫星影像数据。同时,GEE提供基于JavaScript和Python语言的API接口、分析算法与工具,方便用户实现大型数据的处理分析与信息挖掘。

简介

Google Earth Engine提供了多个类别PB级的可供分析的数据以及高性能并行计算服务,二者均可通过API获取及控制,并且GEE集成了基于Web的IDE(Interactive Development Environment)使得快速原型实现和结果可视化成为可能。

Google Earth Engine的数据库保存了大量的公共地理信息数据,包括各类卫星、航拍得到的光学和非光学波长的观测数据、环境变化、天气和气候的预报及后判、地形数据、社会经济学数据等等。所有这些数据都是经过一定预处理且保证不丢失相关信息的可直接使用形式(ready-to-use but information-preserving)。

用户可以使用和分析来自公共目录下的数据,也可以通过API调用数据库操作器来添加自己的私人数据库。该操作器是集成到大型并行处理系统的,能够自动划分算力,提供高通量的分析能力。用户既可以通过一个体量很小的Client Library,也可以通过在该Client Library基础上构建的交互开发环境来直接调用这些API。

Google earth engine编程交互界面

数据目录

Google Earth Engine公共数据目录是一个多维度PB级地理信息数据集。主要包括来自Landsat、Sentinel-1、Sentinel-2的地球观测遥感完整存档数据,以及诸如气候、地表覆盖数据等其他环境、地理信息和社会经济学数据(如表 1所示)。数据量还在以每天近6000 scenes(时延约24小时)的速度持续增加和更新。用户可以直接使用这些公共目录下的数据,也可以通过REST交互界面上传自己的私人数据且可以自由选择是否共享给他人。

Google earth engine 常用数据目录
Dataset Nominal resolution Temporal granularity Temporal coverage Spatial coverage
Landsat
Landsat 8 OLI/TIRS 30 m 16 day 2013–Now Global
Landsat 7 ETM + 30 m 16 day 2000–Now Global
Landsat 5 TM 30 m 16 day 1984–2012 Global
Landsat 4–8 surface reflectance 30 m 16 day 1984–Now Global
Sentinel
Sentinel 1 A/B ground range detected 10 m 6 day 2014–Now Global
Sentinel 2A MSI 10/20 m 10 day 2015–Now Global
MODIS
MOD08 atmosphere Daily 2000–Now Global
MOD09 surface reflectance 500 m 1 day/8 day 2000–Now Global
MOD10 snow cover 500 m 1 day 2000–Now Global
MOD11 temperature and emissivity 1000 m 1 day/8 day 2000–Now Global
MCD12 Land cover 500 m Annual 2000–Now Global
MOD13 Vegetation indices 500/250 m 16 day 2000–Now Global
MOD14 Thermal anomalies & fire 1000 m 8 day 2000–Now Global
MCD15 Leaf area index/FPAR 500 m 4 day 2000–Now Global
MOD17 Gross primary productivity 500 m 8 day 2000–Now Global
MCD43 BRDF-adjusted reflectance 1000/500 m 8 day/16 day 2000–Now Global
MOD44 veg. cover conversion 250 m Annual 2000–Now Global
MCD45 thermal anomalies and fire 500 m 30 day 2000–Now Global
ASTER
L1 T radiance 15/30/90 m 1 day 2000–Now Global
Global emissivity 100 m Once 2000–2010 Global
Other imagery
PROBA-V top of canopy reflectance 100/300 m 2 day 2013–Now Global
EO-1 hyperion hyperspectral radiance 30 m Targeted 2001–Now Global
DMSP-OLS nighttime lights 1 km Annual 1992–2013 Global
USDA NAIP aerial imagery 1 m Sub-annual 2003–2015 CONUS
Topography
Shuttle Radar Topography Mission 30 m Single 2000 60°N–54°S
USGS National Elevation Dataset 10 m Single Multiple United States
USGS GMTED2010 7.5″ Single Multiple 83°N–57°S
GTOPO30 30″ Single Multiple Global
ETOPO1 1′ Single Multiple Global
Landcover
GlobCover 300 m Non-periodic 2009 90°N–65°S
USGS National Landcover Database 30 m Non-periodic 1992–2011 CONUS
UMD global forest change 30 m Annual 2000–2014 80°N–57°S
JRC global surface water 30 m Monthly 1984–2015 78°N–60°S
GLCF tree cover 30 m 5 year 2000–2010 Global
USDA NASS cropland data layer 30 m Annual 1997–2015 CONUS
Weather, precipitation & atmosphere
Global precipitation measurement 6′ 3 h 2014–Now Global
TRMM 3B42 precipitation 15′ 3 h 1998–2015 50°N–50°S
CHIRPS precipitation 3′ 5 day 1981–Now 50°N–50°S
NLDAS-2 7.5′ 1 h 1979–Now North America
GLDAS-2 15′ 3 h 1948–2010 Global
NCEP reanalysis 2.5° 6 h 1948–Now Global
ORNL DAYMET weather 1 km Annual 1980–Now North America
GRIDMET 4 km 1 day 1979–Now CONUS
NCEP global forecast system 15′ 6 h 2015–Now Global
NCEP climate forecast system 12′ 6 h 1979–Now Global
WorldClim 30″ 12 images 1960–1990 Global
NEX downscaled climate projections 1 km 1 day 1950–2099 North America
Population
WorldPop 100 m 5 year Multiple 2010–2015
GPWv4 30″ 5 year 2000–2020 85°N–60°S

Earth Engine运行在一个轻量级的“图像”容器中,使用基于2D网络栅格波段的简单且高度通用的数据模型。单一谱段的像素点需在数据类型、分辨率和投影上保持一致,但是图像可以包含任意数量的谱段且一幅图中的谱段不需要有统一的数据类型和投影。每幅图也具有相关键值对的元数据来存储诸如图像的拍摄位置、拍摄时间和条件等信息。

相关的图像,比如同一个Sensor产出的数据分成一组,构成一个数据集。数据集的高速筛选和分类能力使得用户可以轻松地在数以百万计的图像中搜索和选择满足特定区域、时间和其他标准的数据。

导入Earth Engine的数据都会经过预处理以达到快速并有效查询使用的目标。

首先,图片按照原有的投射和分辨率切成瓦片存储到高效且备份的瓦片数据库中。每幅瓦片尺寸为256×256,这是在载入不需要的数据和预加载额外请求数据之间采取的折衷方案。与传统的“data cube”系统不同,这种数据摄入处理能够做到不丢失数据信息:数据仍然保持其原有的投影、分辨率和bit位,从而避免了数据在重采样到一个固定网格下而产生解构导致可能会不适用于特定的应用。

另外,为了在计算开发过程中能够快速可视化,每幅图都会对应有降低分辨率的瓦片图组成的金字塔模型存储在对应瓦片图数据库中。金字塔中的每一层都是通过对上一层以1/2的比例进行缩减像素采样(向下采样)得来的。在向下采样时,连续值谱段通常采取平均采样模式,而对离散值谱段(例如分类标签)使用最小、模式化、最大或固定采样模式中的一种进行采样。当只需要一幅图的一部分数据在降低分辨率下计算时,只有最合适的金字塔层级下相关的瓦片图才会从数据库中被调取出来供使用。这种指数级缩量方式使得数据能够在不同量级都能够使用而又不用耗费大量的存储,同时也能够满足基于web的地图展示。

系统架构

Earth Engine是建立在Google data center环境下的一系列技术基础之上的。其中包括Borg集群管理系统、Bigtable、Spanner分布式数据库、Colossus(Google File System的接班者)以及用于并行管道计算的FlumeJava框架。另外,Earth Engine还可以与Google Fusion Tables(一种基于web的数据库,支持带有属性的点、线、多边形等几何类型的数据表)进行互操作。

1571105424921

Earth Engine的基本架构如图 所示。Code Editor和第三方Web Apps可以借助客户端库通过Web REST API来向系统发送交互或批量请求。即时请求通过Front End servers进行处理并向Compute Master发送下一级请求。Compute Master则负责在Compute Servers资源池中进行分布式计算的任务规划。Batch Computation的运行逻辑与此基本一致,只不过是通过FlumeJava框架来管理分布规划。支持这两个计算系统的是一系列的数据服务,包括Asset Database(存有每幅图的元数据并具有高效筛选能力)。前面提到的Borg集群管理软件则负责管理系统的每个组件和每项服务在多用户之间的负载均衡。任何单一用户的请求失败只会导致重新发起查询请求而不会对系统造成其他影响。

向Earth Engie发送的请求是基于功能进行整合和求解的。用户通过从Earth Engine 的800多个函数中选取自己需要的组合起来构建自己的数据或数据处理请求,这里面既有用于地理信息处理的简繁不等的数学函数,又有机器学习、图像处理等相关的操作函数。这个算法库使得用户可以通过图像代数方法很轻松地处理图像数据,并且支持高阶的函数比如:map()和iterate()(二者均支持对一系列图像应用随机函数)、reduce()(用来计算按照区域、滑窗、时域、光谱或其他形式组织的统计数据)。下表总结了在客户端库中包含的各种类型的函数和算法。

Function category Examples Mode of operation
Numerical operations
Primitive operations add, subtract, multiply, divide, etc. Per pixel/per feature
Trigonometric operations cos, sin, tan, acos, asin, atan, etc.
Standard functions abs, pow, sqrt, exp., log, erf, etc.
Logical operations eq, neq, gt, gte, lt, lte, and, or
Bit/bitwise operations and, or, xor, not, bit shift, etc.
Numeric casting int, float, double, uint8, etc.
Array/matrix operations
Elementwise operations (numeric operations as above) Per pixel/per feature
Array manipulation Get, length, cat, slice, sort, etc.
Array construction Identity, diagonal, etc.
Matrix operations Product, determinant, transpose, inverse, pseudoinverse, decomposition, etc.
Reduce and accumulate Reduce, accum
Machine learning
Supervised classification and regression Bayes, CART, Random Forest, SVM, Perceptron, Mahalanobis, etc. Per pixel/per feature
Unsupervised Classification K-Means, LVQ, Cobweb, etc.
Other per-pixel image operations
Spectral operations Unmixing, HSV transform, etc. Per pixel
Data masking Unmask, update mask, etc.
Visualization Min/max, color palette, gamma, SLD, etc.
Location Pixel area, pixel coordinates, etc.
Kernel operations
Convolution Convolve, blur, etc. Per image tile
Morphology Min, max, mean, distance, etc.
Texture Entropy, GLCM, etc.
Simple shape kernels Circle, rectangle, diamond, cross, etc.
Standard kernels Gaussian, Laplacian, Roberts, Sobel, etc.
Other kernels Euclidean, Manhattan and Chebyshev distance, arbitrary kernels and combinations
Other Image Operations
Band manipulation Add, select, rename, etc. Per image
Metadata properties Get, set, etc.
Derivative Pixel-space derivative, spatial gradient
Edge detection Canny, Hough transform
Terrain operations Slope, aspect, hillshade, fill minima, etc.
Connected components Components, component size
Image clipping Clip
Resampling Bilinear, bicubic, etc.
Warping Translate, changeProj
Image registration Register, displacement, displace
Other tile-based operations Cumulative cost, medial axis, reduce resolution with arbitrary reducers, etc.
Image aggregations Sample region(s), reduce region(s) with arbitrary reducers
Reducers
Simple Count, distinct, first, etc. Context-dependent
Mathematical sum, product, min, max, etc.
Logical Logical and/or, bitwise and/or
Statistical Mean, median, mode, percentile, standard deviation, covariance, histogram, etc.
Correlation Kendall, Spearman, Pearson, Sen’s slope
Regression Linear regression, robust linear regression
Geometry Operations
Types Point, LineString, Polygon, etc. Per-feature
Measurements Length, area, perimeter, distance, etc.
Constructive operations Intersection, union, difference, etc.
Predicates Intersects, contains, withinDistance, etc.
Other operations Buffer, centroid, transform, simplify, etc.
Table/collection operations
Basic manipulation Sort, merge, size, first, limit, distinct, flatten, remap, etc. Streaming
Property filtering eq, neq, gt, lt, date range, and, or, not, etc.
Spatial filtering Intersects, contains, withinDistance, etc.
Parallel processing Map, reduce, iterate
Joins Simple, inner, grouping, etc.
Vector/raster operations
Rasterization Paint/draw, distance Per tile
Spatial interpolation Kriging, IDW interpolation
Vectorization reduceToVectors Scatter/gather
Other data types
Number, string, list, dictionary, date, daterange, projection, *etc.* Context-dependent

该库的基于图像的函数的大部分是对每个像素进行代数运算,在每个频段或频段之间的基础上运行,涵盖整数和浮点数学运算,逻辑比较,位操作,类型转换,条件替换和用于处理像素数组值化的多维数组操作。另外还包括常见的像素操纵函数,例如表查找,分段线性插值,多项式求解和普遍存在的归一化差异等。该库利用几个预制的机器学习工具包,可以轻松访问20多种类型的监督分类,回归和无监督集群,以及用于准确性求解的混淆矩阵操作。对于机器视觉任务,可以使用常见的基于内核的窗口操作,例如卷积,形态操作,距离和纹理分析,以及简单的基于邻域的操作,例如梯度,斜率,宽高比和连通性。其他功能还包括图像和波段元数据操作,投影和重采样操作,屏蔽和裁剪,图像到图像位移和配准以及遥感应用常用的各种专用工具,包括约束光谱分离,区域增长和成本映射操作等等。

用户可以组合这些库函数以构建希望执行的计算描述。该计算描述最终以有向非循环图(DAG)的形式呈现给Earth Engine,其中每个节点代表一个单一函数或数据访问器的执行,并包含命名函数参数的键/值对。实质上,这是一个纯函数式编程环境,而Earth Engine利用了函数式语言常用的标准技术,例如参考透明度和惰性求值,以实现显著的优化和效率提升。

用户使用客户端库(目前可用Python和JavaScript语言)编写Earth Engine程序,允许用户使用熟悉的程序式编程范式来描述如何处理图表。客户端库为图像、集合和其他数据类型(如数字,字符串,几何和列表)提供代理对象。用户脚本操纵这些代理对象,这些对象记录操作链并将它们组装成能够表达完整计算的DAG。然后将此DAG发送到Earth Engine服务进行求解。

DAG的求解是通过一系列的图表转换实现的。转换后的子图表如果可能的话会立即进行进一步求解来实现贪婪简化,从而避免冗余计算以及任何并行计算不可用的地方。比如,在子图表中表达式3+7会立即被简化成10这种值的形式。图中的其他节点会被扩展,例如当求解一个指向图像集合的节点时,它会被扩展为一个图像序列以便后续处理批量化执行。表示复杂处理操作的节点可以采用下一节中描述的分布式处理的几种策略中的任何一种。

Earth Engine旨在支持快速,交互式探索和分析空间数据,允许用户平移和缩放结果来每次查看图像的一个子集。为此,Earth Engine使用惰性计算模型,该模型仅允许计算满足当前请求所需的特定部分来进行输出。

举个例子,比如说用户可能希望计算两个季节合成图像之间的差异,进而突出显示由于气候或者积雪引起的变化情况。最简单的实现方式就是通过Earth Engine的客户端库拉取两个复合图像的差集。

Google earth engine code拉取复合图像差集

此代码创建了两个过滤的集合,其中一个是11月,12月和1月的所有Landsat 8图像,另一个是6月,7月和8月的所有Landsat 8图像。首先计算每个数据集中每个谱段的时域插值(目的是为了最小化云和云阴影的影响),然后减去所得到的复合图像结果来计算值的变化。该计算描述的DAG表示如图 所示。

Google Earth Engine DAG计算描述

传统(非惰性)计算环境可能会在处理表达式后立即开始计算一个或两个复合材料的像素,这通常需要提前将输入数据集预处理为公共地图投影,分辨率和感兴趣区域。

相反,Earth Engine采用了不同的方法:它推迟计算输出像素,知道它清楚了解到这些结果是在什么样的情况下是必须的。例如,如果结果当前显示在交互式地图上,则地图的缩放级别和视图边界可以动态地决定输出的投影和分辨率,并且可以将像素计算限制为仅可查看的像素。或者,如果结果当前要被用作另一计算的输入,则该计算可以请求所需像素的适当投影,分辨率和边界。此信息用于动态地自动重新采样和重新投影输入数据,从而可以快速可视化结果或在更复杂的计算中使用该表达式,而无需用户预先指定需要哪些像素。默认情况下,使用输入的最近邻重采样来对所请求的输出投影进行重新投影和重采样(从每个输入的下一个最高分辨率金字塔等级中选择像素),以保持频谱完整性。但是,当用户对如何管理这种重投影有偏好时,他们可以选择精确控制投影网格,并可以选择双线性和双三次采样模式。

这种方式更有助于采用交互模式和迭代模式来开展数据挖掘和算法开发。一旦用户完成算法的开发并想要大规模应用,他们可以向Earth Engine提交一个批处理请求来计算得到完整结果并在Earth Engine中实例化为一副图像或者是可供下载的多幅图像、表格甚至是视频文件。

数据分布式处理模型

Earth Engine库中的功能使用多种内置并行化和数据分布模型来实现高性能。每种模型都针对不同的数据访问模式进行了优化。

图像瓦片处理

在遥感中使用的许多光栅处理操作是局部的:任何特定输出像素的计算仅取决于某个固定距离内的输入像素。比如例如频带数学计算或频谱解混,以及诸如卷积或纹理分析的邻域操作等针对但像素的操作。通过将区域细分为区块并独立地计算每个区域,可以很容易地并行处理这些操作。处理每个输出图块通常需要为每个输入仅检索一个或少量图块。这与金字塔输入和合理的缓存相结合,可以在任何要求的比例或投影中快速计算结果。如前所述,输入会根据需要随时重新投影以匹配投影输出的需求。当然,如果用户确定不希望使用向下采样或重新投影输入,则可以在输入的投影和比例中明确指定计算方式。

大多数基于图块的操作都是使用两种策略中的一种在Earth Engine中实现的,具体取决于它们的计算成本。 针对成本较高的操作以及一次性计算整个瓦片具有显著优势的操作,会将结果写入瓦片尺寸相匹配的输出缓冲区。 瓦片通常为256×256像素,以匹配输入预处理的瓦片大小。

对于成本低的单像素运算,是通过在一个可直接相互调用的图表中执行图像处理的界面中按照每次一个像素来进行。这种结果目的是为了充分利用这些操作在Java虚拟机(JVM)环境中执行的优势,该环境具有即时(JIT)编译器,编译器负责提取并编译重复发生的函数调用序列。结果表明,在多数情况下,原始图像操作的任意链式操作都可以像手工编译的代码一样有效地执行。

空间聚合

正如某些类别的计算本质上是局部的,其他类别本质上是非局部的,例如区域或全局统计的计算,光栅到矢量的转换,或者采样图像以训练分类器。 这些操作或它们的一部分通常仍然可以并行执行,但计算最终结果需要将许多子结果聚合在一起。 例如,计算整个图像的平均值可以通过细分图像,在每个部分上并行计算和计数,然后对这些部分和计数求和来得到所需结果。

在Earth Engine中,这些类型的计算使用分散-聚集模型作为分布式进程执行。需要执行聚合的空间区域会被划分为子区域分配给分布式算力资源池中的算力单元以便进行批量求解。 每个算力单元获取或计算所需的输入像素,然后运行所需的累积操作以计算其部分结果。 这些结果将被发送回主计算器进行此计算,该计算将它们组合并将结果转换为最终形式。 例如,当计算平均值时,每个算力单元将计算总和与计数,主运算收集并汇总这些中间结果,并以总和除以总计数得到最终计算结果。

流式聚合

处理大型遥感数据集的另一个常见操作是时间序列分析。应用在空间上的相同统计聚合操作也可以应用于计算整个图像堆栈中随着时间变化的像素级变化情况。 这些操作都是通过瓦片组合来实现的。以前述方式使用延迟图像求解并行计算得到每个瓦片的输出。在每个瓦片图内,针对每个像素都会执行聚合操作。来自输入图像集合的像素数据是批量请求的,并且通过单像素聚合器进行一次性“流式传输”。 与输出瓦片相交的所有输入处理完成后,就会在每个像素处都应用最终转换以生成输出结果。

对于具有小中间状态的聚合(比如计算最小值),该分布模型可以做到快速且有效。但是对于不具有中间状态的聚合来说,可能就会非常耗费内存(比如计算Pearson的相关性,就需要在计算最终结果之前在每个像素上都存储完整的数据序列)。不过,只要瓦片的大小是明显小于完整图像的,那么即使是非常大的数据集合,流式传输仍然可以做到非常快。例如,对于Lansat5、7、8的完整数据集堆栈,包含超过500万张图片,在任意一点都只有少于2000张瓦片的深度,平均来看其实只有500张的深度。

缓存以及常见的子表达式消除

Earth Engine中的许多处理操作的成本和数据密集度都非常高,因此避免冗余计算将回事非常有价值的。例如,在地图上查看结果的单个用户将触发对输出区块的多个独立请求,所有输出区块经常依赖于一个或多个公共的子表达式,例如大空间聚合或监督分类器的训练等。为了避免重新计算先前已经请求的值,使用子图的散列作为高速缓存键将成本较高的的中间结果存储在分布式的高速缓存中。虽然多个用户可能共享缓存中的项目,但两个独立用户独立地进行相同查询的情况并不常见。但是,单个用户在增量算法开发期间重复相同的查询并因此受益于这种缓存机制就是非常常见的了。在单个查询的分布式执行期间,高速缓存还用作共享存储器的形式,存储对应于查询的子图的中间结果。

当对相同计算的后续请求到达时,较早的计算可能已经完成或者仍可能正在进行中。在开始成本高昂的操作之前,会优先去检索缓存并返回之前计算的结果。为了处理早期计算仍在进行中的情况,所有计算都是通过少量计算主服务器发送给分布式算力单元的。这些服务器会在任何给定时刻跟踪群集中正在执行的计算。当新查询到达时,如果它依赖于某些正在进行的计算,该查询将会直接加入原始查询序列中去以等待计算完成。如果计算主机出现失败,则正在进行的计算的处理可能会丢失,这种情况下可能会允许启动冗余计算,但前提是在现有计算任务完成之前就重新请求查询。

应用

Earth Engine正广泛应用于各个领域,涵盖全球森林变化,全球地表水变化、作物产量估算、稻田制图、城市测绘、洪水测绘、火灾恢复和疟疾风险绘图等等不同主题。它还被整合到许多第三方应用中,例如分析物种栖息地范围(Map of Life)、监测气候(Climate Engine)和评估土地利用变化(Collect Earth)等等。Google Earth Engine从项目创立至今,主要面向的用户为世界各地与地理信息研究相关的院校、研究机构和科研人员。