Choosing coordinate system to store geography data for celestial coordinates

Question

I am doing an astronomy project. I want to have the information about our images stored in a spatially enabled database. This, I would think, should be a very easy special case for GIS functions because the sky can be treated as perfectly spherical, and does not require an elliptical treatment like the surface of the earth. Unfortunately, I have not found a way to do this yet and I have been dodging mines with spatial functions that use an elliptical earth. (Pretty much any function that returns meters instead of degrees may be using an elliptical calculation. Luckily, many of the PostGIS functions I've needed appear to have incomplete implementations where the documentation explicitly states that the returned results are for a sphere and not for the ellipsoid. But that may change with future versions, which is a cause for concern.)

Background: I am presently using PostgreSQL with PostGIS and WGS 84 coordinates (SRID=4326). This works fairly well. I am creating a closed POLYGON from the right ascension and declination of the four corners of the image. I have a lot of images (10k or more), covering a large area of the sky. Each images is about 1 degree square. From the set of these images, I am making mosaics from small subsets of 15 to 30 images. Each mosaic is about 1.5 degree square.

Presently, I am storing the geography of the mosaics as a MULTIPOLYGON that consists of all the POLYGONS corresponding to each image that went into the mosaic. [A better solution would be to create a single POLYGON that describes the perimeter of the union of all the individual polygons. I don't know if this can be done in spherical coordinates (i.e., that the geography type). This would also be an interesting answer for me too.] The date line and celestial poles may be included in an image in the dataset so I have been avoiding projecting to planar coordinates to the extent possible.

What coordinate system ought I use for celestial coordinates with PostGIS functions?

I have looked at http://spatialreference.org/ but have not found anything so far. Google has turned up little. I am stumped. Basically, I want to ensure that if a function returns meters as a distance, it is meters along a great circle on a sphere.

---

I am using PostGIS 1.5.2. I have not yet tried PostGIS 2.0. I am curious if the ST_CoveredBy function works with a POLYGON and a MULTIPOLYGON of type geography. If anybody is running 2.0, could you tell me if you get the same error as this:

mydb=# select ST_CoveredBy(ST_GeographyFromText('MULTIPOLYGON(( (10.37795 -69.57926,8.9498 -69.54875,9.0178 -69.21643,10.4242 -69.24648,10.37795 -69.57926),(10.42436 -69.24618,9.01774 -69.2162,     9.08363 -68.88389,10.46914 -68.91344,10.42436 -69.24618)))'),ST_GeographyFromText('POLYGON((10.46915 -68.91315,9.08371 -68.88364,9.14755 -68.5513,10.5125 -68.58038,10.46915 -68.91315))'));
ERROR:  geography_covers: only POLYGON and POINT types are currently supported
CONTEXT:  SQL function "st_coveredby" statement 1

---

I have tried PostGIS 2.0. This function still only works on points and polygons, not more general shapes.

Answer 1

Check out pgsphere, it's specifically designed for handling astronomical data.

http://pgsphere.projects.postgresql.org/

Answer 2

It is possible to store celestial positions in PostGIS - you just need to create your own coordinate system!

PostGIS gets all its coordinate system and projection information from the table spatial_ref_sys which is normally populated when the database is initialized. But there is nothing stopping you adding your own projections - indeed it is practically encouraged.

In common with almost every GIS/spatial database/mapping product out there, PostGIS uses Proj4 for its projection needs, and so you need to put a Proj4 string into the spatial_ref_sys table. A simple spherical SRS in Proj4 form is: +proj=longlat +ellps=sphere +no_defs. PostGIS also requires a WKT version of the projection, but I think that's just used as pretty text.

You will also need to come up with a unique SRID for your new SRS, as well as an "authority", but this can be anything you like.

So to insert a new entry into spatial_ref_sys, just perform this SQL:

insert into spatial_ref_sys values(40000, 'ME', 1, 
'GEOGCS["Normal Sphere (r=6370997)",DATUM["unknown",SPHEROID["sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]',
'+proj=longlat +ellps=sphere +no_defs');

Note that I have chosen 40000 as the SRID - this is the number you use in your celestial object table. The authroity is "ME", but this could be your name, organization, or anything really up to 256 characters. The next number, 1, is just your unique identifier for that entry, relative to the authority. In theory you could refer to this entry as ME:1, but for all PostGIS processing, its the unique SRID that counts. The WKT entry I generated with GDAL and Python:

import osgeo.osr as osr
srs = osr.SpatialReference()
srs.ImportFromProj4('+proj=longlat +ellps=sphere +no_defs')
srs.ExportToWkt()

Now the caveats:

• Right ascension will have to be specified in degrees rather than hour-angles.
• A number of PostGIS functions aren't designed for unprojected data, but it's the same issue if you have terrestrial data in WGS84 long/lat.
• As it stands, the data is geocentric. If you want to do any observational work with it, I suggest using something like PyEphem.
• I've not tried creating any data in this SRS, so YMMV.
• I'm quite interested in this now though, so I may have to play around with importing the Hipparchos catalogue... :)