1

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

2
  • you could try breaking the calculation out into separate steps to see which step is failing, or switch to using a geographic library like Shapely to handle the maths for you.
    – Ian Turton
    Feb 17, 2023 at 17:08
  • Alas, I can't, as all of this logic is passed inside an annotation in my Django queryset
    – Ash_Crow
    Feb 17, 2023 at 17:43

2 Answers 2

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)
1
  • Thanks, that was helpful, because it confirmed the intuition I had after posting the question: it is a float rounding error and depending on the system it tries to calculate the acos of something slightly inferior to 1 (your result) or slightly superior (mine, which results in a math error).
    – Ash_Crow
    Feb 20, 2023 at 9:09
0

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")
)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.