CODEWITHSUNDEEP
postgresqlpostgis
Ashish Gamit
Ashish Gamit

Reputation: 55

Distance calculated wrongly with the ST_Distance function in PostgreSQL database?

We are developing an application that involves calculating the shortest distance for a given point against the polygons stored within the PostgreSQL table.

We are using the ST_Distance function from PostgreSQL database.

When we compare the calculated distance with the Google Earth, there is vast difference between them.

Of course, neither Google nor PostgreSQL could be wrong (or are they?), so we are obviously missing something here.

Any ideas what is going wrong?

I have given below the sample query that we have tested with PostgreSQL along with the screenshots from Google Earth.

SELECT ST_Distance(Place1, Place2) As Place1ToPlace2
    , ST_Distance(Place1, Spot1) As Place1ToSpot1
    , ST_Distance(Place1, Spot2) As Place1ToSpot2
    , ST_Distance(Place2, Spot1) As Place2ToSpot1
    , ST_Distance(Place2, Spot2) As Place2ToSpot2
FROM (SELECT
    ST_PolygonFromText('SRID=4326;POLYGON((-74.0050636293915 40.75265123968514,-74.00500355126653 40.75268991743845,-74.00498169169283 40.75267084386348,-74.00503571044075 40.75263867886528,-74.0050636293915 40.75265123968514))') as Spot1
    ,ST_PolygonFromText('SRID=4326;POLYGON((-74.00503571044075 40.75263867886528,-74.00498225451273 40.75267084385684,-74.00495878551709 40.75265859837483,-74.00501023946696 40.75262521978885,-74.00503571044075 40.75263867886528))') as Spot2
    ,ST_GeogFromText('SRID=4326;POINT(-74.00489 40.752894)') As Place1
    ,ST_GeogFromText('SRID=4326;POINT(-74.004774 40.752846)') As Place2
    ) As foo  ;

It results in following values:

place1toplace2 |place1tospot1 |place1tospot2 |place2tospot1 |place2tospot2 |
---------------|--------------|--------------|--------------|--------------|
11.152362504   |24.608417285  |25.977083731  |26.004190091  |26.011579435  |

Following are the screenshots from Google Earth:

  1. Place1ToPlace2
  2. Place1ToSpot1
  3. Place1ToSpot2
  4. Place2ToSpot1
  5. Place2ToSpot2

Thanks in advance!!

Following are the KMLs exported from Google Earth:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">
<Document>
    <name>Spot 1.kml</name>
    <Style id="inline">
        <LineStyle>
            <color>ff0000ff</color>
            <width>2</width>
        </LineStyle>
        <PolyStyle>
            <fill>0</fill>
        </PolyStyle>
    </Style>
    <Style id="inline0">
        <LineStyle>
            <color>ff0000ff</color>
            <width>2</width>
        </LineStyle>
        <PolyStyle>
            <fill>0</fill>
        </PolyStyle>
    </Style>
    <StyleMap id="inline1">
        <Pair>
            <key>normal</key>
            <styleUrl>#inline</styleUrl>
        </Pair>
        <Pair>
            <key>highlight</key>
            <styleUrl>#inline0</styleUrl>
        </Pair>
    </StyleMap>
    <Placemark>
        <name>Spot 1</name>
        <styleUrl>#inline1</styleUrl>
        <Polygon>
            <tessellate>1</tessellate>
            <outerBoundaryIs>
                <LinearRing>
                    <coordinates>
                        -74.0050636293915,40.75265123968514,0 -74.00500355126653,40.75268991743845,0 -74.00498169169283,40.75267084386348,0 -74.00503571044075,40.75263867886528,0 -74.0050636293915,40.75265123968514,0 
                    </coordinates>
                </LinearRing>
            </outerBoundaryIs>
        </Polygon>
    </Placemark>
</Document>
</kml>


