# [postgis-users] Spheroid Question

Gerry Creager gerry.creager at tamu.edu
Fri Jun 24 17:19:17 PDT 2005

```I shouldn't think without checking references....

of GRS80 parameters.  Equatorial radius is really, 6378137 m, while the

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

```