[postgis-users] Using SRID with meters, but still not getting UNITS in meters with ST_DWithin

[Michael Moore]

On Mon, Apr 11, 2016 at 12:31 PM, Michael Moore <michaeljmoore at gmail.com>
wrote:

>
>
> On Mon, Apr 11, 2016 at 11:50 AM, Paul Ramsey <pramsey at cleverelephant.ca>
> wrote:
>
>> Before I can even begin to address what's going on with the
>> index/query side, I have to ask about the coordinate system:
>>
>> You are using a geocentric system for postal code points?
>>
>> https://epsg.io/4896
>>
>> This seems very odd indeed. If you were flying a satellite or drilling
>> a deep well, maybe you'd use geocentric coordinates, but simple postal
>> code queries? So, perhaps there's a core error to deal with first: why
>> are you in EPSG:4896?
>>
>> P
>>
>>
>> On Mon, Apr 11, 2016 at 11:21 AM, Michael Moore <michaeljmoore at gmail.com>
>> wrote:
>> > I am trying to find all zip codes withing a given range of current zip
>> code.
>> > For example, User is at zip code 95076 and want to know all zip codes
>> within
>> > a 30 mile range. This will always be a short distance, nothing over 50
>> > miles. The source zip code is any place in the USA.  My input table has
>> a
>> > latitude field, a longitude field and a POINT field for each zip code.
>> > I want to supply in input in meters, not degrees. I do NOT want to cast
>> to
>> > geography because this prevents the index from being used.
>> > Here is a working query that uses the index, but ST_DWithin  is not
>> > understanding ".05" as meters:
>> >
>> >
>> > lcd1_dev=>   select t2.city, t2.postalcode,
>> > ST_Distance(t1.geo_position,t2.geo_position)
>> > lcd1_dev->   from      tpostalcoordinate t1
>> > lcd1_dev->   left join tpostalcoordinate t2 on
>> > ST_DWithin(t1.geo_position,t2.geo_position , .05)
>> > lcd1_dev->   where t1.postalcode = '94404'and t1.countrycode2tcountry =
>> 'US'
>> > and t2.countrycode2tcountry= 'US';
>> >      city     | postalcode |     st_distance
>> > --------------+------------+---------------------
>> >  Redwood City | 94065      |  0.0273766323714193
>> >  San Mateo    | 94408      | 0.00504738546179631
>> >  Belmont      | 94002      |  0.0440065904155286
>> >  San Mateo    | 94404      |                   0
>> >  San Mateo    | 94403      |  0.0370314731005901
>> >  San Mateo    | 94407      |  0.0416118372581607
>> >
>> > This shows that I am using SRID 4896
>> > lcd1_dev=> select * from geometry_columns where f_table_name =
>> > 'tpostalcoordinate';
>> >  f_table_catalog | f_table_schema |   f_table_name    |
>> f_geometry_column |
>> > coord_dimension | srid | type
>> >
>> -----------------+----------------+-------------------+-------------------+-----------------+------+-------
>> >  lcd1_dev        | qsn_app        | tpostalcoordinate | geo_position
>>   |
>> > 2 | 4896 | POINT
>> >
>> > 4896 is UNIT "metre" as show here:
>> >  4896 | EPSG      |      4896 |
>> >
>> GEOCCS["ITRF2005",DATUM["International_Terrestrial_Reference_Frame_2005",SPHEROID["GRS
>> > 1980",6378137,298.257222101,AUTHORIT
>> >
>> Y["EPSG","7019"]],AUTHORITY["EPSG","6896"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Geocentric
>> > X",OTH
>> > ER],AXIS["Geocentric Y",OTHER],AXIS["Geocentric
>> > Z",NORTH],AUTHORITY["EPSG","4896"]] | +proj=geocent +ellps=GRS80
>> +units=m
>> > +no_de
>> >
>> > The following (from EXPLAIN) proves that the index is being used in the
>> > working query shown above:
>> >    ->  Index Scan using tpostalcoordinate_pk on tpostalcoordinate t1
>> > (cost=0.42..8.45 rows=1 width=32)
>> >
>> > What do I need to do to get this in meters without losing my index
>> access?
>> >
>> > TIA,
>> > Mike
>> >
>> > _______________________________________________
>> > postgis-users mailing list
>> > postgis-users at lists.osgeo.org
>> > http://lists.osgeo.org/mailman/listinfo/postgis-users
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/postgis-users
>
>
> Here is what we do to create the column on the table:
>
>
>
> 1.    SELECT AddGeometryColumn('qsn_app','tpostalcoordinate',
> 'geo_position',4326,'POINT',2);
>
>
>
> 2.    UPDATE tpostalcoordinate pc set pc.geo_position =
> ST_SetSRID(ST_Point(pc.longitude, pc.latitude), 4326);
>
>
>
> 3.    CREATE INDEX tpostal_geo_position_idx ON tpostalcoordinate USING
> gist(geo_position);
>
>
> I am on EPSG:4896 because I tried about 10 others that were also meters
> and that's just the one I ended up on .  There is one thing you said that
> leads me to believe there is something I don't understand: "...geocentric
> system...". My understanding is that there are two datatypes to choose
> from ; geometry and geography. I have chosen geometry because 1. it should
> be faster due to simpler calculations and 2. It should be faster because it
> can use the index.  Beyond that I thought that I would then need to choose
> a proper SRID. The two criteria for choosing an SRID for me are 1. the
> units are meters and 2. it is based on a datum that is appropriate for the
> USA. Now the fact that you said: "geocentric system" makes me think that
> maybe there is something else I need to be looking at when choosing a
> proper SRID.  A little more on the application: The user wants to know all
> of the schools that are within a specific distance of his current location,
> his zip code. Precision is not important. We get the user's zip code from
> the screen and we have a list of schools and their zip codes.
>
>
> Thanks,
>
> Mike
>

John,
The latitude and longitude are correct and I know that because if I cast my
geometry points to geography, every thing works perfectly: it understand
meters and expected target zip codes are in the result set. The only
problem with with casting to geography is that geography does not use the
index.
Regards,
Mike
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20160411/a0fa51b7/attachment.html>