Geospatial Knowledge Graph: From Raw Geometry to Reasoning¶

# scripts/gen_data.py
import json, random
from pathlib import Path
from shapely.geometry import Point, LineString, Polygon, mapping

random.seed(7)
out = Path("scripts/out")
out.mkdir(parents=True, exist_ok=True)

# Cities (random around SF Bay area)
bbox = (-122.7, 37.2, -121.7, 37.9)
def rand_pt():
    x = random.uniform(bbox[0], bbox[2])
    y = random.uniform(bbox[1], bbox[3])
    return (x, y)

cities = []
for i, name in enumerate(["Alpha","Bravo","Charlie","Delta","Echo","Foxtrot","Golf","Hotel","India","Juliet"], start=1):
    x, y = rand_pt()
    cities.append({"id": f"C{i}", "name": name, "geom": mapping(Point(x, y))})

# Facilities near cities
facilities = []
for i in range(1, 21):
    c = random.choice(cities)
    px, py = c["geom"]["coordinates"]
    # Jitter ~0.01 degrees (~1km)
    fx, fy = px + random.uniform(-0.01, 0.01), py + random.uniform(-0.01, 0.01)
    facilities.append({
        "id": f"F{i}",
        "kind": random.choice(["warehouse", "hospital", "substation"]),
        "geom": mapping(Point(fx, fy))
    })

# Regions (polygons)
regions = []
for i, nm in enumerate(["NorthZone", "SouthZone", "EastPark", "WestReserve"], start=1):
    cx, cy = rand_pt()
    poly = Point(cx, cy).buffer(0.08, resolution=32)  # ~few km radius
    regions.append({
        "id": f"R{i}",
        "name": nm,
        "geom": mapping(Polygon(poly.exterior.coords))
    })

# Roads: connect random pairs of cities
roads = []
pairs = set()
for i in range(1, 20):
    a, b = random.sample(cities, 2)
    key = tuple(sorted([a["id"], b["id"]]))
    if key in pairs:
        continue
    pairs.add(key)
    pa = a["geom"]["coordinates"]
    pb = b["geom"]["coordinates"]
    # Slightly bend the line
    mid = (
        (pa[0] + pb[0]) / 2 + random.uniform(-0.02, 0.02),
        (pa[1] + pb[1]) / 2 + random.uniform(-0.02, 0.02)
    )
    geom = LineString([pa, mid, pb])
    roads.append({
        "id": f"E{i}",
        "a": a["id"],
        "b": b["id"],
        "geom": mapping(geom)
    })

# Export as GeoJSON
for name, data in [("cities", cities), ("facilities", facilities), ("regions", regions), ("roads", roads)]:
    Path(out / f"{name}.geojson").write_text(
        json.dumps({
            "type": "FeatureCollection",
            "features": [{
                "type": "Feature",
                "properties": {k: v for k, v in d.items() if k != 'geom'},
                "geometry": d["geom"]
            } for d in data]
        }, indent=2)
    )

print("Generated scripts/out/*.geojson")

# scripts/load_postgis.py
import json, os, time
import psycopg
from pathlib import Path

dsn = "postgresql://geo:geo@postgis:5432/geo"

schema_sql = """
CREATE EXTENSION IF NOT EXISTS postgis;

-- Drop existing tables
DROP TABLE IF EXISTS cities CASCADE;
DROP TABLE IF EXISTS facilities CASCADE;
DROP TABLE IF EXISTS regions CASCADE;
DROP TABLE IF EXISTS roads CASCADE;

-- Create tables with proper geometry types
CREATE TABLE cities(
    id text primary key, 
    name text, 
    geom geometry(Point, 4326)
);

CREATE TABLE facilities(
    id text primary key, 
    kind text, 
    geom geometry(Point, 4326)
);

CREATE TABLE regions(
    id text primary key, 
    name text, 
    geom geometry(Polygon, 4326)
);

CREATE TABLE roads(
    id text primary key, 
    a text, 
    b text, 
    geom geometry(LineString, 4326)
);

-- Create spatial indexes for performance
CREATE INDEX ON cities USING GIST (geom);
CREATE INDEX ON facilities USING GIST (geom);
CREATE INDEX ON regions USING GIST (geom);
CREATE INDEX ON roads USING GIST (geom);

-- Create additional indexes for common queries
CREATE INDEX ON cities (name);
CREATE INDEX ON facilities (kind);
CREATE INDEX ON regions (name);
"""

