Analyzing NYC Taxis with GeoMesa

%md ## 1. Spatial Search Problems

Spatial data has two fundamental query types:
<p></p>
* K nearest neighbors
* Range queries

%md

Retrieve the k closest points to a given query point

<img src="files/notebooks/images/nyc-taxi/knn.png" width="50%"></img>

%md

e.g. find all points inside a polygon or circle

<img src="files/notebooks/images/nyc-taxi/rangeq.png" width="50%"></img>

%md ## 2. Spatial Indexing

Efficient spatial queries at scale requires putting the geometries into a **spatial index**.

Querying indexed data with geometric predicate pushdown can dramatically improve performance

**R Trees** and **Space Filling Curves** are two common spatial indexing techniques

Almost all spatial data structures share the same principle to enable efficient search: branch and bound. It means arranging data in a tree-like structure that allows discarding branches at once if they do not fit our search criteria (data skipping).

%md

R-trees are tree data structures used for spatial access methods, i.e., for indexing multi-dimensional information such as geographical coordinates, rectangles or polygons.  The key idea of the data structure is to **group nearby objects and represent them with their minimum bounding rectangle** in the next higher level of the tree; the "R" in R-tree is for rectangle. Since all objects lie within this bounding rectangle, a query that does not intersect the bounding rectangle also cannot intersect any of the contained objects.

<img src="files/derek/geospark/rtree.png" width="50%"></img>

%md

Space Filling Curves are technique to colocate related information in the same set of files

i.e. they allow you to **map multidimensional data to one dimension while maintaining locality**

Mapping things down to 1 dimension allows for concepts like range queries, which essentially cuts out scanning entire sections of data since we preserve the locality of the data in higher dimensions with the curve.  

<img src="files/notebooks/images/nyc-taxi/sfc.gif" width="40%"></img>

%md ## 3. GeoMesa: Geospatial Analytics 

There are a number of purpose-built libraries which extend Apache Spark for geospatial analytics. [GeoMesa](https://github.com/locationtech/geomesa) is a a distributed framework especially adept at handling vector data. Others used by Databricks customers include [GeoSpark](http://geospark.datasyslab.org/), [GeoTrellis](https://geotrellis.io/), and [Rasterframes](https://rasterframes.io/). While they often come with a bit of a learning curve, these frameworks often offer multiple language bindings and have much better scaling and performance than non-formalized approaches.  

%md _Set serializers for GeoMesa._

%md _Add GeoMesa User Defined Types (UDT) and Functions (UDF)._

%md ## 4. Load NYC Taxi Data from Delta Table

%md _Since we are dealing with 200M points and some skew, want to increase parallelism._

%md 
* _Get a DataFrame reference to the `nyctaxidelta` table and convert the pickup and dropoff string coordinates into Geometry types using GeoMesa._
* _Also, NYC data should be ~-75 LON and ~40 LAT, so we want to filter out bad points which are present in the provided data._

%md _Let's consider the skew for pickups, notice ~110M pickups are in a single geohash (Manhattan)._

%md ## 5. Load & View NYC Neighborhood Polygon Data (WKT)
<small></small>
* Data from https://data.cityofnewyork.us/Transportation/NYC-Taxi-Zones/d3c5-ddgc, exported as CSV which has the WKT.
* Data was placed in DBFS at `dbfs:/ml/blogs/geospatial/nyc_taxi_zones.wkt.csv`

%md ## 6. Let's execute a spatial join at scale!

We want to associate each pickup (point) with the appropriate NYC taxi zone (polygon).

This _"if point is within polygon"_ query predicate would require comparision of >200M points to ~250 polygons in our example.  (i.e. **worst case 50,000,000,000 expensive comparisons**)

To constain joins we are leveraging the precomputed **geohash** information to significantly "prune" the solution space.

We then evaluate the **geospatial predicate** `st_contains($"polygon", $"pickupPoint")` (to filter out false positives) 

%md _First we figure out `geohash` values with more than 1M pickups.  (used as a "skew hint" to optimize performance below)_

%md _Cache joined for further use in the notebook (still have the ~200M Points)._