proj4 string for perfect sphere

classic Classic list List threaded Threaded
16 messages Options
Reply | Threaded
Open this post in threaded view
|

proj4 string for perfect sphere

Demitri Muna
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

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

Is the wrong SRID the problem in the first place?

Cheers,
Demitri

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

Frank Warmerdam
On Fri, Apr 13, 2012 at 3:57 PM, Demitri Muna
<[hidden email]> 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@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, [hidden email]
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Software Developer
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

Demitri Muna
Hi Frank,

Thanks very much for the help!

On Apr 13, 2012, at 7:26 PM, Frank Warmerdam wrote:

> 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.

Ah, I didn't see that option; that makes much more sense. I chose Mollweide when I came across this:
http://zerospace.org/gis/projections/mollweide.html

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

Ah, that one looks right!

> 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.

Let me step back a bit and explain what I am trying to do (I became focussed on the SRID detail!). I have two lists of coordinates (ra/dec in my field, but lat/lon on a perfect sphere is the same). I'm trying to find out how many of the points from list 1 matches those on list 2 to within some value, e.g. 5 arcsec. I want to load each list into [x], and for each object in the first list perform a search to see if it appears in the second list. I was hoping to take advantage of GIS spatial indexing to make this a fast operation - I'm guessing this is a simple operation in the field.

[x] above is some GIS tool that I'm trying to find. I'm a huge fan of PostgreSQL (and will at some point investigate PostGIS), but my first goal is to do this with a Python script for catalogs < ~100,000. If I have two simple lists, I don't want to have to go into PG, create the schemas, populate, etc. spatialite seemed perfect - create the db file on the fly, populate it from the script, do the search with the GIS tools, then throw the file away. Simple. (At least that's the idea.) It also can be used without requiring someone to install PG+PostGIS, and I can wrap the spatialite complexity in Python classes. If this can be done with another (e.g. command line) tool I'd love to hear about it.

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

I was using the distance operator to test the geometry, which is how I got here. spatialite does seem to have a function:

GreatCircleLength( c Curve ) : Double precision

which is a little unclear to me how to use:

spatialite> GreatCircleLength(GeomFromText('LINESTRING(2.47254166667 25.9237777778, 2.77729166667 -12.1073136111)',4035))
---------------------------------------------------------------------------------------------                
4228993.04708947                                                                                            

Those are not degrees! :)

Thanks again for your help. I'm trying to build a bridge between the astronomy and GIS communities...

Cheers,
Demitri
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

Demitri Muna
Hi,

Quick update. Once I recognized 6371000 as the radius of the Earth in meters, I created my own SRID with these parameters (r = 1 radian as degrees):

"+proj=longlat +a=57.2957795 +b=57.2957795 +no_defs"

This worked perfectly:

> SELECT GreatCircleLength(GeomFromText('LINESTRING(2.47254166667 25.9237777778, 2.77729166667 -12.1073136111)',50000));
38.0322481765926

The correct answer! I would still appreciate your thoughts/advice on what I'm trying to do more generally.

Cheers,
Demitri
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

Fischer, Robert P. (GISS-6110)[COLUMBIA UNIVERSITY]
My thoughts:

1. You can compute distance between two lat/lon points on a sphere using the Haversine formula or Vincenty's formula specialized to the sphere.  They are described at:
   http://en.wikipedia.org/wiki/Great-circle_distance
Code for Haversine is below.  Fix it up as needed.

2. If all you need is an R-Tree, just use an R-Tree.  Don't bring in a database just because it has an R-Tree built inside, that will distort your entire programming effort needlessly.  Here is an R-Tree that works in 2 or 3 dimensions.  It worked for me in just a couple of hours of work:
   https://github.com/nushoin/RTree
NOTE: There are multiple versions of this file floating around the Internet.  Use what works for you.

3. Your real problem is that you want an R-Tree on the surface of a sphere, not on a plane.  Projecting your points to a plane and then using an R-Tree will only work if the points are within a localized area (and you have a good bound on how much your distances might get distored).  The reason it's not a global solution is that any projection has certain areas of the sphere where nearby points get projected to far-away points on the plane.  Because a sphere is periodic and a plane is not.