def load(fc_path, table, cols):
    """Load GeoJSON features into PostGIS"""
    fc = json.loads(Path(fc_path).read_text())
    rows = []
    for f in fc["features"]:
        props = f["properties"]
        geom = json.dumps(f["geometry"])
        rows.append(tuple([props[c] for c in cols] + [geom]))
    return rows

# Load data into PostGIS
with psycopg.connect(dsn) as conn:
    with conn.cursor() as cur:
        cur.execute(schema_sql)
        conn.commit()

    paths = Path("scripts/out")
    with conn.cursor() as cur:
        # Cities
        for rid, name, geom in load(paths / "cities.geojson", "cities", ["id", "name"]):
            cur.execute(
                "INSERT INTO cities(id, name, geom) VALUES (%s, %s, ST_SetSRID(ST_GeomFromGeoJSON(%s), 4326))",
                (rid, name, geom)
            )

        # Facilities
        for rid, kind, geom in load(paths / "facilities.geojson", "facilities", ["id", "kind"]):
            cur.execute(
                "INSERT INTO facilities(id, kind, geom) VALUES (%s, %s, ST_SetSRID(ST_GeomFromGeoJSON(%s), 4326))",
                (rid, kind, geom)
            )

        # Regions
        for rid, name, geom in load(paths / "regions.geojson", "regions", ["id", "name"]):
            cur.execute(
                "INSERT INTO regions(id, name, geom) VALUES (%s, %s, ST_SetSRID(ST_GeomFromGeoJSON(%s), 4326))",
                (rid, name, geom)
            )

        # Roads
        for rid, a, b, geom in load(paths / "roads.geojson", "roads", ["id", "a", "b"]):
            cur.execute(
                "INSERT INTO roads(id, a, b, geom) VALUES (%s, %s, %s, ST_SetSRID(ST_GeomFromGeoJSON(%s), 4326))",
                (rid, a, b, geom)
            )

    conn.commit()

print("Loaded GeoJSON into PostGIS with spatial indexes")

# scripts/build_kg.py
import math
import psycopg
from neo4j import GraphDatabase

PG = "postgresql://geo:geo@postgis:5432/geo"
BOLT = "bolt://neo4j:7687"
AUTH = ("neo4j", "test123")

def km(meters: float) -> float:
    """Convert meters to kilometers"""
    return round(meters / 1000.0, 3)

