DEV Community

EME GUG
EME GUG

Posted on

Stop Parsing Address Strings: Fast Point-in-Polygon Lookups with PostGIS and Shapely

#postgis #geospatial #database #python

Tuần này HN rộ lên tin USPS chính thức công bố polygon ranh giới cho ZIP Code. Với dev Mỹ, đây là chuyện vui. Với dev Việt Nam, nó nhắc tới một nỗi đau quen thuộc: sau đợt sắp xếp lại đơn vị hành chính năm 2025 (còn 34 tỉnh thành, bỏ cấp huyện, phường xã gộp lại hàng loạt), mọi hệ thống lưu địa chỉ dạng text như Phường 7, Quận 3, TP.HCM đều bị lỗi thời chỉ sau một đêm. Team mình từng mất gần hai sprint chỉ để viết regex và bảng mapping cho chuyện này. Bài học rút ra: đừng suy ra khu vực hành chính từ chuỗi địa chỉ. Hãy lưu tọa độ, rồi tra bằng polygon. Bài này chia sẻ cách mình làm point-in-polygon lookup với PostGIS và Shapely, chạy production với vài triệu request mỗi ngày.

Tại sao lưu tọa độ + polygon tốt hơn lưu text

Địa chỉ dạng text có ba vấn đề lớn:

  • Không chuẩn hóa được: Q.3, Quận 3, Q3, District 3 đều là một chỗ. Fuzzy match thì sai với các tên trùng như Phường 1 (có ở hàng chục quận).
  • Ranh giới thay đổi: khi nhà nước gộp phường, text cũ vẫn nằm trong DB nhưng không còn đúng.
  • Không trả lời được câu hỏi không gian: shipper này có nằm trong vùng giao hàng không? Đơn hàng này thuộc kho nào?

Cách làm đúng: lưu lat/lng (lấy từ GPS của app hoặc geocode một lần), và giữ một bảng polygon ranh giới có version. Khi ranh giới đổi, bạn chỉ cần import bộ polygon mới và tra lại, dữ liệu gốc (tọa độ) không bao giờ sai.

flowchart LR
    A[App gửi lat/lng] --> B[API lookup]
    B --> C{Cache hit?}
    C -- Có --> D[Trả kết quả]
    C -- Không --> E[PostGIS ST_Contains]
    E --> F[(Bảng boundaries có GiST index)]
    F --> E
    E --> G[Ghi cache]
    G --> D

Import polygon vào PostGIS bằng ogr2ogr

Dữ liệu ranh giới thường có dạng GeoJSON hoặc Shapefile (nguồn mở như OpenStreetMap qua Overpass, hoặc dữ liệu từ cơ quan nhà nước). Mình dùng PostgreSQL 17 + PostGIS 3.5, chạy bằng image postgis/postgis:17-3.5, và GDAL 3.10 để import.

# Khởi động PostGIS local
docker run -d --name gis -p 5432:5432 \
  -e POSTGRES_PASSWORD=dev postgis/postgis:17-3.5

# Import GeoJSON vào bảng boundaries, ép về SRID 4326
ogr2ogr -f PostgreSQL \
  PG:'host=localhost user=postgres password=dev dbname=postgres' \
  wards_2025.geojson \
  -nln boundaries -nlt MULTIPOLYGON \
  -t_srs EPSG:4326 -lco GEOMETRY_NAME=geom \
  -lco SPATIAL_INDEX=GIST -overwrite

# Kiểm tra geometry hỏng (rất hay gặp với dữ liệu tự vẽ)
psql -h localhost -U postgres -c \
  'SELECT count(*) FROM boundaries WHERE NOT ST_IsValid(geom);'
Enter fullscreen mode Exit fullscreen mode

Nếu có polygon invalid (self-intersection, ring không đóng), sửa ngay bằng UPDATE boundaries SET geom = ST_MakeValid(geom) WHERE NOT ST_IsValid(geom);. Polygon hỏng là nguyên nhân số một khiến ST_Contains trả kết quả sai mà không báo lỗi.

Một mẹo quan trọng: thêm cột version (ví dụ 2025-07) và valid_from. Khi có bộ ranh giới mới, import vào cùng bảng với version mới thay vì overwrite. Như vậy bạn vẫn tra được đơn hàng năm 2024 thuộc phường nào theo ranh giới cũ, rất cần cho báo cáo thuế và đối soát.

Query: ST_Contains, index và cái bẫy hiệu năng

Query cơ bản rất ngắn:

SELECT ward_code, ward_name, province_name
FROM boundaries
WHERE version = '2025-07'
  AND ST_Contains(geom, ST_SetSRID(ST_MakePoint(106.6893, 10.7769), 4326))
LIMIT 1;
Enter fullscreen mode Exit fullscreen mode

Lưu ý thứ tự: ST_MakePoint(lng, lat), kinh độ trước, vĩ độ sau. Mình đã thấy ít nhất ba team đảo ngược cái này, và kết quả là điểm rơi xuống giữa Ấn Độ Dương.

