Skip to content

Commit

Permalink
Merge pull request #505 from jbouffard/feature/point-value
Browse files Browse the repository at this point in the history 
Point Value Lookups
  • Loading branch information
Jacob Bouffard committed Oct 9, 2017
2 parents fd33978 + c0e9c69 commit c8a7cda
Show file tree
Hide file tree
Showing 4 changed files with 529 additions and 4 deletions.
92 changes: 91 additions & 1 deletion ...ark-backend/geotrellis/src/main/scala/geopyspark/geotrellis/SpatialTiledRasterLayer.scala
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ import geotrellis.raster.io.geotiff._
import geotrellis.raster.io.geotiff.compression._
import geotrellis.raster.rasterize._
import geotrellis.raster.render._
import geotrellis.raster.resample.ResampleMethod
import geotrellis.raster.resample.{ResampleMethod, PointResampleMethod, Resample}
import geotrellis.spark._
import geotrellis.spark.costdistance.IterativeCostDistance
import geotrellis.spark.io._
Expand Down Expand Up @@ -46,6 +46,8 @@ import org.apache.spark.SparkContext._
import java.util.ArrayList
import scala.reflect._
import scala.collection.JavaConverters._
import scala.collection.mutable.ArrayBuffer


class SpatialTiledRasterLayer(
val zoomLevel: Option[Int],
Expand Down Expand Up @@ -313,6 +315,94 @@ class SpatialTiledRasterLayer(

def collectKeys(): java.util.ArrayList[Array[Byte]] =
PythonTranslator.toPython[SpatialKey, ProtoSpatialKey](rdd.keys.collect)

def getPointValues(
points: java.util.Map[Long, Array[Byte]],
resampleMethod: PointResampleMethod
): java.util.Map[Long, Array[Double]] = {
val mapTrans = rdd.metadata.layout.mapTransform

val idedKeys: Map[Long, Point] =
points
.asScala
.mapValues { WKB.read(_).asInstanceOf[Point] }
.toMap

val pointKeys =
idedKeys
.foldLeft(Map[SpatialKey, Array[(Long, Point)]]()) {
case (acc, elem) =>
val pointKey = mapTrans(elem._2)

acc.get(pointKey) match {
case Some(arr) => acc + (pointKey -> (elem +: arr))
case None => acc + (pointKey -> Array(elem))
}
}

val matchedKeys =
resampleMethod match {
case r: PointResampleMethod => _getPointValues(pointKeys, mapTrans, r)
case _ => _getPointValues(pointKeys, mapTrans)
}

val remainingKeys = idedKeys.keySet diff matchedKeys.keySet

if (remainingKeys.isEmpty)
matchedKeys.asJava
else
remainingKeys.foldLeft(matchedKeys){ case (acc, elem) =>
acc + (elem -> null)
}.asJava
}

def _getPointValues(
pointKeys: Map[SpatialKey, Array[(Long, Point)]],
mapTrans: MapKeyTransform,
resampleMethod: PointResampleMethod
): Map[Long, Array[Double]] = {
val resamplePoint = (tile: Tile, extent: Extent, point: Point) => {
Resample(resampleMethod, tile, extent).resampleDouble(point)
}

rdd.flatMap { case (k, v) =>
pointKeys.get(k) match {
case Some(arr) =>
val keyExtent = mapTrans(k)
val rasterExtent = RasterExtent(keyExtent, v.cols, v.rows)

arr.map { case (id, point) =>
(id, v.bands.map { resamplePoint(_, keyExtent, point) } toArray)
}
case None => Seq()
}
}.collect().toMap
}

def _getPointValues(
pointKeys: Map[SpatialKey, Array[(Long, Point)]],
mapTrans: MapKeyTransform
): Map[Long, Array[Double]] =
rdd.flatMap { case (k, v) =>
pointKeys.get(k) match {
case Some(arr) =>
val keyExtent = mapTrans(k)
val rasterExtent = RasterExtent(keyExtent, v.cols, v.rows)

arr.map { case (id, point) =>
val (gridCol, gridRow) = rasterExtent.mapToGrid(point)

val values = Array.ofDim[Double](v.bandCount)

cfor(0)(_ < v.bandCount, _ + 1){ index =>
values(index) = v.band(index).getDouble(gridCol, gridRow)
}

(id, values)
}
case None => Seq()
}
}.collect().toMap
}


Expand Down
93 changes: 91 additions & 2 deletions ...rk-backend/geotrellis/src/main/scala/geopyspark/geotrellis/TemporalTiledRasterLayer.scala
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ import geotrellis.raster.io.geotiff._
import geotrellis.raster.io.geotiff.compression._
import geotrellis.raster.rasterize._
import geotrellis.raster.render._
import geotrellis.raster.resample.ResampleMethod
import geotrellis.raster.resample.{ResampleMethod, PointResampleMethod, Resample}
import geotrellis.spark._
import geotrellis.spark.costdistance.IterativeCostDistance
import geotrellis.spark.io._
Expand Down Expand Up @@ -44,10 +44,11 @@ import org.apache.spark.rdd._
import org.apache.spark.SparkContext._

import java.util.ArrayList
import java.time.ZonedDateTime
import java.time.{ZonedDateTime, ZoneId}

import scala.reflect._
import scala.collection.JavaConverters._
import scala.collection.mutable.ArrayBuffer


class TemporalTiledRasterLayer(
Expand Down Expand Up @@ -387,6 +388,94 @@ class TemporalTiledRasterLayer(

def collectKeys(): java.util.ArrayList[Array[Byte]] =
PythonTranslator.toPython[SpaceTimeKey, ProtoSpaceTimeKey](rdd.keys.collect)

def getPointValues(
points: java.util.Map[Long, Array[Byte]],
resampleMethod: PointResampleMethod
): java.util.Map[Long, (Long, Array[Double])] = {
val mapTrans = rdd.metadata.layout.mapTransform

val idedKeys: Map[Long, Point] =
points
.asScala
.mapValues { WKB.read(_).asInstanceOf[Point] }
.toMap

val pointKeys =
idedKeys
.foldLeft(Map[SpatialKey, Array[(Long, Point)]]()) {
case (acc, elem) =>
val pointKey = mapTrans(elem._2)

acc.get(pointKey) match {
case Some(arr) => acc + (pointKey -> (elem +: arr))
case None => acc + (pointKey -> Array(elem))
}
}

val matchedKeys =
resampleMethod match {
case r: PointResampleMethod => _getPointValues(pointKeys, mapTrans, r)
case _ => _getPointValues(pointKeys, mapTrans)
}

val remainingKeys = idedKeys.keySet diff matchedKeys.keySet

if (remainingKeys.isEmpty)
matchedKeys.asJava
else
remainingKeys.foldLeft(matchedKeys){ case (acc, elem) =>
acc + (elem -> null)
}.asJava
}

def _getPointValues(
pointKeys: Map[SpatialKey, Array[(Long, Point)]],
mapTrans: MapKeyTransform,
resampleMethod: PointResampleMethod
): Map[Long, (Long, Array[Double])] = {
val resamplePoint = (tile: Tile, extent: Extent, point: Point) => {
Resample(resampleMethod, tile, extent).resampleDouble(point)
}

rdd.flatMap { case (k, v) =>
pointKeys.get(k.getComponent[SpatialKey]) match {
case Some(arr) =>
val keyExtent = mapTrans(k)
val rasterExtent = RasterExtent(keyExtent, v.cols, v.rows)

arr.map { case (id, point) =>
(id -> (k.instant, v.bands.map { resamplePoint(_, keyExtent, point) } toArray))
}
case None => Seq()
}
}.collect().toMap
}

def _getPointValues(
pointKeys: Map[SpatialKey, Array[(Long, Point)]],
mapTrans: MapKeyTransform
): Map[Long, (Long, Array[Double])] =
rdd.flatMap { case (k, v) =>
pointKeys.get(k.getComponent[SpatialKey]) match {
case Some(arr) =>
val keyExtent = mapTrans(k)
val rasterExtent = RasterExtent(keyExtent, v.cols, v.rows)

arr.map { case (id, point) =>
val (gridCol, gridRow) = rasterExtent.mapToGrid(point)

val values = Array.ofDim[Double](v.bandCount)

cfor(0)(_ < v.bandCount, _ + 1){ index =>
values(index) = v.band(index).getDouble(gridCol, gridRow)
}

(id -> (k.instant, values))
}
case None => Seq()
}
}.collect().toMap
}


Expand Down
136 changes: 135 additions & 1 deletion geopyspark/geotrellis/layer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
import json
import datetime
import shapely.wkb
from shapely.geometry import Polygon, MultiPolygon
from shapely.geometry import Polygon, MultiPolygon, Point
from geopyspark.geotrellis.protobufcodecs import (multibandtile_decoder,
projected_extent_decoder,
temporal_projected_extent_decoder,
Expand Down Expand Up @@ -1480,6 +1480,140 @@ def normalize(self, new_min, new_max, old_min=None, old_max=None):

return TiledRasterLayer(self.layer_type, srdd)

def get_point_values(self, points, resample_method=None):
"""Returns the values of the layer at given points.
Note:
Only points that are contained within a layer will be sampled.
This means that if a point lies on the southern or eastern boundary
of a cell, it will not be sampled.
Args:
points([shapely.geometry.Point] or {k: shapely.geometry.Point}):
Either a list of, or a dictionary whose values are ``shapely.geometry.Point``\s.
If a dictionary, then the type of its keys does not matter.
These points must be in the same projection as the tiles within the layer.
resample_method(str or :class:`~gepyspark.ResampleMethod`, optional): The resampling method
to use before obtaining the point values. If not specified, then ``None`` is used.
Note:
Not all ``ResampleMethod``\s can be used to resample point values.
``ResampleMethod.NEAREST_NEIGHBOR``, ``ResampleMethod.BILINEAR```,
``ResampleMethod.CUBIC_CONVOLUTION``, and ``ResampleMethod.CUBIC_SPLINE``
are the only ones that can be used.
Returns:
The return type will vary depending on the type of ``points`` and the ``layer_type`` of
the sampled layer.
If ``points`` is a ``list`` and the ``layer_type`` is ``SPATIAL``:
[(shapely.geometry.Point, [float])]
If ``points`` is a ``list`` and the ``layer_type`` is ``SPACETIME``:
[(shapely.geometry.Point, datetime.datetime, [float])]
If ``points`` is a ``dict`` and the ``layer_type`` is ``SPATIAL``:
{k: (shapely.geometry.Point, [float])}
If ``points`` is a ``dict`` and the ``layer_type`` is ``SPACETIME``:
{k: (shapely.geometry.Point, datetime.datetime, [float])}
The ``shapely.geometry.Point`` in all of these returns is the original sampled point
given. The ``[float]`` are the sampled values, one for each band. If the ``layer_type``
was ``SPACETIME``, then the timestamp will also be included in the results represented
by a ``datetime.datetime`` instance.
Note:
The sampled values will always be returned as ``float``\s. Regardless of the
``cellType`` of the layer.
If ``points`` was given as a ``dict`` then the keys of that dictionary will be the
keys in the returned ``dict``.
"""

if resample_method:
resample_method = ResampleMethod(resample_method)

point_resampling_methods = [
ResampleMethod.NEAREST_NEIGHBOR,
ResampleMethod.BILINEAR,
ResampleMethod.CUBIC_CONVOLUTION,
ResampleMethod.CUBIC_SPLINE
]

if resample_method not in point_resampling_methods:
raise ValueError(resample_method, "Cannot be used to resample point values")

ided_points = {}
ided_bytes = {}

def to_datetime(instant):
return datetime.datetime.utcfromtimestamp(instant / 1000)

if isinstance(points, list):
list_result = []

types = [type(x) for x in points]

if Point not in types or len(set(types)) > 1:
raise TypeError(set(types), "found, but only shapely Points can sample values")

for point in points:
point_id = id(point)

ided_bytes[point_id] = shapely.wkb.dumps(point)
ided_points[point_id] = point

scala_result = self.srdd.getPointValues(ided_bytes, resample_method)

if self.layer_type == LayerType.SPATIAL:
for key, value in scala_result.items():
list_result.append((ided_points[key], list(value) if value else None))
else:
for key, value in scala_result.items():
list_result.append(
(ided_points[key],
to_datetime(value._1()) if value else None,
list(value._2()) if value else None
)
)

return list_result

elif isinstance(points, dict):
ided_labels = {}
dict_result = {}

types = [type(x) for x in points.values()]

if Point not in types or len(set(types)) > 1:
raise TypeError(set(types), "found, but only shapely Points can sample values")

for key, value in points.items():
point_id = id(value)

ided_points[point_id] = value
ided_bytes[point_id] = shapely.wkb.dumps(value)
ided_labels[point_id] = key

scala_result = self.srdd.getPointValues(ided_bytes, resample_method)

if self.layer_type == LayerType.SPATIAL:
for key, value in scala_result.items():
dict_result[ided_labels[key]] = (ided_points[key], list(value) if value else None)
else:
for key, value in scala_result.items():
dict_result[ided_labels[key]] = (
ided_points[key],
to_datetime(value._1()) if value else None,
list(value._2()) if value else None
)

return dict_result
else:
raise TypeError("Expected a list or dict. Instead got", type(points))


@staticmethod
def _process_polygonal_summary(geometry, operation):
if isinstance(geometry, (Polygon, MultiPolygon)):
Expand Down

0 comments on commit c8a7cda

Please sign in to comment.