[postgis-users] Spheroid Question

[Paul Ramsey]

You are correct Stephen. I did not read the SQL particularly closely, 
but it is indeed exanding 2000 units (which are degrees in this case). 
And your analysis of the fix is also right.

Paul

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