Merge pull request #137 from DataIntegrationGroup/jab-location-transfer-updates

BDMS 111: location updates & other

Author: jirhiker

.github/workflows/tests.yml

@@ -5,7 +5,7 @@ name: Tests
 
 on:
   pull_request:
-    branches: [ "main",'pre-production', 'transfer']
+    branches: ['production', 'staging', 'transfer']
 
 permissions:
   contents: read

api/geospatial.py

@@ -14,7 +14,7 @@
 # limitations under the License.
 # ===============================================================================
 import json
-from typing import Annotated, List, Union
+from typing import Annotated, List
 
 from fastapi import APIRouter, Query, HTTPException
 from fastapi.responses import FileResponse
@@ -100,7 +100,9 @@ def get_feature_collection(
 
     things = get_thing_features(session, thing_type, group)
 
-    def make_feature_dict(thing, geometry, *other):
+    def make_feature_dict(thing, geometry, elevation, *other):
+        geometry = json.loads(geometry)
+        geometry["coordinates"].append(elevation)
         return {
             "type": "Feature",
             "properties": {
@@ -109,7 +111,7 @@ def make_feature_dict(thing, geometry, *other):
                 "name": thing.name,
                 "group": group,
             },
-            "geometry": json.loads(geometry),
+            "geometry": geometry,
         }
 
     features = [make_feature_dict(*item) for item in things]

api/location.py

@@ -30,6 +30,7 @@
 from services.geospatial_helper import make_within_wkt
 from services.query_helper import make_query, order_sort_filter, simple_get_by_id
 from services.crud_helper import model_patcher, model_deleter, model_adder
+from services.location_helper import set_geographic_attributes
 
 from fastapi import APIRouter
 
@@ -48,7 +49,9 @@ async def create_location(
     """
     Create a new sample location in the database.
     """
-    return model_adder(session, Location, location_data, user=user)
+    location = model_adder(session, Location, location_data, user=user)
+    set_geographic_attributes(session, location_data, location)
+    return location
 
 
 @router.patch(
@@ -64,7 +67,9 @@ async def update_location(
     """
     Update a sample location in the database.
     """
-    return model_patcher(session, Location, location_id, location_data, user=user)
+    location = model_patcher(session, Location, location_id, location_data, user=user)
+    set_geographic_attributes(session, location_data, location)
+    return location
 
 
 # @router.get("/shapefile", summary="Get location as shapefile")

constants.py

@@ -15,4 +15,5 @@
 # ===============================================================================
 
 SRID_WGS84 = 4326
+SRID_UTM_ZONE_13N = 26913
 # ============= EOF =============================================

db/group.py

@@ -20,6 +20,7 @@
 from sqlalchemy.orm import relationship, Mapped
 from sqlalchemy.testing.schema import mapped_column
 
+from constants import SRID_WGS84
 from db.base import Base, AutoBaseMixin, ReleaseMixin
 
 
@@ -28,7 +29,7 @@ class Group(Base, AutoBaseMixin, ReleaseMixin):
     description: Mapped[str] = mapped_column(String(255), nullable=True)
     name: Mapped[str] = mapped_column(String(100), nullable=False, unique=True)
     project_area: Mapped[Optional[WKBElement]] = mapped_column(
-        Geometry(geometry_type="MULTIPOLYGON", srid=4326, spatial_index=True)
+        Geometry(geometry_type="MULTIPOLYGON", srid=SRID_WGS84, spatial_index=True)
     )
 
     # Foreign Keys

db/location.py

@@ -21,8 +21,6 @@
 from uuid import UUID
 
 from sqlalchemy import (
-    Column,
-    Integer,
     String,
     ForeignKey,
     DateTime,
@@ -32,6 +30,7 @@
 from sqlalchemy.orm import relationship, Mapped, mapped_column
 from sqlalchemy.ext.associationproxy import association_proxy, AssociationProxy
 
+from constants import SRID_WGS84
 from db.base import Base, AutoBaseMixin, ReleaseMixin
 from db.lexicon import lexicon_term
 
@@ -43,9 +42,12 @@ class Location(Base, AutoBaseMixin, ReleaseMixin):
         String(36), nullable=True, unique=True
     )
     description: Mapped[str] = mapped_column
-    name: Mapped[str] = mapped_column(String(255), nullable=True)
+    # name: Mapped[str] = mapped_column(String(255), nullable=True)
     point: Mapped[WKBElement] = mapped_column(
-        Geometry(geometry_type="POINTZ", srid=4326, spatial_index=True)
+        Geometry(geometry_type="POINT", srid=SRID_WGS84, spatial_index=True)
+    )
+    elevation: Mapped[float] = mapped_column(
+        nullable=False, comment="in meters with vertical datum of NAVD88"
     )
 
     state: Mapped[str] = lexicon_term(nullable=True, default="New Mexico")
@@ -65,7 +67,7 @@ class Location(Base, AutoBaseMixin, ReleaseMixin):
     )
 
     # --- Proxy Definitions ---
-    things: AssociationProxy[list["Thing"]] = association_proxy(
+    things: AssociationProxy[list["Thing"]] = association_proxy(  # noqa: F821
         "thing_associations", "thing"
     )
 
@@ -95,7 +97,9 @@ class LocationThingAssociation(Base, AutoBaseMixin):
 
     # --- Relationship Definitions ---
     location: Mapped["Location"] = relationship(back_populates="thing_associations")
-    thing: Mapped["Thing"] = relationship(back_populates="location_associations")
+    thing: Mapped["Thing"] = relationship(  # noqa: F821
+        back_populates="location_associations"
+    )
 
 
 # ============= EOF =============================================

schemas/location.py

@@ -34,9 +34,10 @@ class CreateLocation(BaseCreateModel):
     Schema for creating a sample location.
     """
 
-    name: str | None = None
+    # name: str | None = None
     notes: str | None = None
     point: str  # point is required and should be in WKT format
+    elevation: float
     release_status: str | None = "draft"
     elevation_accuracy: float | None = None
     elevation_method: str | None = None
@@ -64,9 +65,12 @@ class LocationResponse(BaseResponseModel):
     Response schema for sample location details.
     """
 
-    name: str | None
+    # name: str | None
     notes: str | None
     point: str
+    elevation: float | None
+    horizontal_datum: str = "WGS84"
+    vertical_daum: str = "NAVD88"
     release_status: str | None
     elevation_accuracy: float | None
     elevation_method: str | None
@@ -103,9 +107,10 @@ class UpdateLocation(BaseUpdateModel):
     Schema for updating a location.
     """
 
-    name: str | None = None
+    # name: str | None = None
     notes: str | None = None
     point: str | None = None
+    elevation: float | None = None
     release_status: str | None = None
     elevation_accuracy: float | None = None
     elevation_method: str | None = None

services/geospatial_helper.py

@@ -13,11 +13,8 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 # ===============================================================================
-import json
-
 import shapefile
 from shapely.errors import GEOSException
-from geoalchemy2 import functions as geofunc
 from shapely.io import from_geojson
 
 import constants
@@ -46,7 +43,7 @@ def get_thing_features(
     #     selection_args.append(SpringThing)
 
     sql = (
-        select(Thing, ST_AsGeoJSON(Location.point).label("geojson"))
+        select(Thing, ST_AsGeoJSON(Location.point).label("geojson"), Location.elevation)
         .join(LocationThingAssociation, Thing.id == LocationThingAssociation.thing_id)
         .join(Location, LocationThingAssociation.location_id == Location.id)
     )
@@ -77,7 +74,7 @@ def create_shapefile(things: list, filename: str = "things.shp") -> None:
         shp.field("id", "L")
         shp.field("name", "C")
 
-        for thing, point in things:
+        for thing, point, elevation in things:
             # Assume loc.point is WKT or a Shapely geometry or GeoJSON
             if isinstance(point, str):
                 try:
@@ -88,7 +85,7 @@ def create_shapefile(things: list, filename: str = "things.shp") -> None:
                 geom = to_shape(point)
 
             shp.point(geom.x, geom.y)
-            shp.record(thing.id, thing.name)
+            shp.record(thing.id, thing.name, elevation)
 
 
 def make_within_wkt(sql: Select, wkt: str) -> Select:

services/location_helper.py

@@ -0,0 +1,27 @@
+from shapely.wkt import loads
+from pydantic import BaseModel
+from sqlalchemy.orm import Session
+
+from db.location import Location
+from services.util import (
+    get_state_from_point,
+    get_county_from_point,
+    get_quad_name_from_point,
+)
+
+
+def set_geographic_attributes(
+    session: Session, payload: BaseModel, location: Location
+) -> None:
+    """
+    Set geographic attributes for a location based off of the point. This function
+    is to be used for both POST and PATCH requests.
+    """
+    if payload.point is not None:
+        point = loads(payload.point)
+        longitude = point.x
+        latitude = point.y
+        location.state = get_state_from_point(longitude, latitude)
+        location.county = get_county_from_point(longitude, latitude)
+        location.quad_name = get_quad_name_from_point(longitude, latitude)
+        session.commit()

services/util.py

@@ -0,0 +1,110 @@
+from shapely.ops import transform
+import pyproj
+import httpx
+
+from constants import SRID_WGS84
+
+TRANSFORMERS = {}
+
+
+def transform_srid(geometry, source_srid, target_srid):
+    """
+    geometry must be a shapely geometry object, like Point, Polygon, or MultiPolygon
+    """
+    transformer_key = (source_srid, target_srid)
+    if transformer_key not in TRANSFORMERS:
+        source_crs = pyproj.CRS(f"EPSG:{source_srid}")
+        target_crs = pyproj.CRS(f"EPSG:{target_srid}")
+        transformer = pyproj.Transformer.from_crs(
+            source_crs, target_crs, always_xy=True
+        )
+        TRANSFORMERS[transformer_key] = transformer
+    else:
+        transformer = TRANSFORMERS[transformer_key]
+    return transform(transformer.transform, geometry)
+
+
+def get_tiger_data(
+    lon: float, lat: float, layer: int, outfields: str = "*"
+) -> dict | None:
+    url = f"https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/State_County/MapServer/{layer}/query"
+    params = {
+        "f": "json",
+        "where": "1=1",
+        "geometry": f"{lon},{lat}",
+        "geometryType": "esriGeometryPoint",
+        "inSR": f"{SRID_WGS84}",
+        "spatialRel": "esriSpatialRelIntersects",
+        "outFields": outfields,
+        "returnGeometry": "false",
+    }
+    resp = httpx.get(url, params=params, timeout=30)
+    data = resp.json()
+    if not data.get("features"):
+        return None
+
+    return data["features"][0]["attributes"]
+
+
+def get_state_from_point(lon: float, lat: float) -> str:
+    attrs = get_tiger_data(lon, lat, layer=0, outfields="BASENAME")
+    return attrs["BASENAME"]
+
+
+def get_county_from_point(lon: float, lat: float) -> str:
+    """
+    Look up county for a given longitude/latitude
+    using the US Census TIGERWeb REST API.
+    """
+
+    attrs = get_tiger_data(lon, lat, layer=1, outfields="BASENAME")
+    return attrs["BASENAME"]
+
+
+def get_quad_name_from_point(lon: float, lat: float) -> str:
+    url = "https://carto.nationalmap.gov/arcgis/rest/services/map_indices/MapServer/10/query"
+    params = {
+        "f": "json",
+        "geometry": f"{lon},{lat}",
+        "geometryType": "esriGeometryPoint",
+        "inSR": f"{SRID_WGS84}",
+        "spatialRel": "esriSpatialRelIntersects",
+        "outFields": "CELL_NAME,CELL_MAPCODE",
+        "returnGeometry": "false",
+    }
+
+    resp = httpx.get(url, params=params, timeout=30)
+    data = resp.json()
+
+    if data["features"]:
+        attrs = data["features"][0]["attributes"]
+        return attrs["CELL_NAME"]
+    else:
+        print(f"No quad name found for POINT ({lon} {lat})")
+        return None
+
+
+def get_epqs_elevation_from_point(lon: float, lat: float) -> float:
+    url = "https://epqs.nationalmap.gov/v1/json"
+    params = {
+        "x": lon,
+        "y": lat,
+        "units": "Meters",
+        "wkid": f"{SRID_WGS84}",
+        "includeDate": False,
+    }
+
+    resp = httpx.get(url, params=params)
+    data = resp.json()
+
+    return data["value"]
+
+
+if __name__ == "__main__":
+    x = -106.904107
+    y = 34.068198
+
+    print(get_state_from_point(x, y))
+    print(get_county_from_point(x, y))
+    print(get_quad_name_from_point(x, y))
+    print(get_epqs_elevation_from_point(x, y))