There are two ways to solve this:
 a) Develop an R-Tree algorithm on the sphere --- i.e. for a periodic type of geometry.  This would involve some work, and there's no obvious evidence that anyone has done it.  But I believe it would work.

 b) Use a 3-D R-Tree, based on the Cartesian coordinates of your points in 3-D.  That approach is described here, and apparently works:
   http://blog.cleverelephant.ca/2009/11/postgis-gets-spherical-directors-cut.html


I quote: "How does this magic index work? The "trick" is to be as lazy as possible. I didn't want to write a whole new indexing scheme, and I already had a serviceable R-Tree index. What I needed to do was convert the lat/lon coordinates to a domain where the problems of the singularities at the poles and dateline would go away. The answer is to convert from spherical coordinates (lat/lon) relative to Greenwich into geocentric coordinates (x/y/z) relative to the center of earth. It's easy then to build a 3D R-Tree on the geocentric bounds of the features. Calculating the bounds in 3D is tricky, because the curvature of the features over the spherical surface must be respected, but once that is done, the index performs admirably."

Unfortunately, the author (Paul Ramsey) does not explain how he calculates his bounds in 3D.  I would pester him and find out.  Would be interested in hearing back if you do.

Either way, projections aren't needed or involved.  Hence, you shouldn't be using proj.4.  Using proj.4 on a "no-op" spherical projection doesn't make mush sense.  If all you want is a distance formula, that's really easy without proj.4.  Nor does a database really buy you much either, as far as I can tell.

Or, you could use the GEOGRAPHY data type that Paul Ramsey wrote about implementing in PostgreSQL:
    http://blog.opengeo.org/2009/11/04/postgis-gets-spherical/

If you end up using PostgreSQL in this way, I would not look at is a black box.  I would instead look at it as a convenient package of algorithms, already coded up, that you basically understand.

-- Bob


============================================
D2R = M_PI / 180
inline double sqr(double x) { return x*x; }

/** See: http://www.cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
@return Distance (in degrees) */
extern double haversine_distance(
double lon1_deg, double lat1_deg,
double lon2_deg, double lat2_deg)
{
        // Convert inputs to degrees
        double lon1 = lon1_deg * D2R;
        double lat1 = lat1_deg * D2R;
        double lon2 = lon2_deg * D2R;
        double lat2 = lat2_deg * D2R;

        // Apply the Haversine Formula
        double dlon = lon2 - lon1;
        double dlat = lat2 - lat1;
        double a = sqr(sin(dlat/2)) + cos(lat1) * cos(lat2) * sqr(sin(dlon/2));
        double c = 2 * atan2( sqrt(a), sqrt(1-a) );
        return c * R2D; // Convert to degrees
}


________________________________________
From: [hidden email] [[hidden email]] On Behalf Of Demitri Muna [[hidden email]]
Sent: Friday, April 13, 2012 9:24 PM
To: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] proj4 string for perfect sphere

Hi,

Quick update. Once I recognized 6371000 as the radius of the Earth in meters, I created my own SRID with these parameters (r = 1 radian as degrees):

"+proj=longlat +a=57.2957795 +b=57.2957795 +no_defs"

This worked perfectly:

> SELECT GreatCircleLength(GeomFromText('LINESTRING(2.47254166667 25.9237777778, 2.77729166667 -12.1073136111)',50000));
38.0322481765926

The correct answer! I would still appreciate your thoughts/advice on what I'm trying to do more generally.

Cheers,
Demitri
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

Demitri Muna
Hi Bob,

Thanks very much for the extensive comments. The biggest thing I learned about GIS this week (starting from knowing nothing about it last week) is that calculations are primarily done on an XY plane which are projections of a curved surface. I had naively assumed going in that calculations and indices (e.g. spatialite) were optimized for the surface of a sphere.

On Apr 13, 2012, at 10:16 PM, Fischer, Robert P. (GISS-6110)[COLUMBIA UNIVERSITY] wrote:

