[postgis-users] Spheroid Question

[Stephen Woodbridge]

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;