Merge pull request #12 from pacificclimate/projections

handle projections mismatch between database and front end

README.md

@@ -1,5 +1,5 @@
 # scip-backend
-This application accesses a PostGIS database containing information on predefined areas a user might care about (watersheds, basins, etc) and the locations of various salmon populations. Its role is to provide data to the [SCIP frontend](https://github.com/pacificclimate/scip-frontend).
+This application accesses a PostGIS database containing information on predefined areas a user might care about (watersheds, basins, and salmon conservation units) and metadata about salmon populations. Its role is to provide data to the [SCIP frontend](https://github.com/pacificclimate/scip-frontend).
 
 ## Installation for development
 

pyproject.toml

@@ -17,6 +17,8 @@ flask-sqlalchemy = "^3.0.3"
 sqlalchemy-sqlschema = "^0.1"
 black = "^23.3.0"
 sphinx = "^7.0.1"
+shapely = "^2.0.2"
+geojson = "^3.1.0"
 
 [tool.poetry.group.test.dependencies]
 pytest = "^7.3.1"

scip/api/population.py

@@ -3,6 +3,7 @@
 from sqlalchemy_sqlschema.sql import get_schema
 from sqlalchemy import func
 from scip.api.validators import parse_common_name, parse_subgroup
+from scip.api.projection import to_4326, assume_4326
 
 
 def population(
@@ -43,15 +44,17 @@ def population(
                 Taxon.subgroup.label("subgroup"),
                 ConservationUnit.name.label("name"),
                 ConservationUnit.code.label("code"),
-                func.ST_AsGeoJSON(ConservationUnit.boundary).label("boundary"),
-                func.ST_ASGeoJSON(ConservationUnit.outlet).label("outlet"),
+                func.ST_AsGeoJSON(to_4326(ConservationUnit.boundary)).label("boundary"),
+                func.ST_ASGeoJSON(to_4326(ConservationUnit.outlet)).label("outlet"),
             )
             .join(Population, Population.conservation_unit_id == ConservationUnit.id)
             .join(Taxon, Population.taxon_id == Taxon.id)
         )
 
         if overlap:
-            q = q.filter(ConservationUnit.boundary.ST_Intersects(overlap))
+            q = q.filter(
+                to_4326(ConservationUnit.boundary).ST_Intersects(assume_4326(overlap))
+            )
 
         if name:
             q = q.filter(ConservationUnit.name == name)

scip/api/projection.py

@@ -0,0 +1,34 @@
+# The postGIS database is in the BC Albers projection. The front end uses no
+# particular projection, but we assume it to be EPSG 4326 for convenience
+# functions in this file convert between the two systems.
+
+# use assume_4326 on parameters received from the front end, to add a projection
+# use to_4326 on data fetched from the back end, to convert it to EPSG 4236
+
+# All geoJSON is officially defined as 4326 as of 2016.
+# See https://www.rfc-editor.org/rfc/rfc7946#section-4
+# So assume_4326 should theoretically just be able to convert
+# any WKT string to geoJSON and have that be sufficient to establish
+# the geometry as 4326. However, we got errors doing this, so we are
+# still using the deprecated "crs" attribute with our geoJSON strings.
+# TODO: determine source of those errors, switch to modern geoJSON
+# handling.
+
+from sqlalchemy import func
+import shapely.wkt
+import geojson
+import json
+
+
+def to_4326(geom):
+    """convert a geometry to 4326 to send to the front end"""
+    return func.ST_Transform(geom, 4326)
+
+
+def assume_4326(wkt):
+    """convert a WKT string (no projection) to an equivalent geoJSON string in WPSG 4326"""
+    shape = shapely.wkt.loads(wkt)
+    gj = json.loads(geojson.dumps(shape))
+    gj["crs"] = {"type": "name", "properties": {"name": "epsg:4326"}}
+
+    return json.dumps(gj)

scip/api/region_helpers.py

@@ -17,6 +17,7 @@
 
 from salmon_occurrence import Region, ConservationUnit, Population, Taxon
 from sqlalchemy import func
+from scip.api.projection import to_4326, assume_4326
 
 
 def cu_with_taxon(session, common_name, subgroup):
@@ -44,8 +45,8 @@ def cu_geometry_only(session):
     q = session.query(
         ConservationUnit.name.label("name"),
         ConservationUnit.code.label("code"),
-        func.ST_AsGeoJSON(ConservationUnit.boundary).label("boundary"),
-        func.ST_ASGeoJSON(ConservationUnit.outlet).label("outlet"),
+        func.ST_AsGeoJSON(to_4326(ConservationUnit.boundary)).label("boundary"),
+        func.ST_ASGeoJSON(to_4326(ConservationUnit.outlet)).label("outlet"),
     )
 
     return q
@@ -62,7 +63,9 @@ def build_cu_query(
         q = cu_geometry_only(session)
 
     if overlap:
-        q = q.filter(ConservationUnit.boundary.ST_Intersects(overlap))
+        q = q.filter(
+            to_4326(ConservationUnit.boundary).ST_Intersects(assume_4326(overlap))
+        )
 
     if name:
         q = q.filter(ConservationUnit.name == name)
@@ -101,8 +104,8 @@ def region_geometry_only(session, kind):
         Region.name.label("name"),
         Region.code.label("code"),
         Region.kind.label("kind"),
-        func.ST_AsGeoJSON(Region.boundary).label("boundary"),
-        func.ST_AsGeoJSON(Region.outlet).label("outlet"),
+        func.ST_AsGeoJSON(to_4326(Region.boundary)).label("boundary"),
+        func.ST_AsGeoJSON(to_4326(Region.outlet)).label("outlet"),
     ).filter(Region.kind == kind)
 
     return q
@@ -118,7 +121,7 @@ def build_region_query(
         q = region_geometry_only(session, kind)
 
     if overlap:
-        q = q.filter(Region.boundary.ST_Intersects(overlap))
+        q = q.filter(to_4326(Region.boundary).ST_Intersects(assume_4326(overlap)))
 
     if name:
         q = q.filter(Region.name == name)