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

[Michael Moore]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20160411/af2fb55f/attachment.html>