(本文原文链接在此 How to Build a Simple Geoprocessing App with GeoTrellis and React

作者:凯利·英尼斯 (Kelly Innes)   2018 年 2 月 13 日.              翻译:Jade Wang

这篇博文是使用 GeoTrellis、Akka HTTP和 React 创建GIS处理的网页应用程序的指南。

这是用 “Azavea 10%的研究时间”(Azavea’s 10% research time)完成的研究项目。该研究项目旨在了解 GeoTrellis 以构建地理处理服务,Model My Watershed geoprocessing service 地理处理服务,从而使用 Akka HTTP 代替最初的基于 Spark 作业服务器(Spark Job Server-based)的实现。

在接下来的内容中,我将介绍为 Web 应用程序(a web application)设置开发环境,该环境允许用户在地图上绘制宾夕法尼亚州内一处地理名胜,然后从 GeoTIFF 数据集( GeoTIFF data set)中提取的 GeoTrellis 图层中检索和显示 NLCD 土地利用类型单元统计。

最后,你能通过使用 GeoTrellis设置一个运行的开发环境。

设置项目

你可以在 GitHub 存储库中找到完整的项目:

https://github.com/azavea/geotrellis-collections-api-research

设置前,你确保需要使用 macOS 或 Linux,并且还需要提前安装一些依赖(dependencies)。在此处查看要求列表(list of requirements here)。 简而言之,该应用程序需要以下依赖(dependencies):

要检查是否已安装所有内容,请从 shell 运行以下命令:

$ which make curl docker java sbt spark-submit

并验证您否看到每个列出的目录路径:

/usr/bin/make
/usr/bin/curl
/usr/local/bin/docker
/usr/bin/java
/usr/local/bin/sbt
/usr/local/bin/spark-submit

你还需要确保 Docker 当前正在运行,并且 spark-submit在你的路径上可用。

确认所有内容都已安装并运行后,复制the project repository,将 cd放入其目录,然后运行以下命令来构建所有内容:

$ make

这将启动由单独的 make 规则编排的一系列构建步骤:它将构建一个容器化的 React 应用程序客户端,下载一个 GeoTIFF 数据集,并将 GeoTIFF 处理成一个 GeoTrellis 层。

(注意:整个过程可能需要一点时间。)

构建完成后,你可以通过运行来启动项目的服务器

$ make server

该命令将在端口 7000 上启动 GeoTrellis API 服务器,在端口 9555 上启动 React 应用程序。 一旦运行,您可以在浏览器中访问 localhost:9555以加载客户端:

为了方便,让我们来看看构建步骤和一些代码。

将 GeoTIFF 摄取到 GeoTrellis 图层中

GeoTrellis 可以直接读取 GeoTIFF(GeoTrellis can read GeoTIFFs directly),但要利用其更快的栅格数据类型,需要通过通常称为“摄取”的过程将 GeoTIFF 转换为  TileLayerRDD 。

上面的 make 规则的一部分将运行make ingest,依次是:

  • 使用 cURL 下载 GeoTIFF
  • 创建一个可执行的摄取 jar 文件以通过 spark-submit运行
  • 运行摄取并将图层目录写入ingest/land-cover-data

执行摄取的实际代码位于 LandCoverIngest object中:

object LandCoverIngest {
  val localGeoTiffPath =
    new java.io.File(new java.io.File(".").getCanonicalFile,
      "land-cover-data/geotiff/nlcd_pa.tif").getAbsolutePath
  val localCatalogPath =
    new java.io.File(new java.io.File(".").getCanonicalFile,
      "land-cover-data/catalog").getAbsolutePath
  def main(args: Array[String]): Unit = {
    val conf = new SparkConf()
      .setIfMissing("spark.master", "local[*]")
      .setAppName("Ingest PA Land Cover GeoTiff")
      .set("spark.serializer", classOf[KryoSerializer].getName)
      .set("spark.kryo.registrator", classOf[KryoRegistrator].getName)
    implicit val sc = new SparkContext(conf)
    try {
      val geoTiffRDD =
        HadoopGeoTiffRDD.spatial(new Path(localGeoTiffPath))
      val (_, metadata) =
        geoTiffRDD.collectMetadata[SpatialKey](FloatingLayoutScheme(256))
      val paLandCoverLayer =
        ContextRDD(
          geoTiffRDD
            .tileToLayout(metadata, NearestNeighbor)
            .mapValues { tile => tile.convert(ByteConstantNoDataCellType) },
          metadata.copy(cellType = ByteConstantNoDataCellType))
      val paLandCoverLayerID = LayerId("nlcd-pennsylvania", 0)
      FileLayerWriter(localCatalogPath)
        .write(paLandCoverLayerID, paLandCoverLayer, ZCurveKeyIndexMethod)
    } finally {
        sc.stop()
    }
  }
}

此处LandCoverIngest 代码采用输入 GeoTIFF 路径和输出目录路径,使用 Spark 生成 RDD,然后将其写回为具有 nlcd-pennsylvania id 的 GeoTrellis 层,其中包含一些描述诸如投影、图块等内容的元数据范围和布局,以及单元格数据类型。

在ingest/land-cover-data/catalog/attributes/nlcd-pennsylvania__.__0__.__metadata.json中,你将看到一个以以下几行开头的清单文件:

[
    {
        "name": "nlcd-pennsylvania",
        "zoom": 0
    },
    {
        "header": {
            "format": "file",
            "keyClass": "geotrellis.spark.SpatialKey",
            "valueClass": "geotrellis.raster.Tile",
            "path": "nlcd-pennsylvania/0"
        },
        "metadata": {
            "extent": {
                "xmin": 1261380,
                "ymin": 1962000,
                "xmax": 1782480,
                "ymax": 2295150
            },
            "layoutDefinition": {
                "extent": {
                    "xmin": 1261380,
                    "ymin": 1957230,
                    "xmax": 1783620,
                    "ymax": 2295150
                },
                "tileLayout": {
                    "layoutCols": 68,
                    "layoutRows": 44,
                    "tileCols": 256,
                    "tileRows": 256
                }
            },
            "bounds": {
                "minKey": {
                    "col": 0,
                    "row": 0
                },
                "maxKey": {
                    "col": 67,
                    "row": 43
                }
            },
            "cellType": "int8",
            "crs": "+proj=longlat +datum=WGS84 +no_defs "
        },

在本地运行摄取 ,将层目录写回本地文件系统 ——是有点人为操作的意思。 GeoTrellis 旨在利用分布式计算集群,这也是它使用 Spark 和 Spark 的 RDD 或“弹性分布式数据集” (“Resilient Distributed Dataset”)的部分原因。

这意味着对于较大的摄取工作,可以使用Amazon AWS EMR cluster完成集群执行摄取,对国家级别的和国际级别的数据集的摄取,并创建多个缩放级别的图层,例如在需要生成视觉图块图层的情况下 。

同样,GeoTrellis 还可以使用其适当命名的 S3LayerReader and S3LayerWriter 类,将图层写入 S3 并从 S3 读取图层。

为简化起见,我使用了一个小数据集,在本地执行摄取,并将图层存储在本地。 下一步中,我将逐步设置服务以使用 GeoTrellis FileCollectionLayerReader和 Akka HTTP 查询该图层。

使用 Akka HTTP 创建 GeoTrellis 服务

该项目使用土地覆盖栅格作为其 GeoTrellis & Akka HTTP API 的数据源。最相关的 API 代码片段在这里,查看链接: API code live here,包括

服务器代码依靠 Akka HTTP 和 Spray-json 来接受 POST 请求并返回结果。 它首先声明一些 case 类以匹配输入和输出 JSON:

case class GeoJsonData(geometry: String)
case class ResponseData(response: Map[String, Int])
object RequestResponseProtocol extends DefaultJsonProtocol {
  implicit val requestFormat = jsonFormat1(GeoJsonData)
  implicit val responseFormat = jsonFormat1(ResponseData)
}

接下来,它声明了两条(two routes )支持 CORS 的路由(从浏览器查询 API 所必需的),包括一个用于健康检查的 /ping路由和一个 /panlcdcount路由,它接受地理名胜的多边形区域并返回一个响应,将 NLCD 单元类型值映射到他们在该地理名胜的领域内的计数:

post {
  path("panlcdcount") {
    entity(as[GeoJsonData]) { shape =>
      complete {
        Future {
          getPANLCDCount(shape)
        }
      }
    }
  }
}

getPANLCDCount 函数调用包含在 Future 中以避免阻塞; 该函数创建一个 GeoTrellis 可以从地理名胜区域解析的形状,从与地理名胜区域相交的栅格图层中检索,然后通过 rddCellCount 操作计算和格式化结果:

def getPANLCDCount(aoi: GeoJsonData): ResponseData = {
  val areaOfInterest = createAOIFromInput(aoi.geometry)
  val rasterLayer = fetchLocalCroppedPANLCDLayer(areaOfInterest)
  ResponseData(rddCellCount(rasterLayer, areaOfInterest))
}

rddCellCount 可能是代码中最有趣的部分:

private def rddCellCount(
  rasterLayer: TileLayerCollection[SpatialKey],
  areaOfInterest: MultiPolygon
): Map[String, Int] = {
  val init = () => new LongAdder
  val update = (_: LongAdder).increment()
  val metadata = rasterLayer.metadata
  val pixelCounts: TrieMap[Int, LongAdder] = TrieMap.empty
  rasterLayer.foreach({ case (key: SpatialKey, tile: Tile) =>
    val extent = metadata.mapTransform(key)
    val re = RasterExtent(extent, metadata.layout.tileCols,
      metadata.layout.tileRows)
    Rasterizer.foreachCellByMultiPolygon(areaOfInterest, re) { case (col, row) =>
      val pixelValue: Int = tile.get(col, row).toInt
      val acc = pixelCounts.getOrElseUpdate(pixelValue, init())
      update(acc)
    }
  })
  pixelCounts
    .map { case (k, v) => k.toString -> v.sum.toInt }
    .toMap
}

此操作接收一个瓦片图层集合和一个地理名胜区域,设置一个 LongAdder实例的 TrieMap用作多线程累加器,然后遍历与地理名胜区域相交的每个瓦片的每个单元格以确定单元格本身是否相交,形状和单元格值是什么:代表 NLCD 土地利用类型的整数。使用 TrieMap和 LongAdder使服务能够并行处理来自多个栅格图层的瓦片。

并发累加器在这里有点多余,因为这个项目只查询单个图层。 然而,对于更复杂的地理处理服务,使用 TrieMap 和 LongAdder 使开发人员能够通过在图层集合上调用 .par来使用 Scala 的并行集合( Parallel Collections)。 例如,在 Model My Watershed’s geoprocessing service 的地理处理服务中,我们使用并行集合来同时检索土地利用类型和土壤类型图层的像元计数。 除了像元计数操作之外,Model My Watershed’s geoprocessing service 还有其他几个操作,包括用于查找平均栅格像元值以及最大和最小像元值的操作,其中大部分旨在处理多个栅格图层。

随着 API 服务器的运行,客户端可以发出这样的请求:

curl --request POST \
  --url http://localhost:7000/panlcdcount \
  --header 'content-type: application/json' \
  --data '{"geometry":"{\"type\":\"Polygon\",\"coordinates\":[[[-75.915527,40.010787],[-75.695801,40.430224],[-75.322266,40.329796],[-75.487061,39.926588],[-75.915527,40.010787]]]}"}'

并收到如下所示的响应(美化后):

{
    "response": {
        "23": 42796,
        "11": 16202,
        "33": 9054,
        "22": 14163,
        "21": 139348,
        "43": 110900,
        "32": 4059,
        "82": 94575,
        "42": 59245,
        "81": 644752,
        "92": 5078,
        "41": 770469,
        "91": 2563,
        "85": 6853
    }
}

在响应中,每个键代表根据 NLCD Land Cover Class Definitions的土地覆盖类型; 每个值代表风景名胜区域内该土地覆盖类型的像元计数。 您可以在此处查看 GeoTIFF 的原始元数据( GeoTIFF’s original metadata here)。

使用 React 和 Leaflet 构建前端客户端

作为前端客户端,我使用我在 “Getting Started with React and Leaflet” 中描述的模式以略复杂的方式创建了一个 React 应用程序。客户端使用通过 Webpack 构建的 React、Redux 和 Leaflet 以及  Leaflet.draw使用户能够创建地理名胜的区域,使用 axios 向地理处理服务发出 API 请求,并使用 Victory将结果显示在饼图中。

您可以在此处查看客户端应用程序源(the client application source here);一些点包括:

Leaflet.draw 的设置方式类似于我在 React/Leaflet tutorial教程中描述的模式。在这种情况下,我添加了一个 Leaflet.draw 处理程序并设置地图以侦听 Leaflet 地图上的 drawstart和 drawcreated事件,在单击应用程序标题中的铅笔图标时触发 drawstart

并将 submitAreaOfInterest函数dispatch到 Redux 以提交绘制完成后绘制图层:

componentDidMount() {
    const {
        map: {
            leafletElement: leafletMap,
        },
        props: {
            dispatch,
        },
    } = this;
    esri.basemapLayer('Imagery').addTo(leafletMap);
    leafletMap.on('draw:drawstart', () => {
        dispatch(clearData());
        dispatch(clearAPIError());
    });
    leafletMap.on('draw:created', this.onCreate);
    this.polygonDrawHandler = new L.Draw.Polygon(leafletMap);
}
componentWillReceiveProps({ drawingActive }) {
    if (drawingActive) {
        this.polygonDrawHandler.enable();
    } else {
        this.polygonDrawHandler.disable();
    }
}
onCreate({ layer }) {
    this.props.dispatch(submitAreaOfInterest(layer.toGeoJSON()));
}

submitAreaOfInterest验证绘制的形状是否在宾夕法尼亚州的边界内,格式化 API 请求的形状,然后通过 axios 发出请求,同时更新 Redux 状态以指示请求已启动、失败或成功:

export function submitAreaOfInterest(aoi) {
    cancelPriorRequest();
    return (dispatch, getState) => {
        dispatch(startSubmitAreaOfInterest(aoi));
        const { geometry: paGeom } = pennsylvaniaBoundaries.features[0];
        if (!turfContains(paGeom, aoi)) {
            const errorMessage = 'Drawn shape must be within Pennsylvania';
            return dispatch(failSubmitAreaOfInterest(errorMessage));
        }
        const { appPage: { selectedApiEndpoint } } = getState();
        axios.post(<code data-enlighter-language="raw" class="EnlighterJSRAW">${apiServerURL}${selectedApiEndpoint}</code>,
            JSON.stringify({ geometry: JSON.stringify(aoi.geometry) }),
            {
                headers: {
                    'Content-Type': 'application/json',
                },
                cancelToken: new CancelToken((c) => { cancelAxiosRequest = c; }),
            })
             .then(({ data }) => dispatch(completeSubmitAreaOfInterest(data)))
             .catch(() => dispatch(failSubmitAreaOfInterest('API error')));
        return null;
    };
}

一旦数据成功返回,它就会被添加到 Redux 状态,然后通过应用程序传播到 NLCDChart component ,该组件使用 Victory 格式化数据以进行显示。

结果如下:

了解有关 GeoTrellis 的更多信息

现在你可以运行和导航应用程序,了解有关 GeoTrellis、React 和 Leaflet 的更多信息进行开发。例如,你可以:

  • 摄取不同的 GeoTIFF 作为源栅格图层
  • 重写地理处理 API 以执行平均值、最大值或其他地图代数运算map algebra
  • 更新地理处理 API 以并行处理多个栅格图层
  • 渲染视觉图块render visual tiles
  • 调整服务器以替代格式返回结果
  • 在不同类型的图表different kind of chart中显示结果

如果您想了解有关在 Web 应用程序中使用 GeoTrellis 的更多信息,请查看以下资源:

GeoTrellis landsat 教程项目GeoTrellis landsat tutorial project

本教程将介绍如何使用 GeoTrellis 将 landsat 图像处理为 NDVI 切片,并将它们作为 PNG 切片用于在 Leaflet 地图上显示。

使用 GeoTrellis、Lambda 和 API Gateway 提供切片Serving Tiles with GeoTrellis, Lambda, and API Gateway

这篇博文是“Azavea 10% 的研究时间“( Azavea’s 10% research time)的另一个成果,它展现了使用 AWS Lambda 通过 GeoTrellis 创建无服务器切片服务器的可行性。

如何使用 Docker 在 GeoNotebook 中运行 GeoPySparkHow to Run GeoPySpark in a GeoNotebook with Docker

为了便于不熟悉 Scala 的开发人员更了解GeoTrellis,GeoTrellis 最近添加了名为 GeoPySpark 的 Python 绑定( added Python bindings called GeoPySpark)。这篇博文演示了 GIS 开发人员如何使用 Docker 和 Jupyter notebook 快速设置交互式 GeoTrellis 环境。

如何在地理空间应用程序中构建异步工作流How to Build Asynchronous Workflows in a Geospatial Application

这篇博文解释了如何使用 Celery在 Web 应用程序的请求-响应周期内编排复杂的(并且可能需要长时间运行)地理处理操作。

                              **[Kelly Innes(https://www.azavea.com/blog/author/kinnes/)**凯利·英尼斯

                                             Former Software Engineer前软件工程师

03-05 14:38