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.

  • Isn't this similar to what wcs2kml does ? If so, maybe you could adapt some of the code for your uses. code.google.com/p/wcs2kml Commented Oct 9, 2010 at 16:24
  • I ran into this USGS presentation, "PLANETARY GIS 101" and briefly saw there's a few slides on projections, maybe It'll help you out.
    – jonatr
    Commented Oct 29, 2010 at 6:31
  • Rather than creating multipolygons, why not create multiple polygons that share a grouping id?
    – raphael
    Commented Oct 5, 2014 at 18:13

2 Answers 2


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


  • This is very good stuff. Unfortunately, it doesn't seem to support any "smultipoly" geometry class. Thanks for the great heads-up on this project though.
    – SO Stinks
    Commented Oct 19, 2010 at 5:55
  • 2
    When I try to follow this link I get "Forbidden You don't have permission to access / on this server."
    – PolyGeo
    Commented Dec 15, 2018 at 6:24

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

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... :)
  • 2
    +1. You can get a good start on mapping the skies by loading a version of the HYG database, multiplying the right ascension by 15 and subtracting 180 to convert to a standard GIS "longitude," and using any perfectly spherical datum you like. For display and mapping, gnomonic and orthographic projections are fairly standard.
    – whuber
    Commented Jan 27, 2014 at 18:45
  • @whuber : and the latitude? is the DEC?
    – Magno C
    Commented Jan 29, 2014 at 16:09
  • @MagnoC Yes, that's correct. The fields are described on the website I linked to: just scroll down a little. To check it I dumped the "small" version (only 31K stars) into a 3D viewing program, converted to cartesian coordinates (on a unit celestial sphere, ignoring distance), and plotted them: looks good.
    – whuber
    Commented Jan 29, 2014 at 16:11
  • @whuber :"RA, Dec: The star's right ascension and declination, for epoch 2000.0. Stars present only in the Gliese Catalog, which uses 1950.0 coordinates, have had these coordinates precessed to 2000". I'ts not so clear about Lat/Lon. So, LON = (RA*15) - 180 and LAT = DEC ?
    – Magno C
    Commented Jan 29, 2014 at 16:17
  • @MagnoC I have found en.wikipedia.org/wiki/Equatorial_coordinate_system to be helpful in sorting this out.
    – whuber
    Commented Jan 29, 2014 at 16:19

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.