Getting started with GEOS

bpolaszek commented 3 months ago

Hello there,

Sorry for asking for support here, but in an attempt to get rid of MySQL I/Os, I wanted to give GEOS a try but I get very unexpected results.

Here's an example distance from Biarritz (FR) to Strasbourg (FR):

use Brick\Geo\Engine\GEOSEngine;
use Brick\Geo\Engine\PDOEngine;
use Brick\Geo\Point;
use PDO;

$pdoEngine = new PDOEngine(new PDO('mysql:host=localhost'));
$geosEngine = new GEOSEngine();
$from = Point::fromText('POINT(-1.50331 43.48997)', 4326); // Biarritz Lon+Lat
$to = Point::fromText('POINT(7.74420 48.58754)', 4326); // Strasbourg Lon+Lat
$dbDistance = $pdoEngine->distance($from, $to); // 1168767.43 meters -> 1168 km, seems accurate
$geosDistance = $geosEngine->distance($from, $to); // 10.56 meters => 0.01 km 🤨

Is there something I didn't get in the first place?

Thank you! Ben

nono303 commented 3 months ago

Hi @bpolaszek See #29

bpolaszek commented 3 months ago

Hi @nono303,

Thanks for your answer. So, if I understand correctly, GEOSEngine returns a distance in degrees, whereas MySQL returns a distance in meters. I'm very bad at maths (I just need a fast and reliable way to calculate a distance in meters between 2 points on Earth) so I asked ChatGPT to provide me with a function to translate those degrees to a distance in meters, and it came up with:

function degreesToRadians(float $degrees): float {
    return $degrees * M_PI / 180.0;
}

function degreesToMeters(float $degrees, float $latitude): float {
    $earthRadius = 6371000; // in meters
    // Convert latitude to radians
    $latInRadians = degreesToRadians($latitude);
    // Approximate distance in meters for 1 degree latitude/longitude
    $metersPerDegree = $earthRadius * cos($latInRadians) * M_PI / 180;
    return $degrees * $metersPerDegree;
}

$midLatitude = ($from->y() + $to->y()) / 2;
$geosDistanceInMeters = degreesToMeters($geosDistance, $midLatitude); // 815065.50511093

This gives us ~815km VS ~1168km with MySQL, which is a pretty big difference. I assume that calculations are done in different ways and that accuracy may vary, but I would expect a difference of a few kilometers, not hundreds. Is there something I could do to optimize that? Did I forget or misuse something, somewhere?

Thank you, Ben

nono303 commented 3 months ago

here is my "dirty" Distance.php class, you may find some interest

BenMorel commented 3 months ago

Hi @bpolaszek,

So, if I understand correctly, GEOSEngine returns a distance in degrees, whereas MySQL returns a distance in meters.

Correct.

I'm very bad at maths (I just need a fast and reliable way to calculate a distance in meters between 2 points on Earth) so I asked ChatGPT to provide me with a function to translate those degrees to a distance in meters, and it came up with: (...) This gives us ~815km VS ~1168km with MySQL, which is a pretty big difference. I assume that calculations are done in different ways and that accuracy may vary, but I would expect a difference of a few kilometers, not hundreds. Is there something I could do to optimize that? Did I forget or misuse something, somewhere?

The calculation ChatGPT gave you returns a number of a meters per degree on the x-axis (lon), given a latitude. However, you have another number of meters per degree on the y-axis (lat), which does not depend on the longitude.

A function for calculating this one would be:

function degreesToMeters(float $degrees): float {
    $earthRadius = 6371000;
    $metersPerDegree = $earthRadius * M_PI / 180;
    return $degrees * $metersPerDegree;
}

Which would return 1174 km between your 2 points. As you can see, the actual result (~1168 km~ actually 911 km, see my next post below), lies somewhere between the numbers returned by the two functions.

And as you can see, attempting to convert degrees to meters is bound to failure.

I think your options are:

if you only need distances between 2 points, you can calculate them using the Haversine or Vincenty's formulas; the PHP script uploaded by @nono303 above contains implementations of them (I didn't check them, though)
if you need more complicated calculations, but on geometries in well defined geographical areas, you can project your points on a flat plane whose coordinates are in meters (I think this is the case for Lambert 93 in France, but I would have to check)
finally, you can keep using MySQL :confused:

bpolaszek commented 3 months ago

Hello @BenMorel @nono303,

Alright, thanks for your explanations! Both Haversine & Vincenty's formulaes return 911km, which is still far from what MySQL computes, but I'll cope with this.

Sorry for bothering! Have a nice day, Ben

BenMorel commented 3 months ago

@bpolaszek Actually Haversine & Vincenty's formulas are right, the actual distance is 911 km: https://boulter.com/gps/distance/?from=43.48997%2C-1.50331&to=48.58754%2C7.74420&units=k

If you invert lat & lon though, you get 1168 km.

Indeed, MySQL uses the ESPG dataset, which defines SRID 4326 as lat/lon, not lon/lat: https://dev.mysql.com/blog-archive/axis-order-in-spatial-reference-systems/

This code gives the correct distance:

$pdoEngine = new PDOEngine(new PDO('mysql:host=localhost'));
$from = Point::xy(43.48997, -1.50331, 4326); // Biarritz Lat+Lon
$to = Point::xy(48.58754, 7.74420, 4326); // Strasbourg Lat+Lon
$dbDistance = $pdoEngine->distance($from, $to); // 911847

TL;DR: you need to swap lat and lon.

Also, I checked the Lambert 93 projection, and indeed results are in meters.

If I convert your GPS coordinates into Lambert 93, I get:

Biarritz: (-1.50331, 43.48997) => (335771.674, 6276042.545)
Strasbourg: (7.74420, 48.58754) => (1049686.559, 6842433.713)

Then, GEOS can give you the correct answer:

$geosEngine = new GEOSEngine();
$from = Point::xy(335771.674, 6276042.545, 2154); // Biarritz Lambert 93
$to = Point::xy(1049686.559, 6842433.713, 2154); // Strasbourg Lambert 93
$distance = $geosEngine->distance($from, $to); // 911303

Unfortunately, php-geos does not support transform(), so you cannot use it to convert GPS coordinates to Lambert 93 coordinates, and need to find another way to convert coordinates before storing them!

bpolaszek commented 3 months ago

Hi @BenMorel,

Thanks a lot for following up on that topic! That's both fascinating and confusing 😅 On my own MySQL implementation, I correctly got 910389m but I forgot that when I switched to Brick\Geo\Engine\PDOEngine. Somehow the engine inverted Lat and Lon where it should not (IMO), as WKT POINT(a b) is supposed to be POINT(Lon Lat) as I did.

I'll cope with Vincenty's formula PHP implementation @nono303 provided!

Thank you very very much for your precious help! 🙏 Ben

brick / geo

Getting started with GEOS #50