How to sort for non-geographical (X, Y, Z) distances using GeoDjango / PostGIS?

Question

I'm currently using GeoDjango in a "Star searching and sorting" database that provides information about and will simulate information on star and planetary systems.

I'm using GeoDjango 1) because I like it and use it elsewhere, and 2) because I eventually want to use various "geo" searching features like distance/lines/polygon transforms for complex and cool volumetric querying that I can't find elsewhere.

I have the system up and running (https://github.com/jaycrossler/procyon) and it currently uses star Galactic Coordinates (already transformed from Right Ascension/Declination). There are 150k stars currently in the database, and I'm considering increasing to a few million.

After a star is in the database, I build a new table that has a GeoDjango PointField, then populate it with the Galactic X, Y, Z coordinates (which are in parsecs, and such are mostly in the range of -500 to 500). Right now, I've set the SRID to 900913 (so that I'll have a good range of coordinates that won't roll over around the world line)... but when I search on the nearby stars and order by distance, I'm only getting returns that are in a line, rather than being truly near based on X,Y,Z distance.

location = models.PointField(dim=3, blank=True, null=True, srid=900913)
objects = models.GeoManager()

I think this is because every search I make is ultimately getting wrapped onto the surface of a sphere, and that's both inefficient and screwing with results (though if it makes searching only take one line of code, I'm cool with it).

The current searching I'm using in Django is:

    origin = self.location

    distance = 1000
    close_by_stars = StarModel.objects.filter(location__distance_lte=(origin, D(m=distance))).distance(origin).order_by('distance')
    for s in close_by_stars[:200]:
      #export results

But the results returned are not what I expect (I would think they'd clump around one star, not be in a line), visualized:

So, the big question is: 1) Should I use a SRID such as 900913 (Spherical Mercator) or 2) Is there a SRID that isn't mapped to the surface of a a planet, so that I can just search on X, Y, Z distances without it rolling over a -180 into a +180 (or whatever equivalent based on projection system)? I tried using SRID=0, but geodjango than pukes and doesn't allow that.

Answer 1

I've got a fix, which I'm sharing here to potentially help others.

I think the issue is that the 'location__distance_lte=' function gets transformed in POSTGIS to the geo function ': ST_DWithin

From looking at pg_log/latest, I see the SQL command:

SELECT (ST_Distance("modelname"."location",
ST_GeomFromEWKB('\x01010000a031bf0d0021e527d53e2963409f3c2cd49ab2404081b22957787d5540'::bytea))) AS "distance",
"modelname"."info", "modelname"."location" 
FROM "modelname" WHERE ST_Distance("modelname"."location",
ST_GeomFromEWKB('\x01010000a031bf0d0021e527d53e2963409f3c2cd49ab2404081b22957787d5540'::bytea))
<= 10.0 ORDER BY "distance" ASC

So, when searching by X, Y, Z and looking for the nearest points, it only searches within 2D space - and looks for the ones within X, Y distance... not the ones within Z.

There is an ST_3DDWithin (http://postgis.net/docs/ST_3DDWithin.html) but unfortunately Django doesn't know about it: https://github.com/django/django/blob/master/django/contrib/gis/db/backends/postgis/operations.py#L154.

Instead of overriding the django source I could use the raw sql method: https://docs.djangoproject.com/en/dev/topics/db/sql/#django.db.models.Manager.raw. But, then the orm would lose a lot of it's benefit.

But instead, I decided to go a bit simpler/more complicated. I kept the same search (which returns basically a "cylinder" instead of a sphere of results). Then, in the function, I loop through the results and parse out the ones that aren't within the spherical results:

origin = item.location
origin_array = numpy.array((origin.x, origin.y, origin.z))

close_by_stars = star_model.objects.filter(location__distance_lte=(origin, D(m=distance))).distance(origin).order_by('distance')

star_list = []
for s in close_by_stars:
    location_array = numpy.array((s.location.x, s.location.y, s.location.z))

    dist = numpy.linalg.norm(origin_array - location_array)
    if dist > distance:
        continue

    star_handle = dict()
    star_handle['data'] = s.data ...