[postgis-users] Spheroid Question

[Stephen Woodbridge]

Gerry,

Excellent point. And 6378137 m * 2 * pi = 40075016 and I mis-spoke below 
  as I meant to say circumference and I mis-typed it as 40075160, but 
the  111319 m per degree is still a good number for approximating. 
Remember you are setting this for high level rejection and the accurate 
math is be done with the *_spheroid() functions.

The key is to not use any of these as magic formula or numbers or you 
will wind up using them in a wrong way eventually.

-Steve

Gerry Creager wrote:
> I shouldn't think without checking references....
> 
> See http://dgfi2.dgfi.badw-muenchen.de/geodis/REFS/grs80.html for a list 
> of GRS80 parameters.  Equatorial radius is really, 6378137 m, while the 
> polar radius is 6356752.3141 m.
> 
> gerry
> 
> Stephen Woodbridge wrote:
> 
>> Dan,
>>
>> I understood that you probably wanted '2000' to mean meters, but just 
>> because that is your intent, does not make the code nessacarily 
>> execute that way. If the code interpreted 2000 as decimal degrees 
>> instead of 2000 meters your SQL would still work and return the same 
>> results. In fact your SQL will work and return the same results if you 
>> remove
>>
>> the_geom && expand(setSRID(makepoint(-72.607912, 
>> 42.058052),4326),2000)  AND
>>
>> from the SQL statement. So why bother adding the phrase. The phrase is 
>> supposed to reject objects that are outside the expanded box from 
>> consideration from the clause that follows which is more expensive.
>>
>> Since your point is in SRID:4326 which is in decimal degrees, when you 
>> expand the point's bbox you have to expand it in the units of that 
>> projection, in this case, decimal degrees. So your example code is 
>> expanding it by +-2000 decimal degrees (not meters as you expected). 
>> It will reject NO objects and ALL the object have to go through the 
>> expensive test, AND you get the correct result, AND it costs a lot of 
>> time and required you to change the distance calc to be less costly 
>> and less accurate one.
>>
>> So the radius of the earth is about 40075160 meters (from my high 
>> school geography book) and there are 360 degrees in a ball or 
>> approximately 111319 meters per degree or 0.179664 degrees in 2000 
>> meters.
>>
>> So if you go back to your original SQL and change it to say:
>>
>> the_geom && expand(setSRID(makepoint(-72.607912, 
>> 42.058052),4326),0.179664)  AND
>>
>> And rerun it, it should be MUCH faster and give you the same results.
>>
>> I hope this all makes sense, and I hope that someone will speak up if 
>> I have this wrong, but I think it is correct.
>>
>> Best regards,
>>   -Stephen Woodbridge
>>
>> Dan Phillips wrote:
>>
>>> Hey Steve,
>>>
>>> The '2000' represents 2000 meters & the query returns the correct 
>>> results
>>> for both distance_spheroid() & distance_sphere(). As previously noted,
>>> distance_sphere() is much faster. The units returned are in meters.
>>>
>>> Dan
>>>
>>>
>>>
>>> -----Original Message-----
>>> From: postgis-users-bounces at postgis.refractions.net
>>> [mailto:postgis-users-bounces at postgis.refractions.net]On Behalf Of
>>> Stephen Woodbridge
>>> Sent: Friday, June 24, 2005 4:18 PM
>>> To: PostGIS Users Discussion
>>> Subject: Re: [postgis-users] Spheroid Question
>>>
>>>
>>> Dan,
>>>
>>> I am be wrong on my understanding, but it looks like you might be doing
>>> a full table scan. Since SRID:4326 is in decimal degrees, and you are
>>> expanding you point by +-2000 decimal degrees which is likely
>>> overlapping every object in the table. If you want to select items in
>>> 2km of your point you might want to convert 2000 into its approximate
>>> degrees for the && clause.
>>>
>>> Also what units do the *_spheroid() return in?
>>>
>>> -Steve W.
>>>
>>> Dan Phillips wrote:
>>>
>>>> Okay, this works. It's just really slow (~40 minutes). Any 
>>>> suggestions on
>>>> speeding this up?
>>>>
>>>> SELECT
>>>>  CRID,
>>>>  distance_spheroid(centroid(the_geom),
>>>>  setSRID(makepoint(-72.607912, 42.058052),4326),
>>>>  'SPHEROID["GRS_1980",6378137,298.257222101]') as distance
>>>> FROM crrts
>>>> WHERE
>>>>  the_geom && expand(setSRID(makepoint(-72.607912, 
>>>> 42.058052),4326),2000)
>>>>  AND distance_spheroid(centroid(the_geom),
>>>>  setSRID(makepoint(-72.607912, 42.058052),4326),
>>>>  'SPHEROID["GRS_1980",6378137,298.257222101]')<2000 ;
>>>>
>>>> Thanks,
>>>>
>>>> Dan
>>>>
>>>> -----Original Message-----
>>>> From: postgis-users-bounces at postgis.refractions.net
>>>> [mailto:postgis-users-bounces at postgis.refractions.net]On Behalf Of
>>>> strk at refractions.net
>>>> Sent: Friday, June 24, 2005 11:39 AM
>>>> To: PostGIS Users Discussion
>>>> Subject: Re: [postgis-users] Spheroid Question
>>>>
>>>>
>>>> On Fri, Jun 24, 2005 at 11:36:27AM -0400, Dan Phillips wrote:
>>>>
>>>>
>>>>> Thanks Jeff & Paul. That was the problem.
>>>>>
>>>>> Now I'm getting "ERROR:  Operation on two geometries with different 
>>>>> SRIDs"
>>>>>
>>>>> I assume this is because I'm not setting the SRID for my 
>>>>> makepoint(). How
>>>>
>>>>
>>>>
>>>> do
>>>>
>>>>
>>>>> I do this to make it match the SRID (4326) of my 'crrts' table?
>>>>
>>>>
>>>>
>>>>
>>>> Use the setSRID() function:
>>>> setSRID(makepoint())
>>>> --strk;
>>
>>
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at postgis.refractions.net
>> http://postgis.refractions.net/mailman/listinfo/postgis-users
> 
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>