Re: Geometry vs Geography (what to use)

On Tue, Apr 5, 2016 at 10:07 AM, Lee Hachadoorian <
Lee(dot)Hachadoorian+L(at)gmail(dot)com> wrote:

> Mike,
>
> ST_DWithin takes a parameter as distance in the units of the SRID. Since
> 4326 is decimal degrees (latitude and longitude), you are manually
> converting from degrees to miles. This conversion will only come close to
> working if you are covering a pretty small area and if you look up the
> conversion based on the latitude, and will *always* be hobbled that that a
> degree longitude is not the same length as a degree latitude (although they
> are close at the equator).
>
> The standard way to accomplish this is to either project the coordinates
> to an appropriate planar coordinate reference system, or to use the
> geography data type. If you use geography, the units are in meters.
>
> You could ask your DBAs to change the SRID of the table, but you can also
> transform them as part of your query.
>
> Option 1 (using SRID for UTM 10N, see below)
>
> SELECT pc1.city AS city,
> pc1.postalcode AS postalcode
> FROM tpostalcoordinate pc1
> WHERE ST_DWithin(ST_Transform(geo_position, 32610),
> (SELECT ST_Transform(pc2.geo_position, 32610)
> FROM tpostalcoordinate pc2
> WHERE pc2.postalcode = '95050'
> AND pc2.countrycode2tcountry = 'US'),
> 80467.2)
>
> I've selected SRID 32610 as the appropriate UTM zone based on the searched
> ZIP code 95050. If you need to do these calculations over a different area,
> you would have to use a different UTM zone, or pick a wide-area CRS
> (something that would apply to the continental US), with some loss of
> accuracy.
>
> Option 2 (using geography)
>
> SELECT pc1.city AS city,
> pc1.postalcode AS postalcode
> FROM tpostalcoordinate pc1
> WHERE ST_DWithin(geo_position::geography),
> (SELECT pc2.geo_position::geography
> FROM tpostalcoordinate pc2
> WHERE pc2.postalcode = '95050'
> AND pc2.countrycode2tcountry = 'US'),
> 80467.2)
>
> Neither query has been tested.
>
> Best,
> --Lee
>
> On Tue, Apr 5, 2016 at 11:28 AM, Michael Moore <michaeljmoore(at)gmail(dot)com>
> wrote:
>
>>
>>
>> On Mon, Apr 4, 2016 at 6:09 PM, Steve Midgley <science(at)misuse(dot)org> wrote:
>>
>>> On Mon, Apr 4, 2016 at 5:20 PM, Michael Moore <michaeljmoore(at)gmail(dot)com>
>>> wrote:
>>>
>>>> I am converting from Oracle to postgres. We have an application that
>>>> takes a value for MILES as an input and returns a list of SCHOOLS that are
>>>> within range. This all works fine in Oracle, no problem. In Oracle the
>>>> datatype of the table field is SDO_GEOMETRY. Our DBAs have set up the same
>>>> table in Postgres with datatype geometry(Point,4326). Here is where the
>>>> problem comes in. I am trying to use the ST_DWithin function and it kind of
>>>> works except that one of the input parameters to this function is expressed
>>>> in RADIANS.
>>>>
>>>> SELECT pc1.city AS city,
>>>> pc1.postalcode AS postalcode
>>>> FROM tpostalcoordinate pc1
>>>> WHERE ST_DWithin(geo_position,
>>>> (SELECT pc2.geo_position
>>>> FROM tpostalcoordinate pc2
>>>> WHERE pc2.postalcode = '95050'
>>>> AND pc2.countrycode2tcountry = 'US'),
>>>> 50 * .01539)
>>>> So, in the above query I am selecting all the cities with a 50 mile
>>>> radius. It WORKS ... kind of. The thing I don't like is my screwy way of
>>>> converting radians to miles. ( miles * .01539 ). It's just a best guess
>>>> that gives me results 'almost' the same as Oracle's. The second thing I
>>>> don't like is the 'almost' the same. I'm am guessing the difference in the
>>>> result set is due to the planar calculation vs the spheroid
>>>> calculations. So my questions are
>>>>
>>>> 1. Should I be using GEOGRAPHY(POINT,4326) instead of GEOMETRY(POINT,4326)I
>>>> assume this would allow me to express the distance meters and it would do spheroid
>>>> calculations which should give me results more consistent with Oracle's?
>>>> 2. Any advice. Is there something else I should be doing? What did
>>>> I miss?
>>>>
>>>> Just reading this article
>>> <http://gis.stackexchange.com/questions/26082/what-is-the-difference-between-geometric-and-geographic-columns>[1] makes
>>> me think that you might get better spherical calculations for lat/long
>>> using the Geography data type.
>>>
>>> Regarding ST_DWithin itself, I'm confused about your convertion to
>>> radians. The docs I'm looking at say that the third parameter to that
>>> function is meters, not radians. (And radians aren't a measure of distance
>>> anyway, as I'm sure you know, but which is also confusing me in your
>>> message). So are you just converting 50 miles to meters?
>>>
>>> Also, maybe your conversion factor is accurate, and it is possible that
>>> the *Oracle side* is using a linear projection (treating lat/long like a
>>> cartesian grid), and you are doing proper spherical math on the Pg side,
>>> and that's why your results are close but not quite right?
>>>
>>> I hope this helps your work.. Best,
>>> Steve
>>>
>>> [1]
>>> http://gis.stackexchange.com/questions/26082/what-is-the-difference-between-geometric-and-geographic-columns
>>>
>>> Steve,
>> I think I read somewhere in the documents that Radians was the parameter
>> for Geometry and Meters is the parameter for Geography. Even with Radians,
>> if the radius of the circle is known then so too is it for 1 radian. I'm
>> guessing there is a default circle size for a geometry(Point,4326)
>> datatype. I'll keep digging and look at the docs you suggested.
>> Regards,
>> Mike
>>
>>
>
>
> --
> Lee Hachadoorian
> Asst Professor of Geography, Dartmouth College
> http://freecity.commons.gc.cuny.edu/
>

Lee,
I tried casting to geography, but I get this:
ERROR: GetProj4StringSPI: Cannot find SRID (8307) in spatial_ref_sys
********** Error **********
So, I discovered that "select * from spatial_ref_sys;" gives no results,
meaning that the table is empty. I'll be talking with our DBAs about this.

That being as it may, I read on somebody's blog that casting to geography
can really slow things down so my plan is to add a new column like this:
alter table tpostalcoordinate add column geography_position
geography(POINT,4326) ;
then I will populate it like this:
UPDATE tpostalcoordinate set geography_position = ST_SetSRID(ST_Point(
longitude, latitude), 4326);
and build an index like:
CREATE INDEX tpostal_geo_geography_idx ON tpostalcoordinate USING
gist(geography_position);

I'll let every know how it goes.