Python实现st_project与PostGIS st_project结果不一致问题排查

Let's break down why your Python implementation isn't matching PostGIS's output, and fix the issues step by step:

Key Differences to Investigate #

1. Coordinate Order Mismatch (Most Likely Culprit) #

geopy's Point uses the order (latitude, longitude), while PostGIS's ST_MakePoint requires (longitude, latitude). If you flipped these in your PostGIS query, you're calculating a point in a completely different location, leading to unrelated results.

2. Distance Parameter/Calculation Issues #

Your geopy output shows a longitude jump of ~0.09 degrees, which translates to ~9 kilometers—way farther than the 50 meters you intended. This suggests:

• You might have a typo in your distance parameter (e.g., accidentally using kilometers=50 instead of meters=50).
• You're using the deprecated VincentyDistance (geopy has marked this as obsolete in favor of GeodesicDistance, which uses more accurate WGS84-based calculations).

3. Azimuth Unit or Definition #

Both tools define azimuth as degrees clockwise from true north, but double-check:

• In geopy, you correctly converted 90° (east) to radians with radians(90).
• In PostGIS, ensure you're also passing radians (not raw degrees) to ST_Project's third parameter.

4. Ellipsoid Model Consistency #

Both geopy and PostGIS default to WGS84, but if you've overridden the ellipsoid in either tool, calculations will diverge.

Fixed Implementation #

Corrected geopy Code #

Use the recommended GeodesicDistance for accurate, up-to-date calculations:

import geopy
from geopy.distance import geodesic
from geopy.units import radians

# geopy Point: (latitude, longitude)
start = geopy.Point(27.725778916300465, 35.201911926269524)
# 50 meters east (90° azimuth, converted to radians)
result = geodesic(meters=50).destination(start, radians(90))
print(result)

This will output a point with a tiny longitude increase (~0.0005 degrees) matching the 50-meter eastward movement.

Matching PostGIS Query #

Ensure you use the correct coordinate order (lon, lat) and parameters:

SELECT ST_AsText(
    ST_Project(
        ST_SetSRID(ST_MakePoint(35.201911926269524, 27.725778916300465), 4326),
        50,  -- Distance in meters
        radians(90)  -- Eastward azimuth in radians
    )
);

Why Your Original geopy Code Failed #

Your original output suggests the distance calculation wasn't using 50 meters—likely due to a parameter bug in the deprecated VincentyDistance or a hidden typo. Switching to geodesic eliminates this issue and aligns with PostGIS's underlying calculations.

内容的提问来源于stack exchange，提问作者Chems Bezzaz