Merge pull request #23 from GlobalFishingWatch/22-dateline

Handle fishing events crossing the dateline

CHANGES.md

@@ -18,6 +18,10 @@ Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to
 * [#18](https://github.com/GlobalFishingWatch/pipe-events/issues/18): Include
   the `encountered_vessel_id` in the `event_info` for encounters. 
 
+* [#22](https://github.com/GlobalFishingWatch/pipe-events/issues/22): Fix a bug in 
+  fishing events where events crossing the dateline had the wrong min/max/avg 
+  for longitude
+
 ## v0.3.1 2018-11-29
 
 ### Changed

assets/fishing-events.sql.j2

@@ -1,5 +1,8 @@
 #standardSQL
 
+# Include some utility functions
+{% include 'util.sql.j2' %}
+
 #
 # Fishing Events
 #
@@ -30,6 +33,7 @@ WITH
     SELECT
       seg_id,
       timestamp,
+      ST_GEOGPOINT(lon, lat) as point,
       lat,
       lon,
       ifnull(nnet_score,
@@ -41,8 +45,6 @@ WITH
       AND '{{ end_yyyymmdd }}'
       AND lat > -90
       AND lat < 90
-      AND lon > -180
-      AND lon < 180
   ),
 
   #
@@ -179,14 +181,15 @@ WITH
   SELECT
     vessel_id,
     seg_id,
-    AVG(lat) AS lat_mean,
-    AVG(lon) AS lon_mean,
+    ST_CENTROID(ST_UNION_AGG(point)) as centroid,  # compute centoid of all the lat/lon pairs in the event
     event_start,
     MAX(timestamp) as event_end,
     MIN(lat) AS lat_min,
     MAX(lat) AS lat_max,
     MIN(lon) AS lon_min,
     MAX(lon) AS lon_max,
+    MIN(anti_lon(lon)) as anti_lon_min,   # Also get min/max for the anti_longitude (180 degrees opposite) to deal wiht dateline crossing
+    MAX(anti_lon(lon)) as anti_lon_max,
     COUNT(*) AS message_count,
     STRING_AGG(CONCAT(CAST(lon AS string), ' ', CAST(lat AS string)), ', '
       ORDER BY timestamp) AS points_wkt
@@ -202,6 +205,28 @@ WITH
     event_start
   HAVING
     TIMESTAMP_DIFF(event_end, event_start, SECOND) > {{ min_duration }}
+  ),
+
+  # Correct lon_min and lon_max for crossing the dateline (anti-meridian)
+  # And extract lat and lon from the centriod
+  #
+  fishing_event_dateline as (
+    SELECT
+      * except (centroid, lon_min, lon_max, anti_lon_min, anti_lon_max),
+
+      # Get the lat and lon from the computed centroid
+      geopoint_to_struct(centroid).lat as lat_mean,
+      geopoint_to_struct(centroid).lon as lon_mean,
+
+      # determine which direction around the globe is shorter - across the equator (eg -1.0 to 1.0), across the
+      # dateline (eg -179.0 to 179.0) or neither (eg 10.0 to 12.0).  Use this to select which values to use for
+      # min and max longitude
+
+      IF ( (lon_max - lon_min) <= (anti_lon_max - anti_lon_min), lon_min, anti_lon(anti_lon_max) ) as lon_min,
+      IF ( (lon_max - lon_min) <= (anti_lon_max - anti_lon_min), lon_max, anti_lon(anti_lon_min) ) as lon_max
+
+    FROM
+      fishing_event
   )
 
 
@@ -223,4 +248,4 @@ SELECT
   TO_JSON_STRING(STRUCT( message_count )) AS event_info,
   ST_GEOGFROMTEXT(CONCAT( 'MULTIPOINT', " (", points_wkt, ')')) AS event_geography
 FROM
-  fishing_event
+  fishing_event_dateline

assets/util.sql.j2

@@ -0,0 +1,38 @@
+#standardSQL
+
+# Utility Functions
+
+
+# Get the longitude of the antipode for a given longitude
+#
+# This funtion transforms the longitude to the meridian on the opposite side of the world
+# This is useful for doing operations with data that spans the anti-meridian
+#
+# Apply the function again to the transformed value to return to the original meridian
+#
+# Example
+#   anti_lon(0.0) = 0.0
+#   anti_lon(90.0) = -90.0
+#   anti_lon(-90.0) = 90.0
+#   anti_lon(-179.0) = 1.0
+#   anti_lon(1.0) = -179.0
+#
+CREATE TEMPORARY FUNCTION anti_lon ( lon FLOAT64 )
+AS (
+  IF (lon < 0, 180.0 + lon, (180 - lon) * -1)
+);
+
+
+# Convert a lat,lon pont in a GEOGRAPHY, such as created with ST_GEOGPOINT(lon, lat)
+# Retrns a STRUCT(lon, lat)
+#
+# It seems ridiculous that we have to convert to json and then parse it to do this, but bigquery
+# does not provide any other way to get the lat/lon out of a GEOGRAPHY
+
+CREATE TEMPORARY FUNCTION geopoint_to_struct (pt GEOGRAPHY)
+AS (
+  STRUCT(
+    CAST(JSON_EXTRACT_SCALAR(ST_ASGEOJSON(pt), "$['coordinates'][0]") AS FLOAT64) as lon,
+    CAST(JSON_EXTRACT_SCALAR(ST_ASGEOJSON(pt), "$['coordinates'][1]") AS FLOAT64) as lat
+  )
+);