<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">
<Document>
    <name>Spot 2.kml</name>
    <Style id="inline">
        <LineStyle>
            <color>ff0000ff</color>
            <width>2</width>
        </LineStyle>
        <PolyStyle>
            <fill>0</fill>
        </PolyStyle>
    </Style>
    <StyleMap id="inline0">
        <Pair>
            <key>normal</key>
            <styleUrl>#inline1</styleUrl>
        </Pair>
        <Pair>
            <key>highlight</key>
            <styleUrl>#inline</styleUrl>
        </Pair>
    </StyleMap>
    <Style id="inline1">
        <LineStyle>
            <color>ff0000ff</color>
            <width>2</width>
        </LineStyle>
        <PolyStyle>
            <fill>0</fill>
        </PolyStyle>
    </Style>
    <Placemark>
        <name>Spot 2</name>
        <styleUrl>#inline0</styleUrl>
        <Polygon>
            <tessellate>1</tessellate>
            <outerBoundaryIs>
                <LinearRing>
                    <coordinates>
                        -74.00503571044075,40.75263867886528,0 -74.00498225451273,40.75267084385684,0 -74.00495878551709,40.75265859837483,0 -74.00501023946696,40.75262521978885,0 -74.00503571044075,40.75263867886528,0 
                    </coordinates>
                </LinearRing>
            </outerBoundaryIs>
        </Polygon>
    </Placemark>
</Document>
</kml>

!EDIT 2!

I have tried one more thing. I have taken the point from Spot1 that visually looks nearest to the Place1 and marked it as Place3.

When we look it on the map, it looks like Place3 is nearer to Place1 than Place2, but when checking with the query, distance to Place3 gives higher value.

I have checked with following query:

SELECT ST_Distance(Place1, Place2) As Place1ToPlace2
,ST_Distance(Place1, Place3) As Place1ToPlace3
FROM (SELECT
    ST_GeogFromText('SRID=4326;POINT(-74.00489 40.752894)') As Place1
    ,ST_GeogFromText('SRID=4326;POINT(-74.004774 40.752846)') As Place2
    ,ST_GeogFromText('SRID=4326;POINT(-74.00500355126653 40.75268991743845)') As Place3
    ) As foo;

And it gives following result:

place1toplace2 |place1toplace3 |
---------------|---------------|
11.152362504   |24.608417285   |

While on the Map: Place1ToPlace3

Following is the KML for Place1,Place2 and Place3 (I have removed the style related tags from KML)

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">
<Document>
    <name>Distance Comparison.kml</name>
    <Folder>
        <name>Distance Comparison</name>
        <open>1</open>
        <Style>
            <ListStyle>
                <listItemType>check</listItemType>
                <bgColor>00ffffff</bgColor>
                <maxSnippetLines>2</maxSnippetLines>
            </ListStyle>
        </Style>
        <Placemark>
            <name>Place 1 -  40°45&apos;9.78&quot;N  74° 0&apos;18.07&quot;W (40.752894, -74.00489)</name>
            <open>1</open>
            <LookAt>
                <longitude>-74.00500758183839</longitude>
                <latitude>40.75269419172616</latitude>
                <altitude>0</altitude>
                <heading>-0.0008536233435993688</heading>
                <tilt>29.8433509629012</tilt>
                <range>47.16429940085073</range>
                <gx:altitudeMode>relativeToSeaFloor</gx:altitudeMode>
            </LookAt>
            <styleUrl>#msn_1</styleUrl>
            <Point>
                <gx:drawOrder>1</gx:drawOrder>
                <coordinates>-74.00501944444444,40.75271666666666,0</coordinates>
            </Point>
        </Placemark>
        <Placemark>
            <name>Place 2 - 40°45&apos;9.54&quot;N 74° 0&apos;17.70&quot;W (40.752846, -74.004774)</name>
            <open>1</open>
            <LookAt>
                <longitude>-74.00500758183839</longitude>
                <latitude>40.75269419172616</latitude>
                <altitude>0</altitude>
                <heading>-0.0008536233435993688</heading>
                <tilt>29.8433509629012</tilt>
                <range>47.16429940085073</range>
                <gx:altitudeMode>relativeToSeaFloor</gx:altitudeMode>
            </LookAt>
            <styleUrl>#msn_2</styleUrl>
            <Point>
                <gx:drawOrder>1</gx:drawOrder>
                <coordinates>-74.00491666666667,40.75265,0</coordinates>
            </Point>
        </Placemark>
        <Placemark>
            <name>Place 3 -  40°45&apos;9.68&quot;N  74° 0&apos;18.01&quot;W (40.75268991743845 -74.00500355126653)</name>
            <open>1</open>
            <styleUrl>#msn_3</styleUrl>
            <Point>
                <coordinates>-74.00500277777778,40.75268888888889,0</coordinates>
            </Point>
        </Placemark>
    </Folder>