# Extract spatial relationships from PostGIS
with psycopg.connect(PG) as pg, GraphDatabase.driver(BOLT, auth=AUTH) as neo:
    cur = pg.cursor()

    # 1) Extract basic entities
    cur.execute("SELECT id, name, ST_X(geom), ST_Y(geom) FROM cities")
    cities = cur.fetchall()

    cur.execute("SELECT id, kind, ST_X(geom), ST_Y(geom) FROM facilities")
    facilities = cur.fetchall()

    cur.execute("SELECT id, name, ST_AsText(geom), ST_X(ST_Centroid(geom)), ST_Y(ST_Centroid(geom)) FROM regions")
    regions = cur.fetchall()

    cur.execute("""
        SELECT r.id, r.a, r.b, ST_Length(ST_Transform(r.geom, 3857)) AS meters
        FROM roads r
    """)
    roads = cur.fetchall()

    # 2) Compute spatial relationships
    cur.execute("""
        SELECT f.id, r.id
        FROM facilities f
        JOIN regions r ON ST_Within(f.geom, r.geom)
    """)
    fac_within = cur.fetchall()

    cur.execute("""
        SELECT c.id AS city, f.id AS fac, ST_DistanceSphere(c.geom, f.geom) AS meters
        FROM cities c, facilities f
        WHERE ST_DWithin(c.geom, f.geom, 2000.0)
    """)
    near_cf = cur.fetchall()

    # 3) Materialize graph in Neo4j
    with neo.session() as s:
        s.run("MATCH (n) DETACH DELETE n")  # Reset graph

        # Create constraints and indexes
        s.run("CREATE CONSTRAINT city_id IF NOT EXISTS FOR (c:City) REQUIRE c.id IS UNIQUE")
        s.run("CREATE CONSTRAINT fac_id IF NOT EXISTS FOR (f:Facility) REQUIRE f.id IS UNIQUE")
        s.run("CREATE CONSTRAINT region_id IF NOT EXISTS FOR (r:Region) REQUIRE r.id IS UNIQUE")
        s.run("CREATE INDEX city_loc IF NOT EXISTS FOR (c:City) ON (c.location)")

        # Create nodes
        for cid, name, x, y in cities:
            s.run("""
                MERGE (c:City {id: $id}) 
                SET c.name = $name, c.location = point({longitude: $x, latitude: $y})
            """, id=cid, name=name, x=float(x), y=float(y))

        for fid, kind, x, y in facilities:
            s.run("""
                MERGE (f:Facility {id: $id}) 
                SET f.kind = $kind, f.location = point({longitude: $x, latitude: $y})
            """, id=fid, kind=kind, x=float(x), y=float(y))

        for rid, name, wkt, cx, cy in regions:
            s.run("""
                MERGE (r:Region {id: $id})
                SET r.name = $name,
                    r.centroid = point({longitude: $cx, latitude: $cy}),
                    r.wkt = $wkt
            """, id=rid, name=name, wkt=wkt, cx=float(cx), cy=float(cy))

        # Create relationships
        for fid, rid in fac_within:
            s.run("""
                MATCH (f:Facility {id: $fid}), (r:Region {id: $rid})
                MERGE (f)-[:WITHIN]->(r)
            """, fid=fid, rid=rid)

        for cid, fid, meters in near_cf:
            s.run("""
                MATCH (c:City {id: $cid}), (f:Facility {id: $fid})
                MERGE (c)-[n:NEAR]->(f)
                SET n.km = $km
            """, cid=cid, fid=fid, km=km(float(meters)))

        # Create road connections (bidirectional)
        for rid, a, b, meters in roads:
            s.run("""
                MATCH (c1:City {id: $a}), (c2:City {id: $b})
                MERGE (c1)-[e:CONNECTED_TO {road_id: $rid}]->(c2)
                SET e.km = $km
            """, a=a, b=b, rid=rid, km=km(float(meters)))

            s.run("""
                MATCH (c1:City {id: $a}), (c2:City {id: $b})
                MERGE (c2)-[e:CONNECTED_TO {road_id: $rid}]->(c1)
                SET e.km = $km
            """, a=a, b=b, rid=rid, km=km(float(meters)))

print("Knowledge graph built in Neo4j")