Về hiệu năng, GiST index giúp PostGIS lọc theo bounding box trước rồi mới kiểm tra chính xác. Nhưng polygon ranh giới hành chính thường có hàng chục nghìn vertex (đường bờ biển, sông), nên bước kiểm tra chính xác vẫn chậm. Giải pháp là ST_Subdivide:

CREATE TABLE boundaries_sub AS
SELECT ward_code, version, ST_Subdivide(geom, 256) AS geom
FROM boundaries;

CREATE INDEX ON boundaries_sub USING GIST (geom);
ANALYZE boundaries_sub;
Enter fullscreen mode Exit fullscreen mode

Với bộ khoảng 3.300 phường xã, việc chia nhỏ polygon thành các mảnh tối đa 256 vertex giúp p99 latency của mình giảm từ khoảng 40ms xuống dưới 2ms. Bạn tra trên bảng boundaries_sub, lấy ward_code, rồi join về bảng gốc để lấy tên.

Khi không muốn gọi DB: in-memory lookup với Shapely 2.0

Với batch job (gán phường cho 10 triệu đơn hàng cũ) hoặc service cần latency cực thấp, mình load polygon vào RAM và dùng STRtree của Shapely 2.0. API mới của Shapely 2.x hỗ trợ vectorized query, nhanh hơn vòng lặp Python rất nhiều.

# pip install shapely==2.0.6 numpy
import json
import numpy as np
from shapely import STRtree, points
from shapely.geometry import shape

with open('wards_2025.geojson', encoding='utf-8') as f:
    features = json.load(f)['features']

geoms = [shape(ft['geometry']) for ft in features]
codes = np.array([ft['properties']['ward_code'] for ft in features])
tree = STRtree(geoms)

def lookup_batch(lngs, lats):
    pts = points(np.column_stack([lngs, lats]))
    # predicate='within' = điểm nằm trong polygon
    pt_idx, poly_idx = tree.query(pts, predicate='within')
    result = np.full(len(pts), None, dtype=object)
    result[pt_idx] = codes[poly_idx]
    return result

print(lookup_batch([106.6893, 105.8542], [10.7769, 21.0285]))
Enter fullscreen mode Exit fullscreen mode

Trên MacBook M2, đoạn này xử lý khoảng 1 triệu điểm trong vài giây. Kết quả None nghĩa là điểm nằm ngoài mọi polygon: thường là tọa độ GPS lỗi, hoặc điểm nằm trên biển sát bờ. Đừng bỏ qua chúng, hãy log lại để kiểm tra.

flowchart TD
    A[Cần tra point-in-polygon] --> B{Dữ liệu đổi thường xuyên?}
    B -- Có --> C[PostGIS + version column]
    B -- Không --> D{Batch hay realtime?}
    D -- Batch --> E[Shapely STRtree vectorized]
    D -- Realtime latency thấp --> F[Shapely in-memory + reload định kỳ]
    C --> G[ST_Subdivide + GiST index]

Những cái bẫy mình đã dính

  • Điểm nằm đúng trên ranh giới: ST_Contains trả false cho điểm nằm trên cạnh. Nếu cần chắc chắn có kết quả, dùng ST_Intersects kèm ORDER BY ward_code LIMIT 1 để kết quả ổn định.
  • Polygon chồng lấn: dữ liệu từ nhiều nguồn có thể chồng lên nhau vài mét. Chạy ST_Overlaps giữa các polygon một lần sau mỗi lần import để phát hiện.
  • Đơn vị khoảng cách: SRID 4326 tính bằng độ, không phải mét. Muốn tìm điểm trong bán kính 500m, cast sang geography: ST_DWithin(geom::geography, pt::geography, 500).
  • Cache theo grid: làm tròn tọa độ đến 4 chữ số thập phân (khoảng 11m) để làm cache key trong Redis. Hit rate của mình đạt khoảng 70% vì người dùng tập trung ở khu đông dân.

Kết luận

Tin ZIP Code polygon của USPS là lời nhắc tốt: dữ liệu địa lý nên được biểu diễn bằng hình học, không phải bằng chuỗi ký tự. Nếu hệ thống của bạn đang vật lộn với đợt đổi địa giới hành chính, đây là các bước nên làm ngay:

  1. Bắt đầu lưu lat/lng cho mọi địa chỉ mới, kể cả khi chưa dùng đến.
  2. Import ranh giới vào PostGIS bằng ogr2ogr, luôn chạy ST_IsValid và ST_MakeValid sau khi import.
  3. Thêm cột version cho bảng ranh giới, không bao giờ overwrite dữ liệu cũ.
  4. Dùng ST_Subdivide + GiST index khi polygon phức tạp, latency sẽ giảm hàng chục lần.
  5. Batch job thì dùng Shapely 2.0 STRtree với vectorized query thay vì gọi DB từng điểm.
  6. Kiểm tra thứ tự lng, lat trong mọi chỗ tạo point, viết unit test với một tọa độ đã biết như Bến Thành hoặc Hồ Gươm.

Lần tới khi ranh giới thay đổi, bạn chỉ cần import file mới và chạy lại một câu query, thay vì mất hai sprint sửa regex như team mình ngày trước.

Top comments (0)