Math error when trying to calculate the distance between two coordinates in python

Question

I'm not sure if this is the right place but I am trying to develop a Django site that (among other things) shows public projects within 30 km of a given French commune. It mostly works fine but for some municipalities (for example Marseille) I have an error when trying to display the page, which after investigation seems to be related not to Django or PostgreSQL but to the actual coordinates I have in the database (those coordinates are from a public dataset provided by the French government)

As you can see below, I tried to isolate the problem : it occurs when trying to do the distance calculation for the projects located inside Marseille (line 172). If I change the latitude by a few meters (line 173), it is now able to correctly calculate that the distance is zero:

In [166]: paris.longitude
Out[166]: 2.347

In [167]: paris.latitude
Out[167]: 48.8589

In [168]: marseille.longitude = 5.3806

In [169]: marseille.latitude = 43.2803

In [170]: marseille.save()

In [171]: acos(cos(radians(marseille.latitude))
     ...:     * cos(radians(paris.latitude))
     ...:     * cos(radians(paris.longitude) - radians(marseille.longitude))
     ...:     + sin(radians(marseille.latitude))
     ...:     * sin(radians(paris.latitude))
     ...: ) * 6371
Out[171]: 662.847013493097

In [172]: acos(cos(radians(marseille.latitude))
     ...:     * cos(radians(marseille.latitude))
     ...:     * cos(radians(marseille.longitude) - radians(marseille.longitude))
     ...:     + sin(radians(marseille.latitude))
     ...:     * sin(radians(marseille.latitude))
     ...: ) * 6371
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[172], line 1
----> 1 acos(cos(radians(marseille.latitude))
      2    * cos(radians(marseille.latitude))
      3    * cos(radians(marseille.longitude) - radians(marseille.longitude))
      4    + sin(radians(marseille.latitude))
      5    * sin(radians(marseille.latitude))
      6 ) * 6371

ValueError: math domain error

In [173]: marseille.latitude = 43.2804

In [174]: marseille.save()

In [175]: acos(cos(radians(marseille.latitude))
     ...:     * cos(radians(marseille.latitude))
     ...:     * cos(radians(marseille.longitude) - radians(marseille.longitude))
     ...:     + sin(radians(marseille.latitude))
     ...:     * sin(radians(marseille.latitude))
     ...: ) * 6371
Out[175]: 0.0

I'm not a mathematician, nor used to manipulate geographic information. What is the cause of the problem? Can the same error occur when calculating the distance between two different cities? Can I fix it without having to manually change the coordinates for an unknown number of the 36000 French municipalities?

Note : it is a simplification of the actual code used, which is here: https://github.com/MTES-MCT/aides-territoires/blob/7742d9f69aa60255e674542b9faa595b61b3ed20/src/projects/forms.py#L441

Answer 1

Just running this in python it seems to calculate the distance without error, though the marseille to marseille distance is 0.0013 rather than 0.

import numpy
from math import acos, cos, sin, radians

paris = [48.8589, 2.347]

marseille = [43.2803, 5.3806]


p2m = acos(cos(radians(marseille[0]))
     * cos(radians(paris[0]))
     * cos(radians(paris[1]) - radians(marseille[1]))
     + sin(radians(marseille[0]))
     * sin(radians(paris[0]))
     )* 6371
print(p2m)

m2m = acos(cos(radians(marseille[0]))
    * cos(radians(marseille[0]))
    * cos(radians(marseille[1]) - radians(marseille[1]))
    + sin(radians(marseille[0]))
    * sin(radians(marseille[0]))
    ) * 6371
print(m2m)

Answer 2

So, after more checks (and a bit of arithmetic), calculating the distance between two identical points should simplify to calculating the arccosine of 1, which is 0.

Due to a float rounding error, and depending on the system, it sometimes tries to calculates the arccosine of something like 0.999997 or 1.000003 instead of 1, which in turn causes a math error if the value is >1.

Fortunately, Django database functions have a way to cap the value, so I solved the issue this way:

from django.db.models import F
from django.db.models.functions import ACos, Cos, Radians, Sin, Round, Least

qs = (
    qs.annotate(
        distance=Round(
            ACos(
                Least(
                    Cos(Radians(search_perimeter.latitude))
                    * Cos(Radians(F("organization__perimeter__latitude")))
                    * Cos(
                        Radians(F("organization__perimeter__longitude"))
                        - Radians(search_perimeter.longitude)
                    )
                    + Sin(Radians(search_perimeter.latitude))
                    * Sin(Radians(F("organization__perimeter__latitude"))),
                    1.0,
                )
            )
            * 6371,
            precision=2,
        )
    )
    .filter(distance__lte=30)
    .order_by("distance")
)