# Optimized ETL for large datasets
class OptimizedGeospatialETL:
    def __init__(self, postgis_dsn, neo4j_uri, neo4j_auth):
        self.pg_dsn = postgis_dsn
        self.neo4j_uri = neo4j_uri
        self.neo4j_auth = neo4j_auth
        self.batch_size = 1000

    def batch_load_relationships(self, relationships, rel_type):
        """Batch load relationships for performance"""
        with GraphDatabase.driver(self.neo4j_uri, auth=self.neo4j_auth) as driver:
            with driver.session() as session:
                # Use UNWIND for batch operations
                query = f"""
                UNWIND $relationships AS rel
                MATCH (a {{id: rel.from}}), (b {{id: rel.to}})
                MERGE (a)-[r:{rel_type}]->(b)
                SET r += rel.properties
                """
                session.run(query, relationships=relationships)

    def spatial_index_optimization(self):
        """Optimize spatial queries with proper indexing"""
        with psycopg.connect(self.pg_dsn) as conn:
            with conn.cursor() as cur:
                # Create composite indexes for common queries
                cur.execute("""
                    CREATE INDEX IF NOT EXISTS idx_cities_location 
                    ON cities USING GIST (geom);
                """)
                cur.execute("""
                    CREATE INDEX IF NOT EXISTS idx_facilities_kind_location 
                    ON facilities (kind) INCLUDE (geom);
                """)
                cur.execute("""
                    CREATE INDEX IF NOT EXISTS idx_roads_ab 
                    ON roads (a, b);
                """)

    def incremental_updates(self, timestamp):
        """Handle incremental updates to the knowledge graph"""
        with psycopg.connect(self.pg_dsn) as conn:
            with conn.cursor() as cur:
                # Find entities modified since timestamp
                cur.execute("""
                    SELECT id, 'city' as type, name, ST_X(geom), ST_Y(geom)
                    FROM cities 
                    WHERE updated_at > %s
                    UNION ALL
                    SELECT id, 'facility' as type, kind, ST_X(geom), ST_Y(geom)
                    FROM facilities 
                    WHERE updated_at > %s
                """, (timestamp, timestamp))

                modified_entities = cur.fetchall()

                # Update Neo4j with only changed entities
                self.update_neo4j_entities(modified_entities)

# Monitoring for geospatial knowledge graphs
class GeospatialKGMonitor:
    def __init__(self, postgis_dsn, neo4j_uri, neo4j_auth):
        self.pg_dsn = postgis_dsn
        self.neo4j_uri = neo4j_uri
        self.neo4j_auth = neo4j_auth

    def check_data_quality(self):
        """Validate data quality across both systems"""
        issues = []

        # Check PostGIS data quality
        with psycopg.connect(self.pg_dsn) as conn:
            with conn.cursor() as cur:
                cur.execute("SELECT COUNT(*) FROM cities WHERE geom IS NULL")
                null_geoms = cur.fetchone()[0]
                if null_geoms > 0:
                    issues.append(f"Found {null_geoms} cities with null geometry")

                cur.execute("SELECT COUNT(*) FROM facilities WHERE NOT ST_IsValid(geom)")
                invalid_geoms = cur.fetchone()[0]
                if invalid_geoms > 0:
                    issues.append(f"Found {invalid_geoms} facilities with invalid geometry")

        # Check Neo4j data quality
        with GraphDatabase.driver(self.neo4j_uri, auth=self.neo4j_auth) as driver:
            with driver.session() as session:
                result = session.run("MATCH (n) RETURN count(n) as total_nodes")
                total_nodes = result.single()["total_nodes"]

                result = session.run("MATCH ()-[r]->() RETURN count(r) as total_rels")
                total_rels = result.single()["total_rels"]

                if total_nodes == 0:
                    issues.append("Neo4j has no nodes")
                if total_rels == 0:
                    issues.append("Neo4j has no relationships")

        return issues

    def performance_metrics(self):
        """Collect performance metrics"""
        metrics = {}

        # PostGIS performance
        with psycopg.connect(self.pg_dsn) as conn:
            with conn.cursor() as cur:
                cur.execute("SELECT COUNT(*) FROM cities")
                metrics['cities_count'] = cur.fetchone()[0]

                cur.execute("SELECT COUNT(*) FROM facilities")
                metrics['facilities_count'] = cur.fetchone()[0]

                cur.execute("SELECT COUNT(*) FROM roads")
                metrics['roads_count'] = cur.fetchone()[0]

        # Neo4j performance
        with GraphDatabase.driver(self.neo4j_uri, auth=self.neo4j_auth) as driver:
            with driver.session() as session:
                result = session.run("MATCH (n) RETURN count(n) as total_nodes")
                metrics['neo4j_nodes'] = result.single()["total_nodes"]

                result = session.run("MATCH ()-[r]->() RETURN count(r) as total_rels")
                metrics['neo4j_relationships'] = result.single()["total_rels"]

        return metrics