> 1. You can compute distance between two lat/lon points on a sphere using the Haversine formula or Vincenty's formula specialized to the sphere.  They are described at:
>   http://en.wikipedia.org/wiki/Great-circle_distance
> Code for Haversine is below.  Fix it up as needed.

Thanks for the code. I've been using the law of cosines and wasn't aware of the problem at small separations (though this seems to trade that problem for one at large separations).

> 2. If all you need is an R-Tree, just use an R-Tree.  Don't bring in a database just because it has an R-Tree built inside, that will distort your entire programming effort needlessly.  Here is an R-Tree that works in 2 or 3 dimensions.  It worked for me in just a couple of hours of work:
>   https://github.com/nushoin/RTree
> NOTE: There are multiple versions of this file floating around the Internet.  Use what works for you.

I was not familiar with this - I will look into it. You are right; the database isn't crucial for what I'm doing.

I will also look at the other links you mentioned. Since all my source points are in ra/dec (not 3D cartesian), there will be the computation hit to convert them back and forth. So I agree, proj.4 doesn't give me much.

As for indexing on a sphere, I've used this for a few years now with great success:

http://code.google.com/p/q3c/

Some details on how it works are here (though the first link is a little old):

http://www.sai.msu.su/~megera/wiki/SkyPixelization
http://lnfm1.sai.msu.ru/~math/docs/adass_proceedings2005.pdf

My GIS wanderings were to try to find a similar solution for smaller datasets without the overhead of creating a new schema, importing the data, etc. to take advantage of the optimized indexing. It might be in the end that implementing q3c into sqlite or through a Python interface is the best solution.

Thanks again for your suggestions.

Cheers,
Demitri


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

Fischer, Robert P. (GISS-6110)[COLUMBIA UNIVERSITY]
> As for indexing on a sphere, I've used this for a few years now with
> great success:
>
> http://code.google.com/p/q3c/
>
> My GIS wanderings were to try to find a similar solution for smaller
> datasets without the overhead of creating a new schema, importing the
> data, etc. to take advantage of the optimized indexing. It might be in
> the end that implementing q3c into sqlite or through a Python
> interface is the best solution.

Have you tried creating temporary schemas/tables in PostgreSQL with q3c?  Is that fast enough?  It sounds like that would be by far the easiest solution for you.

> Thanks very much for the extensive comments. The biggest thing I
> learned about GIS this week (starting from knowing nothing about it
> last week) is that calculations are primarily done on an XY plane
> which are projections of a curved surface. I had naively assumed going
> in that calculations and indices (e.g. spatialite) were optimized for
> the surface of a sphere.

Yes, that's what I surmised too.  The good news is, the technique works quite well for local computations, and you often choose whatever projection you like.  For example, area of a general closed curve on the surface of a sphere can be computed this way, using any old equal area projection centered near your polygon.

> > 1. You can compute distance between two lat/lon points on a sphere
> >   using the Haversine formula or Vincenty's formula specialized to
> >   the sphere.  They are described at:
> >   http://en.wikipedia.org/wiki/Great-circle_distance Code for
> >   Haversine is below.  Fix it up as needed.
>
> Thanks for the code. I've been using the law of cosines and wasn't
> aware of the problem at small separations (though this seems to trade
> that problem for one at large separations).

Vincenty's formula solves that problem.

> I will also look at the other links you mentioned. Since all my source
> points are in ra/dec (not 3D cartesian), there will be the computation
> hit to convert them back and forth.

Is this signficant for your problem?
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

support.mn
In reply to this post by Demitri Muna
Hello,

Demitri Muna [[hidden email]] kirjoitti:
> Hi Bob,
>
> Thanks very much for the extensive comments. The biggest thing I learned about GIS this week (starting from knowing nothing about it last week) is that calculations are primarily done on an XY plane which are projections of a curved surface. I had naively assumed going in that calculations and indices (e.g. spatialite) were optimized for the surface of a sphere.
>