</Document>
</kml>

Upvotes: 2

Views: 1424

Answers (2)

mlinth
mlinth

Reputation: 3118

Your features in PostGIS don't appear to be where they are in Google... Here's what I got by visualising it.

Points and polygons

Now there are two more issues. First, the difference between geometry and geography. Second, distance calculations.

st_distance on geographies returns values in metres. That's why the output of your first query is basically what you'd expect, given the location of the features (because your points were geographies, Postgres cast the polygons to geography for you). However, st_distance on geometries returns values in the units of the underlying projection, in this case degrees. That explains the output of your second query.

On to distance calculation. Trying to be concise here, but the biggest thing affecting your calculation is the underlying projection. From my understanding, the geography type is better for larger areas - you'd get more accurate results using geometries in a suitable projection for the area you're analysing. The "suitable projection" depends on what parts of the world you're looking at and what your purpose is.

Also, don't assume Google's figures are correct - I did a quick search but didn't find anything on how they calculate their distances. If they use a different method you'll get a different answer...

EDIT Here is the query with the coordinates as they appear in your screenshots. Note the differences in the coordinate values! Also, note that the degrees-minutes-seconds displayed in your screenshots do not match with the decimal degree values in brackets!

SELECT ST_Distance(Place1, Place2) As Place1ToPlace2
    , ST_Distance(Place1, Spot1) As Place1ToSpot1
    , ST_Distance(Place1, Spot2) As Place1ToSpot2
    , ST_Distance(Place2, Spot1) As Place2ToSpot1
    , ST_Distance(Place2, Spot2) As Place2ToSpot2
FROM (SELECT
    ST_PolygonFromText('SRID=4326;POLYGON((-74.0050636293915 40.75265123968514,-74.00500355126653 40.75268991743845,-74.00498169169283 40.75267084386348,-74.00503571044075 40.75263867886528,-74.0050636293915 40.75265123968514))') as Spot1
    ,ST_PolygonFromText('SRID=4326;POLYGON((-74.00503571044075 40.75263867886528,-74.00498225451273 40.75267084385684,-74.00495878551709 40.75265859837483,-74.00501023946696 40.75262521978885,-74.00503571044075 40.75263867886528))') as Spot2
    ,ST_GeogFromText('SRID=4326;POINT( -74.005019 40.752717)') As Place1
    ,ST_GeogFromText('SRID=4326;POINT(-74.004917 40.752650)') As Place2
    ) As foo  ;

The values returned are 11.38223433, 3.27827391, 5.99175215, 5.93327383, 3.65564537. They are within cm of the results from Google.

Upvotes: 1

Tim Biegeleisen
Tim Biegeleisen

Reputation: 521239

I potentially major problem I see is that you are mixing geometries and geographies in your call to ST_Distance. From a glance at Postgres' API, the ST_Distance function takes either geometries or geographies, but not both.

Please try the following query. Here, I have replaced your calls to ST_GeogFromText, which returns a geography, with ST_GeometryFromText, which returns a geometry.

SELECT ST_Distance(Place1, Place2) AS Place1ToPlace2,
       ST_Distance(Place1, Spot1)  AS Place1ToSpot1
       ST_Distance(Place1, Spot2)  AS Place1ToSpot2
       ST_Distance(Place2, Spot1)  AS Place2ToSpot1
       ST_Distance(Place2, Spot2)  AS Place2ToSpot2
FROM
(
    SELECT ST_PolygonFromText('SRID=4326;POLYGON((-74.0050636293915 40.75265123968514,-74.00500355126653 40.75268991743845,-74.00498169169283 40.75267084386348,-74.00503571044075 40.75263867886528,-74.0050636293915 40.75265123968514))') AS Spot1,
           ST_PolygonFromText('SRID=4326;POLYGON((-74.00503571044075 40.75263867886528,-74.00498225451273 40.75267084385684,-74.00495878551709 40.75265859837483,-74.00501023946696 40.75262521978885,-74.00503571044075 40.75263867886528))') AS Spot2,
           ST_GeometryFromText('SRID=4326;POINT(-74.00489 40.752894)') AS Place1
           ST_GeometryFromText('SRID=4326;POINT(-74.004774 40.752846)') AS Place2
) AS foo;

Upvotes: 0

Related Questions