# ❌ WRONG: No geometry validation
def bad_geometry_handling():
    # Loading without validation
    cur.execute("INSERT INTO cities (id, name, geom) VALUES (%s, %s, %s)", 
                (id, name, geom_json))
    # No validation of geometry validity!

# ✅ CORRECT: Proper geometry validation
def good_geometry_handling():
    # Validate geometry before insertion
    cur.execute("""
        INSERT INTO cities (id, name, geom) 
        VALUES (%s, %s, ST_SetSRID(ST_GeomFromGeoJSON(%s), 4326))
        WHERE ST_IsValid(ST_GeomFromGeoJSON(%s))
    """, (id, name, geom_json, geom_json))

# ❌ WRONG: Inefficient spatial queries
def bad_spatial_query():
    # No spatial index usage
    cur.execute("""
        SELECT * FROM cities c, facilities f 
        WHERE ST_Distance(c.geom, f.geom) < 1000
    """)
    # This will be very slow!

# ✅ CORRECT: Efficient spatial queries
def good_spatial_query():
    # Use spatial indexes with ST_DWithin
    cur.execute("""
        SELECT c.id, f.id, ST_Distance(c.geom, f.geom) as distance
        FROM cities c, facilities f 
        WHERE ST_DWithin(c.geom, f.geom, 1000)
    """)

# ✅ DO: Complete spatial data management checklist
def spatial_data_checklist():
    checklist = {
        "geometry_quality": [
            "Validate geometry validity (ST_IsValid)",
            "Check coordinate system consistency",
            "Ensure proper spatial indexing",
            "Monitor geometry complexity"
        ],
        "spatial_operations": [
            "Use ST_DWithin for distance queries",
            "Leverage spatial indexes (GIST)",
            "Optimize coordinate transformations",
            "Batch spatial operations"
        ],
        "graph_consistency": [
            "Validate relationship properties",
            "Check for orphaned nodes",
            "Ensure bidirectional relationships",
            "Monitor graph connectivity"
        ],
        "performance": [
            "Profile spatial queries",
            "Optimize batch operations",
            "Monitor memory usage",
            "Scale incrementally"
        ]
    }
    return checklist

# ✅ DO: Production deployment checklist
class ProductionChecklist:
    @staticmethod
    def infrastructure_checklist():
        return {
            "postgis": [
                "Proper memory configuration",
                "Spatial index optimization",
                "Connection pooling",
                "Backup and recovery"
            ],
            "neo4j": [
                "Memory heap configuration",
                "Graph Data Science plugins",
                "Index optimization",
                "Clustering setup"
            ],
            "monitoring": [
                "Performance metrics",
                "Data quality checks",
                "Error alerting",
                "Capacity planning"
            ]
        }

    @staticmethod
    def data_quality_checks():
        return {
            "geometry_validation": "ST_IsValid(geom)",
            "coordinate_consistency": "ST_SRID(geom) = 4326",
            "relationship_integrity": "MATCH (n)-[r]->(m) WHERE n IS NULL OR m IS NULL",
            "spatial_index_usage": "EXPLAIN ANALYZE for spatial queries"
        }

# Spatial relationship computation
def compute_spatial_relationships():
    # PostGIS: Compute relationships
    cur.execute("""
        SELECT c.id, f.id, ST_DistanceSphere(c.geom, f.geom) as meters
        FROM cities c, facilities f
        WHERE ST_DWithin(c.geom, f.geom, 2000)
    """)

    # Neo4j: Materialize relationships
    for city_id, facility_id, meters in cur.fetchall():
        session.run("""
            MATCH (c:City {id: $city_id}), (f:Facility {id: $facility_id})
            MERGE (c)-[r:NEAR]->(f)
            SET r.km = $km
        """, city_id=city_id, facility_id=facility_id, km=meters/1000)