A small notice for a beginner: If you want to stay accurate and avoid lots of trouble with projections generated errors never leave the sphere (or ellipsoid) just keep doing all calculations there - it is much easier (you don't need projections) and you will stay accurate for ever ;)

regards: Janne.


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

Mikael Rittri

Janne wrote:

> never leave the sphere (or ellipsoid) just keep doing all calculations
> there - it is much easier (you don't need projections)

"Much easier" are not the words I would use. "Much harder" fits
better. That's one reason projections are useful.

Spherical geometry isn't too bad for simple problems (distance,
azimuth etc.), but just think of point-in-polygon test, area of polygon,
and so on.

Regards,
Mikael Rittri
Carmenta
Sweden
http://www.carmenta.com

-----Original Message-----
From: [hidden email] [mailto:[hidden email]] On Behalf Of [hidden email]
Sent: den 17 april 2012 11:54
To: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] proj4 string for perfect sphere

Hello,

Demitri Muna [[hidden email]] kirjoitti:
> Hi Bob,
>
> Thanks very much for the extensive comments. The biggest thing I learned about GIS this week (starting from knowing nothing about it last week) is that calculations are primarily done on an XY plane which are projections of a curved surface. I had naively assumed going in that calculations and indices (e.g. spatialite) were optimized for the surface of a sphere.
>

A small notice for a beginner: If you want to stay accurate and avoid lots of trouble with projections generated errors never leave the sphere (or ellipsoid) just keep doing all calculations there - it is much easier (you don't need projections) and you will stay accurate for ever ;)

regards: Janne.


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

Fischer, Robert P. (GISS-6110)[COLUMBIA UNIVERSITY]
I think it's hard to draw generalizations here, and these issues should be handled on a case-by-case basis.  The relevant issue is, WHAT ERRORS does working with a projection entail, and HOW BIG are they?  The answer varies for different problems.  In some cases, working with a projection is 100% accurate.  In that case, the projection can be seen as a way to derive computational geometry algorithms on a sphere.

As an aside, I was really surprised to see how little support CGAL offered for anything on a sphere.  Although the theory for all this is probably sitting in a geometry book somewhere (I even have an idea of where), the on-the-ground practice really lags behind.

-- Bob
________________________________________
From: [hidden email] [[hidden email]] On Behalf Of Mikael Rittri [[hidden email]]
Sent: Tuesday, April 17, 2012 7:12 AM
To: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] proj4 string for perfect sphere

Janne wrote:

> never leave the sphere (or ellipsoid) just keep doing all calculations
> there - it is much easier (you don't need projections)

"Much easier" are not the words I would use. "Much harder" fits
better. That's one reason projections are useful.

Spherical geometry isn't too bad for simple problems (distance,
azimuth etc.), but just think of point-in-polygon test, area of polygon,
and so on.

Regards,
Mikael Rittri
Carmenta
Sweden
http://www.carmenta.com
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

support.mn
In reply to this post by Demitri Muna
if you need any accuracy and robustness you'll never leave the
sphere (or ellipsoid). It is much easier to spend a little more time
with the general spherical calculation than to fix it afterwards for
150+ projections for each one with maybe a different method.

If you are afraid of the ellipsoid just use the sphere and add some
correction tables or similar methods since the difference is
usually very very small. Another trick is to project it to a single
standard projection which works all over the world and then
do the plane calculations if they are too hard to be done on the
sphere.

Nobody likes the situation where the calculations are made on any
distorted projection plane .. the error can be huge in some cases!
At least the routine should be robust to detect such distortions.. but
this adds a lot of programming effort.. so it is still better not to leave
the sphere (or ellipsoid) which is the general case and works
everywhere and with every projection.

regards: Janne.

--------------------------

Mikael Rittri [[hidden email]] kirjoitti:

>
> Janne wrote:
>
> > never leave the sphere (or ellipsoid) just keep doing all calculations
> > there - it is much easier (you don't need projections)
>
> "Much easier" are not the words I would use. "Much harder" fits
> better. That's one reason projections are useful.
>
> Spherical geometry isn't too bad for simple problems (distance,
> azimuth etc.), but just think of point-in-polygon test, area of polygon,
> and so on.
>
> Regards,
> Mikael Rittri
> Carmenta
> Sweden
> http://www.carmenta.com
>
> -----Original Message-----
> From: [hidden email] [mailto:[hidden email]] On Behalf Of [hidden email]
> Sent: den 17 april 2012 11:54
> To: PROJ.4 and general Projections Discussions
> Subject: Re: [Proj] proj4 string for perfect sphere
>
> Hello,
>
> Demitri Muna [[hidden email]] kirjoitti:
> > Hi Bob,
> >
> > Thanks very much for the extensive comments. The biggest thing I learned about GIS this week (starting from knowing nothing about it last week) is that calculations are primarily done on an XY plane which are projections of a curved surface. I had naively assumed going in that calculations and indices (e.g. spatialite) were optimized for the surface of a sphere.
> >
>
> A small notice for a beginner: If you want to stay accurate and avoid lots of trouble with projections generated errors never leave the sphere (or ellipsoid) just keep doing all calculations there - it is much easier (you don't need projections) and you will stay accurate for ever ;)
>
> regards: Janne.
>
>
> _______________________________________________
> Proj mailing list
> [hidden email]
> http://lists.maptools.org/mailman/listinfo/proj
> _______________________________________________
> Proj mailing list
> [hidden email]
> http://lists.maptools.org/mailman/listinfo/proj
>

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

support.mn
In reply to this post by Demitri Muna
"Fischer, Robert P. (GISS-6110)[COLUMBIA UNIVERSITY]" [[hidden email]] kirjoitti:
> Although the theory for all this is probably sitting in a geometry book somewhere..

I find this web site very informative for a beginner:

http://aa.quae.nl/en/index.html

this shows all the most basic calculations with examples

regards: Janne.



_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

support.mn
In reply to this post by Demitri Muna
> "Fischer, Robert P. (GISS-6110)[COLUMBIA UNIVERSITY]" [[hidden email]] kirjoitti:
> > Although the theory for all this is probably sitting in a geometry book somewhere..
>
> I find this web site very informative for a beginner:
>
> http://aa.quae.nl/en/index.html
>
> this shows all the most basic calculations with examples
>
> regards: Janne.
>

still better link might be

http://aa.quae.nl/en/reken/grootcirkel.html#1

regards: Janne.

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

Fischer, Robert P. (GISS-6110)[COLUMBIA UNIVERSITY]

Although the theory for all this is probably sitting in a geometry book somewhere..
I was referring to this book:
http://www.amazon.com/Computational-Geometry-Surfaces-Performing-Cylinder/dp/1402002025

Computational Geometry on Surfaces: Performing Computational Geometry on the Cylinder, the Sphere, the Torus, and the Cone [Hardcover]

Clara I. Grima (Author), Alberto Márquez (Author)


I find this web site very informative for a beginner:

http://aa.quae.nl/en/index.html
http://aa.quae.nl/en/reken/grootcirkel.html#1

this shows all the most basic calculations with examples

Thanks, a lot of nice formulas there.  Unfortunately, it is not entirely accurate.  For example, Section 7, "Great Circle on a Topographic Map".  The Gnomonic Projection is a counter-example to the statement: " It is impossible to make a map of the world on which all great circles run straight."

if you need any accuracy and robustness you'll never leave the
sphere (or ellipsoid). It is much easier to spend a little more time
with the general spherical calculation than to fix it afterwards for
150+ projections for each one with maybe a different method.

This is not always the case.  Let me be specific.

Suppose you have a closed curve on a plane, and you want to know the area of it.  Suppose you have described the curve through a parameterization, that is functions x(z) and y(z), where z ranges from 0 to 1.  Then you can use Green's Theorem (Surveyor's Formula) to integrate and find the area of this curve.

