Merge branch 'release/4.0.0'

CHANGELOG.md

@@ -1,3 +1,21 @@
+## 4.0.0
+
+- Add a new `/multi` endpoint that can take a set of shapes,
+  a MultiLine string of streamLines, and a set of operations,
+  to return a result containing all the operations performed
+  for all the shapes. This is useful for the new sub-basin
+  modeling, in which a number of adjacent shapes have to be
+  processed together. Since the most expensive part of this
+  process is fetching the tiles from S3, by fetching them
+  once for the entire set of shapes and reusing the fetched
+  tiles, we reduce the time taken to process a large number
+  of shapes by almost an order of magnitude.
+
+  This new endpoint supports `RasterGroupedCount`,
+  `RasterAverage`, `RasterGroupedAverage`, and
+  `RasterLinesJoin`. It does not support `RasterSummary`
+  because its output type is different.
+
 ## 3.1.0
 
 - Add RasterSummary operation that, given a shape and a list

api/src/main/resources/application.conf

@@ -5,6 +5,9 @@ geoprocessing {
 }
 
 akka.http {
+  server {
+    request-timeout = 59 s
+  }
   parsing {
     max-content-length = 50m
   }

api/src/main/scala/ErrorHandler.scala

@@ -9,6 +9,7 @@ sealed trait GeoprocessingException extends Exception
 case class MissingTargetRasterException() extends GeoprocessingException()
 case class MissingVectorCRSException() extends GeoprocessingException()
 case class MissingVectorException() extends GeoprocessingException()
+case class MissingStreamLinesException() extends GeoprocessingException()
 case class InvalidOperationException(val message: String) extends GeoprocessingException()
 case class UnknownCRSException(val crs: String) extends GeoprocessingException()
 
@@ -30,6 +31,10 @@ trait ErrorHandler {
       println(s"Input error: Missing required vector")
       complete(HttpResponse(BadRequest, entity = "Missing required vector"))
     }
+    case MissingStreamLinesException() => {
+      println(s"Input error: Missing required streamLines")
+      complete(HttpResponse(BadRequest, entity = "Missing required streamLines"))
+    }
     case UnknownCRSException(crs) => {
       println(s"Unknown CRS error: $crs")
       complete(HttpResponse(BadRequest, entity = "Unknown CRS"))

api/src/main/scala/Geoprocessing.scala

@@ -12,7 +12,83 @@ import geotrellis.vector._
 
 import geotrellis.spark._
 
+import cats.implicits._
+
 trait Geoprocessing extends Utils {
+
+  @throws(classOf[MissingStreamLinesException])
+  @throws(classOf[MissingTargetRasterException])
+  def getMultiOperations(input: MultiInput): Future[Map[HucID, Map[OperationID, Map[String, Double]]]] = {
+
+    val hucs: Map[HucID, MultiPolygon] =
+      input.shapes.map(huc => huc.id -> normalizeHuc(huc)).toMap
+
+    val aoi = hucs.values.unionGeometries.asMultiPolygon.get
+
+    val rasterIds =
+      (input.operations.flatMap(_.targetRaster) ++ input.operations.flatMap(_.rasters)).distinct
+
+    val cachedOpts: Map[Operation, Rasterizer.Options] =
+      input.operations.map(op => op -> getRasterizerOptions(op.pixelIsArea)).toMap
+
+    val futureRasters: Map[RasterID, Future[RasterLayer]] =
+      rasterIds.map { rid =>
+        rid -> Future(cropSingleRasterToAOI(rid, 0, aoi))
+      }.toMap
+
+    val tabular: Future[List[(HucID, OperationID, Map[String,Double])]] = Future.sequence {
+      (input.shapes, input.operations).mapN { case (huc, op) =>
+        val shape = hucs(huc.id)
+        val opts = cachedOpts(op)
+
+        val futureLayers: Future[List[RasterLayer]] =
+          op.rasters.map(futureRasters(_)).sequence
+
+        val futureTargetLayer: Future[Option[RasterLayer]] =
+          op.targetRaster.map(futureRasters(_)).sequence
+
+        for {
+          layers <- futureLayers
+          targetLayer <- futureTargetLayer
+        } yield {
+          val results = op.name match {
+            case "RasterGroupedCount" =>
+              rasterGroupedCount(layers, shape, opts).mapValues(_.toDouble)
+
+            case "RasterGroupedAverage" =>
+              targetLayer match {
+                case Some(tl) =>
+                  if (layers.isEmpty) rasterAverage(tl, shape, opts)
+                  else rasterGroupedAverage(layers, tl, shape, opts)
+                case None =>
+                  throw new MissingTargetRasterException
+              }
+
+            case "RasterLinesJoin" =>
+              input.streamLines match {
+                case Some(mls) => {
+                  val lines = (parseMultiLineString(mls) & shape).asMultiLine.get
+                  rasterLinesJoin(layers, Seq(lines)).mapValues(_.toDouble)
+                }
+                case None =>
+                  throw new MissingStreamLinesException
+              }
+          }
+
+          (huc.id, op.label, results)
+        }
+      }
+    }
+
+    val nested: Future[Map[HucID, Map[OperationID, Map[String, Double]]]] = tabular.map { list =>
+      list.groupBy { case (a, _, _) => a }.mapValues {
+        grouped => grouped.map {case (_, b, c) => (b, c) }.toMap
+      }
+    }
+
+    nested
+  }
+
   /**
     * For an InputData object, return a histogram of raster grouped count results.
     *

api/src/main/scala/Utils.scala

@@ -74,6 +74,19 @@ trait Utils {
       .get
   }
 
+  /**
+    * Converts a HUC shape into a reprojected and normalized MultiPolygon
+    * Assumes input is in LatLng and rasters are in ConusAlbers
+    * @param   huc    HUC containing shape string
+    * @return         A MultiPolygon
+    */
+  def normalizeHuc(huc: HUC): MultiPolygon = {
+    parseGeometry(huc.shape, LatLng, ConusAlbers)
+      .buffer(0)
+      .asMultiPolygon
+      .get
+  }
+
   /**
     * Given an optional boolean value, if it is true, returns Rasterizer
     * Options treating the raster pixels as areas. If false, returns Options
@@ -148,6 +161,12 @@ trait Utils {
     }
   }
 
+  /**
+    * Convenience flavor of the above with defaults
+    */
+  def parseMultiLineString(geoJson: String): MultiLine =
+    parseMultiLineString(geoJson, LatLng, ConusAlbers)
+
   /**
     * Given a sequence of MultiLines and an area of interest, crops the lines
     * to the area of interest and returns a sequence containing the cropped lines.

api/src/main/scala/WebServer.scala

@@ -26,12 +26,37 @@ case class ResultInt(result: Map[String, Int])
 case class ResultDouble(result: Map[String, Double])
 case class ResultSummary(result: Seq[Map[String, Double]])
 
+// HUCs have an id and a shape. The shape is GeoJSON, but we've transmitted
+// them as Strings in the past so we continue to do so here.
+case class HUC (
+  id: HucID,
+  shape: GeoJSONString // GeoJSON Polygon or MultiPolygon
+)
+
+case class Operation (
+  name: String, // RasterGroupedCount, RasterGroupedAverage, RasterLinesJoin
+  label: OperationID,
+  rasters: List[RasterID],
+  targetRaster: Option[RasterID],
+  pixelIsArea: Option[Boolean]
+)
+
+case class MultiInput (
+  shapes: List[HUC],
+  streamLines: Option[GeoJSONString], // GeoJSON MultiLineString
+  operations: List[Operation]
+)
+
 object PostRequestProtocol extends DefaultJsonProtocol {
   implicit val inputFormat = jsonFormat10(InputData)
   implicit val postFormat = jsonFormat1(PostRequest)
   implicit val resultFormat = jsonFormat1(ResultInt)
   implicit val resultDoubleFormat = jsonFormat1(ResultDouble)
   implicit val resultSummaryFormat = jsonFormat1(ResultSummary)
+
+  implicit val hucFormat = jsonFormat2(HUC)
+  implicit val operationFormat = jsonFormat5(Operation)
+  implicit val multiInputFormat = jsonFormat3(MultiInput)
 }
 
 object WebServer extends HttpApp with App with LazyLogging with Geoprocessing with ErrorHandler {
@@ -63,6 +88,11 @@ object WebServer extends HttpApp with App with LazyLogging with Geoprocessing wi
               }
             }
           }
+        } ~
+        path("multi") {
+          entity(as[MultiInput]) { input =>
+            complete(getMultiOperations(input))
+          }
         }
       }
     }

api/src/main/scala/package.scala

@@ -0,0 +1,11 @@
+package org.wikiwatershed.mmw
+
+import geotrellis.spark._
+
+package object geoprocessing {
+  type HucID = String
+  type OperationID = String
+  type RasterID = String
+  type RasterLayer = TileLayerCollection[SpatialKey]
+  type GeoJSONString = String
+}

examples/MapshedJob_DRB_Medium.json

@@ -1862,7 +1862,7 @@
       "{\"type\":\"LineString\",\"coordinates\":[[-75.2233796296351,39.9095833333373],[-75.2234722222277,39.9096759259299],[-75.2235648148203,39.9097685185225],[-75.2236574074128,39.9098611111151],[-75.2237500000055,39.9099537037077],[-75.223842592598,39.9100462963003],[-75.2239351851906,39.9100462963003],[-75.2240277777832,39.9100462963003],[-75.2241203703758,39.9100462963003],[-75.2242129629684,39.9100462963003],[-75.224305555561,39.9100462963003],[-75.2243981481536,39.9101388888929],[-75.2244907407462,39.9101388888929],[-75.2245833333388,39.9102314814855]]}"
     ],
     "rasters": [
-      "nlcd-2011-30m-epsg5070-0.10.0"
+      "nlcd-2011-30m-epsg5070-512-int8"
     ],
     "rasterCRS": "ConusAlbers",
     "polygonCRS": "LatLng",

examples/MapshedJob_DRB_Small.json

@@ -822,10 +822,10 @@
     ],
     "vectorCRS": "LatLng",
     "rasters": [
-      "nlcd-2011-30m-epsg5070-0.10.0"
+      "nlcd-2011-30m-epsg5070-512-int8"
     ],
     "rasterCRS": "ConusAlbers",
     "operationType": "RasterLinesJoin",
     "zoom": 0
   }
-}
+}

examples/MapshedJob_NHD.json

@@ -18,10 +18,10 @@
     ],
     "vectorCRS": "LatLng",
     "rasters": [
-      "nlcd-2011-30m-epsg5070-0.10.0"
+      "nlcd-2011-30m-epsg5070-512-int8"
     ],
     "rasterCRS": "ConusAlbers",
     "operationType": "RasterLinesJoin",
     "zoom": 0
   }
-}
+}

examples/MapshedJob_RasterAverage.json

@@ -3,7 +3,7 @@
     "polygon": [
       "{\"type\":\"MultiPolygon\",\"coordinates\":[[[[-75.1626205444336,39.95580659996906],[-75.25531768798828,39.94514735903112],[-75.22785186767578,39.89446035777916],[-75.1461410522461,39.88761144548104],[-75.09309768676758,39.91078961774283],[-75.09464263916016,39.93817189499188],[-75.12039184570312,39.94435771955196],[-75.1626205444336,39.95580659996906]]]]}"
     ],
-    "targetRaster": "us-ssugro-kfactor-30m-epsg5070",
+    "targetRaster": "us-ssugro-kfactor-30m-epsg5070-512",
     "rasters": [],
     "rasterCRS": "ConusAlbers",
     "polygonCRS": "LatLng",

examples/MapshedJob_RasterGroupedAverage.json

@@ -3,9 +3,9 @@
     "polygon": [
       "{\"type\":\"MultiPolygon\",\"coordinates\":[[[[-75.1626205444336,39.95580659996906],[-75.25531768798828,39.94514735903112],[-75.22785186767578,39.89446035777916],[-75.1461410522461,39.88761144548104],[-75.09309768676758,39.91078961774283],[-75.09464263916016,39.93817189499188],[-75.12039184570312,39.94435771955196],[-75.1626205444336,39.95580659996906]]]]}"
     ],
-    "targetRaster": "us-ssugro-kfactor-30m-epsg5070",
+    "targetRaster": "us-ssugro-kfactor-30m-epsg5070-512",
     "rasters": [
-      "nlcd-2011-30m-epsg5070-0.10.0"
+      "nlcd-2011-30m-epsg5070-512-int8"
     ],
     "rasterCRS": "ConusAlbers",
     "polygonCRS": "LatLng",

examples/MapshedJob_RasterGroupedCount.json

@@ -4,12 +4,12 @@
       "{\"type\":\"MultiPolygon\",\"coordinates\":[[[[-75.1626205444336,39.95580659996906],[-75.25531768798828,39.94514735903112],[-75.22785186767578,39.89446035777916],[-75.1461410522461,39.88761144548104],[-75.09309768676758,39.91078961774283],[-75.09464263916016,39.93817189499188],[-75.12039184570312,39.94435771955196],[-75.1626205444336,39.95580659996906]]]]}"
     ],
     "rasters": [
-      "nlcd-2011-30m-epsg5070-0.10.0",
-      "ssurgo-hydro-groups-30m-epsg5070-0.10.0"
+      "nlcd-2011-30m-epsg5070-512-int8",
+      "ssurgo-hydro-groups-30m-epsg5070-512-int8"
     ],
     "rasterCRS": "ConusAlbers",
     "polygonCRS": "LatLng",
     "operationType": "RasterGroupedCount",
     "zoom": 0
   }
-}
+}

examples/MapshedJob_RasterGroupedSum.json

@@ -3,9 +3,9 @@
     "polygon": [
       "{\"type\":\"MultiPolygon\",\"coordinates\":[[[[-75.1626205444336,39.95580659996906],[-75.25531768798828,39.94514735903112],[-75.22785186767578,39.89446035777916],[-75.1461410522461,39.88761144548104],[-75.09309768676758,39.91078961774283],[-75.09464263916016,39.93817189499188],[-75.12039184570312,39.94435771955196],[-75.1626205444336,39.95580659996906]]]]}"
     ],
-    "targetRaster": "us-ssugro-kfactor-30m-epsg5070",
+    "targetRaster": "us-ssugro-kfactor-30m-epsg5070-512",
     "rasters": [
-      "nlcd-2011-30m-epsg5070-0.10.0"
+      "nlcd-2011-30m-epsg5070-512-int8"
     ],
     "rasterCRS": "ConusAlbers",
     "polygonCRS": "LatLng",

examples/MapshedJob_RasterSummary.json

@@ -4,7 +4,7 @@
       "{\"type\":\"MultiPolygon\",\"coordinates\":[[[[-75.1626205444336,39.95580659996906],[-75.25531768798828,39.94514735903112],[-75.22785186767578,39.89446035777916],[-75.1461410522461,39.88761144548104],[-75.09309768676758,39.91078961774283],[-75.09464263916016,39.93817189499188],[-75.12039184570312,39.94435771955196],[-75.1626205444336,39.95580659996906]]]]}"
     ],
     "rasters": [
-      "ned-nhdplus-30m-epsg5070",
+      "ned-nhdplus-30m-epsg5070-512",
       "us-percent-slope-30m-epsg5070-512"
     ],
     "rasterCRS": "ConusAlbers",