如何使用PostgreSQL或Python匹配经纬度日志与对应六边形网格ID

Got it, let's solve this problem of mapping each GPS log entry to its corresponding hexagonal grid ID. I'll walk through both PostgreSQL (using PostGIS for geospatial processing) and Python implementations, since those are the two tools you asked for.

First, you'll need the PostGIS extension installed in your PostgreSQL database—it's essential for handling spatial data like points and polygons. Here's a step-by-step approach:

1. Convert raw data to spatial types

   • For the log table (data1), we’ll turn the lat/lon columns into a PostGIS POINT geometry (using the standard WGS84 coordinate system, SRID 4326).
   • For the hexagonal grid table (data2), we’ll take the 6 corner points and convert them into a valid POLYGON geometry. Make sure the points are ordered clockwise or counter-clockwise—this is critical for the polygon to be recognized correctly.

2. Spatial join to match points to grids
   Use ST_Contains to check which grid polygon contains the GPS point, or ST_Within to verify if the point lies inside the polygon—both work, pick whichever feels more intuitive to you.

Here's the full SQL code example:

-- Enable PostGIS if it's not already active
CREATE EXTENSION IF NOT EXISTS postgis;

-- Create a temporary view for log data with spatial point geometry
WITH log_spatial AS (
    SELECT 
        log_id,
        ST_SetSRID(ST_MakePoint(lon, lat), 4326) AS geom
    FROM data1
),
-- Create a temporary view for grid data with hexagonal polygon geometry
grid_spatial AS (
    SELECT 
        grid_id,
        ST_SetSRID(
            ST_MakePolygon(
                ST_MakeLine(ARRAY[
                    ST_MakePoint(point1_lon, point1_lat),
                    ST_MakePoint(point2_lon, point2_lat),
                    ST_MakePoint(point3_lon, point3_lat),
                    ST_MakePoint(point4_lon, point4_lat),
                    ST_MakePoint(point5_lon, point5_lat),
                    ST_MakePoint(point6_lon, point6_lat),
                    ST_MakePoint(point1_lon, point1_lat) -- Close the polygon by repeating the first point
                ])
            ), 4326
        ) AS geom
    FROM data2
)
-- Join the two views to get the matching grid IDs for each log
SELECT 
    ls.log_id,
    gs.grid_id
FROM log_spatial ls
JOIN grid_spatial gs ON ST_Contains(gs.geom, ls.geom)
ORDER BY ls.log_id;

Note: Double-check the order of your grid's corner points. If the polygon is invalid (self-intersecting), ST_Contains won’t work. You can test validity with ST_IsValid(geom).

For Python, we’ll use geopandas (which wraps shapely for spatial operations) since it makes handling spatial dataframes straightforward. First, install the required packages if you haven’t already:

pip install geopandas shapely pandas

Here's the step-by-step code:

import geopandas as gpd
from shapely.geometry import Point, Polygon
import pandas as pd

# 1. Load your datasets (replace paths with your actual data sources)
df_logs = pd.read_csv("data1.csv")
df_grids = pd.read_csv("data2.csv")

# 2. Convert log data to spatial points (Shapely uses (lon, lat) for Point coordinates)
def create_point(row):
    return Point(row["lon"], row["lat"])

df_logs["geometry"] = df_logs.apply(create_point, axis=1)
gdf_logs = gpd.GeoDataFrame(df_logs, crs="EPSG:4326")

# 3. Convert grid data to hexagonal polygons
def create_hexagon(row):
    # Collect 6 corner points in order (ensure no crossing lines!)
    points = [
        (row["point1_lon"], row["point1_lat"]),
        (row["point2_lon"], row["point2_lat"]),
        (row["point3_lon"], row["point3_lat"]),
        (row["point4_lon"], row["point4_lat"]),
        (row["point5_lon"], row["point5_lat"]),
        (row["point6_lon"], row["point6_lat"]),
        (row["point1_lon"], row["point1_lat"])  # Close the polygon
    ]
    return Polygon(points)

df_grids["geometry"] = df_grids.apply(create_hexagon, axis=1)
gdf_grids = gpd.GeoDataFrame(df_grids, crs="EPSG:4326")

# 4. Perform spatial join to match logs to grids
matched = gpd.sjoin(gdf_logs, gdf_grids, how="left", predicate="within")

# 5. Output the required columns
result = matched[["log_id", "grid_id"]]
print(result)
# Save to CSV if needed: result.to_csv("matched_logs_grids.csv", index=False)

Note: If you get no matches or incorrect ones, verify the order of your grid's corner points. You can plot the polygons with gdf_grids.plot() to visually check for issues.

内容的提问来源于stack exchange，提问作者H.S.