[Proj] proj4 string for perfect sphere

Frank Warmerdam warmerdam at pobox.com
Fri Apr 13 18:26:02 EST 2012


On Fri, Apr 13, 2012 at 3:57 PM, Demitri Muna
<thatsanicehatyouhave at me.com> wrote:
> Hi,
>
> I'm an astronomer and I'm exploring the possibility of using GIS tools (specifically spatialite at the moment) with astronomical objects. I'm not conversant in GIS at the moment (or really the whole field), but I've made a little progress.
>
> It seems that I need an appropriate SRID for what I'm doing. I want a perfect sphere, I'm working exclusively in degrees (a lat/lon analogy to ra/dec seems fine), and all of my points can be considered on the surface of the sphere. What SRID would be appropriate for this description?
>
> I came up with this (probably laughably simplistic):
>
> +proj=moll +ellps=sphere

Demitra,

I'm not not sure why you are using +proj=moll (Mollweide) if you are just
going to be working with lat/long degrees.  Instead you should use latlong.

eg.
+proj=latlong +ellps=sphere

However, if you want an existing SRID, you could use EPSG:4035 which is:

warmerdam at gdal:~$ gdalsrsinfo EPSG:4035
PROJ.4 : '+proj=longlat +a=6371000 +b=6371000 +no_defs '

OGC WKT :
GEOGCS["Unknown datum based upon the Authalic Sphere",
    DATUM["Not_specified_based_on_Authalic_Sphere",
        SPHEROID["Sphere",6371000,0,
            AUTHORITY["EPSG","7035"]],
        AUTHORITY["EPSG","6035"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9108"]],
    AUTHORITY["EPSG","4035"]]

> In spatialite using that SRID (or none at all for that matter), I did this:
>
> CREATE TABLE data (
> id INTEGER PRIMARY KEY,
> name TEXT NOT NULL,
> coordinate BLOB
> );
>
> INSERT INTO data (name, coordinate) VALUES ('NGC23',GeomFromText('POINT(2.47254166667 25.9237777778)'));
> INSERT INTO data (name, coordinate) VALUES ('NGC17',GeomFromText('POINT(2.77729166667 -12.1073136111)'));
>
> These are degree values, and when I did this:
>
> SELECT name, DISTANCE( GeomFromText('POINT(2.77729166667 -12.1073136111)'), coordinate ) FROM data;
>
> I get:
>
> NGC23       38.0323123776791
> NGC17       0.0
>
> which apparently is the distance on a flat x-y cartesian plane, not the angular (or great circle) distance I'm looking for (38.0322481853 degrees).

If you want great circle distances (geodescics) you would need
something apart from the default distance operator.   If you just
wanted to do it at the commandline you could use the geod command from
PROJ.4.  If you want to do it in a database, I'd normally recommend
using the special "GEOGRAPHY" geometry type in PostGIS which is
oriented to calculations on a sphere, particularly using geodesics.

I'm not sure if there is an appropriate function in spatialite.

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Software Developer


More information about the Proj mailing list