Now suppose you want to apply this to the sphere.  That is, you have a curve on the surface of the sphere described by theata(z) and phi(z).  All you need to do is choose any old equal area projection.  The projection is described as x(theta,phi) and y(theta,phi).  Combine the two functions to get a curve on the plane, parameterized by z: x(theta(z), phi(z)), y(theta(z), phi(z)).  Now you can use Green's Theorem to integrate and find the area of this curve on the plane.  And since you chose an area-preserving projection, your answer will be exact.

This is essentially a derivation of Green's Theorem on a curved surface.

Another example where planar geometry can yield exact answers is the Cubed Sphere grid.  It uses the Gnomonic projection to create a regular grid made up of great circle paths.

Another trick is to project it to a single
standard projection which works all over the world and then
do the plane calculations if they are too hard to be done on the
sphere.
Unfortunately, there is no such thing.  Every projection has areas of the sphere where it fails miserably.  Not just places of high distortion (eg, Mercator near the poles), but also places where points that are close to each other on the globe end up far apart on the projection (eg, Mercator at the International Date Line).

-- Bob


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

Demitri Muna
In reply to this post by Fischer, Robert P. (GISS-6110)[COLUMBIA UNIVERSITY]
Hi,

Just to tie up a few loose ends, comments below. I've decided that GIS code, for the moment, won't give me significantly more functionality for what I need than other options that don't require the conceptual overhead, but thanks to everyone for all the comments and help. It's nice to learn about something new. :)

On Apr 16, 2012, at 6:24 PM, Fischer, Robert P. (GISS-6110)[COLUMBIA UNIVERSITY] wrote:

> Have you tried creating temporary schemas/tables in PostgreSQL with q3c?  Is that fast enough?  It sounds like that would be by far the easiest solution for you.

q3c is very fast; I haven't needed to even consider this option. My initial aim was to come up with an alternate tool that didn't require the overhead of a PostgreSQL installation for smaller data sets. But yes, for larger data sets, it's the only tool I use.

>> Thanks very much for the extensive comments. The biggest thing I
>> learned about GIS this week (starting from knowing nothing about it
>> last week) is that calculations are primarily done on an XY plane
>> which are projections of a curved surface. I had naively assumed going
>> in that calculations and indices (e.g. spatialite) were optimized for
>> the surface of a sphere.
>
> Yes, that's what I surmised too.  The good news is, the technique works quite well for local computations, and you often choose whatever projection you like.  For example, area of a general closed curve on the surface of a sphere can be computed this way, using any old equal area projection centered near your polygon.

I see now how 2D projections are very beneficial for geographical concerns, but less so for astronomical. That said, I do see some potential uses for visualization, so I might be back at some point. :)

>> Thanks for the code. I've been using the law of cosines and wasn't
>> aware of the problem at small separations (though this seems to trade
>> that problem for one at large separations).
>
> Vincenty's formula solves that problem.

Very glad to have learned that!

>> I will also look at the other links you mentioned. Since all my source
>> points are in ra/dec (not 3D cartesian), there will be the computation
>> hit to convert them back and forth.
>
> Is this signficant for your problem?

For computing a correlation function with a data set of half a million point sources, yes. But for smaller data sets, no. (And as I mentioned above, I was looking for a tool for the latter.)

Cheers,
Demitri
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: proj4 string for perfect sphere

support.mn
In reply to this post by Demitri Muna
Bob Fischer [[hidden email]] kirjoitti:

>
> > Another trick is to project it to a single
> > standard projection which works all over the world and then
> > do the plane calculations if they are too hard to be done on the
> > sphere.
> Unfortunately, there is no such thing.  Every projection has areas of
> the sphere where it fails miserably.  Not just places of high distortion
> (eg, Mercator near the poles), but also places where points that are
> close to each other on the globe end up far apart on the projection (eg,
> Mercator at the International Date Line).
>

I meant more some projection which can be made to work all over the world
by first setting some initial parameters (for example transverse mercator
projection with some polar switches [axis swapping at poles]) to do some
detail work not to be forced to stay on the sphere (or ellipsoid). But this of
course limits your work to detail level only.

regards: Janne.

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj