<div dir="ltr">Okay, let's see if I finally get this. The reason I need to start with SRID 4326 is because the Coordinate system has an Axes of latitude and longitude and I only have latitude and longitude to start with. My starting latitude and longitude were based on the WSG84 Datum. looking at Coordinate system for <a href="https://epsg.io/4326">https://epsg.io/4326</a>. 4326 is also based on WSG84, so it's a match. Unfortunately for me, 4326 is ellipsoidal and as such<div>the Unit is degrees. So in order do do my calculations in meters, I need to "translate" my point onto a flat map. 2163 represents such a map as as would be expected with a flat map, the unit is in meters. </div><div>Also there is the little issue of the "covered area" that needs to be taken into consideration but I understand that. <br>Thanks!</div><div>Mike </div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Apr 12, 2016 at 12:50 PM, Paul Ramsey <span dir="ltr"><<a href="mailto:pramsey@cleverelephant.ca" target="_blank">pramsey@cleverelephant.ca</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">This:<br>
<span class=""><br>
"while lat-long SRIDs (like 4326) say the units are "meters""<br>
<br>
</span>is not correct. Long-lat SRS's pretty clearly state their units are<br>
degrees. More to the point, the units are evaluated in a spherical<br>
context, not a planar one, so applying planar calculations to them<br>
will yield junk answers. You want to apply planar calculations to<br>
planar coordinates, and spherical calculations to spherical<br>
coordinates.<br>
<br>
This is why we have the geometry (planar) and geography (spherical) bifurcation.<br>
<br>
In my previous email I gave you a recipe to run your problem in<br>
geography (spherical) space. Since you've sort of indicated that your<br>
data are USA only, I'll also give you a planar (geometry) recipe,<br>
using a suitable planar projection for the USA.<br>
<br>
We will use the USA National Atlas Albers projection<br>
(<a href="https://epsg.io/2163" rel="noreferrer" target="_blank">https://epsg.io/2163</a>) for this example. It won't provide meter-level<br>
accuracy over such a wide area, but presumable for something as coarse<br>
as a "postal codes in radius" query, that won't matter.<br>
<br>
ALTER TABLE tpostalcoordinate ADD COLUMN geom geometry(point, 2136);<br>
UPDATE tpostalcoordinate SET geom =<br>
ST_Transform(ST_SetSRID(ST_MakePoint(longitude, latitude), 4326),2163);<br>
CREATE INDEX tpostalcoordinate_geom_ix ON tpostalcoordinate USING GIST (geom);<br>
<span class="">SELECT a.*<br>
FROM tpostalcoordinate a<br>
JOIN tpostalcoordinate b<br>
</span> ON ST_DWithin(a.geom, b.geom, 10000)<br>
<span class=""> WHERE<br>
a.countrycode2tcountry = 'US' AND<br>
b.countrycode2tcountry = 'US' AND<br>
b.postalcode = '94404';<br>
<br>
</span>Note the differences with the geography approach:<br>
<br>
- use use the geometry type, and we populate the column with data<br>
transformed into our working planar projection<br>
- then we're free to just use ST_DWithin(geom, geom, radius) and<br>
everything magically works fine.<br>
<br>
ATB,<br>
<br>
P<br>
<div class="HOEnZb"><div class="h5"><br>
<br>
On Tue, Apr 12, 2016 at 12:40 PM, Michael Moore <<a href="mailto:michaeljmoore@gmail.com">michaeljmoore@gmail.com</a>> wrote:<br>
> Lee, Paul,<br>
> Thanks so much for your time.<br>
> Lee,<br>
> My understanding is that there are about 42,000 zips in the US, So a cross<br>
> join should result in about 1.7 Billion entries into the<br>
> postal_zone_distance table. Maybe that's not so bad, I'll have to run the<br>
> numbers. Where as this might be the fastest possible solution, it's probably<br>
> over-kill for my meager needs. I understand what you mean by needing to know<br>
> first the SRID of the Latitude and Longitude from my source table. It means<br>
> very little without knowing the Datum from which they were derived. Thanks<br>
> to some youtube videos I think I grasp this.<br>
><br>
> You said:<br>
>><br>
>> ...while lat-long SRIDs (like 4326) say the units are "meters",<br>
>> ST_Distance and other distance and area functions always just work in<br>
>> Cartesian coordinates, meaning they treat long and lat as X and Y in a<br>
>> Cartesian plan and they take parameters and return answers in degrees, which<br>
>> isn't really meaningful.<br>
><br>
><br>
> I believe you and my experience seems to confirm it, however I can't seem<br>
> to reconcile this with the documentation for ST_Distance which says ...<br>
><br>
> <a href="http://postgis.net/docs/ST_DWithin.html" rel="noreferrer" target="_blank">http://postgis.net/docs/ST_DWithin.html</a><br>
>><br>
>> For Geometries: The distance is specified in units defined by the spatial<br>
>> reference system of the geometries. For this function to make sense, the<br>
>> source geometries must both be of the same coordinate projection, having the<br>
>> same SRID.<br>
><br>
><br>
> Is the documentation wrong??? or am I just not interpreting it correctly.<br>
><br>
> If it IS wrong, for god sake I hope they fix it as I have now just wasted<br>
> too much of everybody's time on and issue that could have been cleared up by<br>
> better documentation.<br>
><br>
> Again, thanks to everybody who replied, you are most generous.<br>
> Regards,<br>
> Mike<br>
><br>
><br>
><br>
> On Tue, Apr 12, 2016 at 9:26 AM, Lee Hachadoorian<br>
> <<a href="mailto:Lee.Hachadoorian%2BL@gmail.com">Lee.Hachadoorian+L@gmail.com</a>> wrote:<br>
>><br>
>> Michael,<br>
>><br>
>> There are many ways to accomplish this, and the most efficient approach is<br>
>> going to depend on what else you are using this data for. If literally the<br>
>> *only* thing you are using the geometries for is to calculate distances<br>
>> between postal zones, I would consider creating a table like:<br>
>><br>
>> CREATE TABLE postal_zone_distance (<br>
>> source_zone text,<br>
>> dest_zone text,<br>
>> distance_m double precision,<br>
>> PRIMARY KEY (source_zone, dest_zone)<br>
>> );<br>
>><br>
>> and populating it *once* with distances between all postal codes in the<br>
>> US, so that you can run future queries using an attribute select (WHERE<br>
>> distance_m < 50000) instead of a spatial select. Rerun whenever postal zones<br>
>> change, which shouldn't be all that often. You can populate it with<br>
>> something like:<br>
>><br>
>> INSERT INTO postal_zone_distance<br>
>> SELECT a.postalcode, b.postalcode, ST_Distance(a.geo_position::geography,<br>
>> b.geo_position::geography)<br>
>> FROM tpostalcoordinate a CROSS JOIN tpostalcoordinate b;<br>
>><br>
>> Will take a while to run, but future queries will be faster.<br>
>><br>
>> If you intend to use this for more than just these distances, there are<br>
>> several options, all of which are valid, but again it depends what you want<br>
>> to use it for. The main issue between using Paul's suggestion to just define<br>
>> columns as geography vs keeping in geometry is whether you might visualize<br>
>> the postal zones using a GIS. If you will be using QGIS or ArcGIS to connect<br>
>> to PostGIS to map the zones, you will need to keep in geometry. In this case<br>
>> you can, as Lucas suggest, define a functional index on the geography cast,<br>
>> while keeping the underlying data in geometry.<br>
>><br>
>> Some other notes:<br>
>><br>
>> * To get to your original question about units of measurement, while<br>
>> lat-long SRIDs (like 4326) say the units are "meters", ST_Distance and other<br>
>> distance and area functions always just work in Cartesian coordinates,<br>
>> meaning they treat long and lat as X and Y in a Cartesian plan and they take<br>
>> parameters and return answers in degrees, which isn't really meaningful.<br>
>> There are special functions like ST_Distance_Spheroid which take lat-long<br>
>> geometry and return distance in meters on the spheroid. And ST_Distance on<br>
>> geography type also returns meters.<br>
>><br>
>> * Be careful about your choice of SRID. The fact that you tried to<br>
>> transform (or assign?) this to EPSG 4896 is a red flag because because one<br>
>> shouldn't just "try" SRIDs. Your data come to you in an SRID, and you should<br>
>> only transform it for a particular purpose. In fact you should make sure<br>
>> that you know the datum of your original lat-long geometries. Considering<br>
>> that the data are for the US, it is quite possible that your data are in<br>
>> EPSG 4269 (lat-long North American Datum 1983, used for much US government<br>
>> data including by the Census Bureau), rather than being EPSG 4326 (World<br>
>> Geodetic System 1984, the datum used by GPS satellites).<br>
>><br>
>> I highly recommend Obe & Hsu's PostGIS in Action<br>
>> (<a href="https://www.manning.com/books/postgis-in-action-second-edition" rel="noreferrer" target="_blank">https://www.manning.com/books/postgis-in-action-second-edition</a>), esp. Ch 3<br>
>> on Spatial Reference System Considerations.<br>
>><br>
>> Best,<br>
>> --Lee<br>
>><br>
>> --<br>
>> Lee Hachadoorian<br>
>> Assistant Professor of Instruction, Geography & Urban Studies<br>
>> Assistant Director, Professional Science Master's in GIS<br>
>> Temple University<br>
>><br>
>> On Mon, Apr 11, 2016 at 11:00 PM, Paul Ramsey <<a href="mailto:pramsey@cleverelephant.ca">pramsey@cleverelephant.ca</a>><br>
>> wrote:<br>
>>><br>
>>> Go back to basics, use a geography column from the start:<br>
>>><br>
>>> ALTER TABLE tpostalcoordinate ADD COLUMN geog geography(point, 4326);<br>
>>> UPDATE tpostalcoordinate SET geog =<br>
>>> geography(ST_SetSRID(ST_MakePoint(longitude, latitude), 4326));<br>
>>> CREATE INDEX tpostalcoordinate_gix ON tpostalcoordinate USING GIST<br>
>>> (geog);<br>
>>> SELECT a.*<br>
>>> FROM tpostalcoordinate a<br>
>>> JOIN tpostalcoordinate b<br>
>>> ON ST_DWithin(a.geog, b.geog, 10000)<br>
>>> WHERE<br>
>>> a.countrycode2tcountry = 'US' AND<br>
>>> b.countrycode2tcountry = 'US' AND<br>
>>> b.postalcode = '94404';<br>
>>><br>
>>> A geocentric system is one where every point is referenced relative to<br>
>>> the center of the earth, so they all have an x,y,z. It's not something<br>
>>> you want to use.<br>
>>><br>
>>> P<br>
>>><br>
>>><br>
>>> On Mon, Apr 11, 2016 at 12:31 PM, Michael Moore <<a href="mailto:michaeljmoore@gmail.com">michaeljmoore@gmail.com</a>><br>
>>> wrote:<br>
>>> ><br>
>>> ><br>
>>> > On Mon, Apr 11, 2016 at 11:50 AM, Paul Ramsey<br>
>>> > <<a href="mailto:pramsey@cleverelephant.ca">pramsey@cleverelephant.ca</a>><br>
>>> > wrote:<br>
>>> >><br>
>>> >> Before I can even begin to address what's going on with the<br>
>>> >> index/query side, I have to ask about the coordinate system:<br>
>>> >><br>
>>> >> You are using a geocentric system for postal code points?<br>
>>> >><br>
>>> >> <a href="https://epsg.io/4896" rel="noreferrer" target="_blank">https://epsg.io/4896</a><br>
>>> >><br>
>>> >> This seems very odd indeed. If you were flying a satellite or drilling<br>
>>> >> a deep well, maybe you'd use geocentric coordinates, but simple postal<br>
>>> >> code queries? So, perhaps there's a core error to deal with first: why<br>
>>> >> are you in EPSG:4896?<br>
>>> >><br>
>>> >> P<br>
>>> >><br>
>>> >><br>
>>> >> On Mon, Apr 11, 2016 at 11:21 AM, Michael Moore<br>
>>> >> <<a href="mailto:michaeljmoore@gmail.com">michaeljmoore@gmail.com</a>><br>
>>> >> wrote:<br>
>>> >> > I am trying to find all zip codes withing a given range of current<br>
>>> >> > zip<br>
>>> >> > code.<br>
>>> >> > For example, User is at zip code 95076 and want to know all zip<br>
>>> >> > codes<br>
>>> >> > within<br>
>>> >> > a 30 mile range. This will always be a short distance, nothing over<br>
>>> >> > 50<br>
>>> >> > miles. The source zip code is any place in the USA. My input table<br>
>>> >> > has<br>
>>> >> > a<br>
>>> >> > latitude field, a longitude field and a POINT field for each zip<br>
>>> >> > code.<br>
>>> >> > I want to supply in input in meters, not degrees. I do NOT want to<br>
>>> >> > cast<br>
>>> >> > to<br>
>>> >> > geography because this prevents the index from being used.<br>
>>> >> > Here is a working query that uses the index, but ST_DWithin is not<br>
>>> >> > understanding ".05" as meters:<br>
>>> >> ><br>
>>> >> ><br>
>>> >> > lcd1_dev=> select t2.city, t2.postalcode,<br>
>>> >> > ST_Distance(t1.geo_position,t2.geo_position)<br>
>>> >> > lcd1_dev-> from tpostalcoordinate t1<br>
>>> >> > lcd1_dev-> left join tpostalcoordinate t2 on<br>
>>> >> > ST_DWithin(t1.geo_position,t2.geo_position , .05)<br>
>>> >> > lcd1_dev-> where t1.postalcode = '94404'and<br>
>>> >> > t1.countrycode2tcountry =<br>
>>> >> > 'US'<br>
>>> >> > and t2.countrycode2tcountry= 'US';<br>
>>> >> > city | postalcode | st_distance<br>
>>> >> > --------------+------------+---------------------<br>
>>> >> > Redwood City | 94065 | 0.0273766323714193<br>
>>> >> > San Mateo | 94408 | 0.00504738546179631<br>
>>> >> > Belmont | 94002 | 0.0440065904155286<br>
>>> >> > San Mateo | 94404 | 0<br>
>>> >> > San Mateo | 94403 | 0.0370314731005901<br>
>>> >> > San Mateo | 94407 | 0.0416118372581607<br>
>>> >> ><br>
>>> >> > This shows that I am using SRID 4896<br>
>>> >> > lcd1_dev=> select * from geometry_columns where f_table_name =<br>
>>> >> > 'tpostalcoordinate';<br>
>>> >> > f_table_catalog | f_table_schema | f_table_name |<br>
>>> >> > f_geometry_column |<br>
>>> >> > coord_dimension | srid | type<br>
>>> >> ><br>
>>> >> ><br>
>>> >> > -----------------+----------------+-------------------+-------------------+-----------------+------+-------<br>
>>> >> > lcd1_dev | qsn_app | tpostalcoordinate | geo_position<br>
>>> >> > |<br>
>>> >> > 2 | 4896 | POINT<br>
>>> >> ><br>
>>> >> > 4896 is UNIT "metre" as show here:<br>
>>> >> > 4896 | EPSG | 4896 |<br>
>>> >> ><br>
>>> >> ><br>
>>> >> > GEOCCS["ITRF2005",DATUM["International_Terrestrial_Reference_Frame_2005",SPHEROID["GRS<br>
>>> >> > 1980",6378137,298.257222101,AUTHORIT<br>
>>> >> ><br>
>>> >> ><br>
>>> >> > Y["EPSG","7019"]],AUTHORITY["EPSG","6896"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Geocentric<br>
>>> >> > X",OTH<br>
>>> >> > ER],AXIS["Geocentric Y",OTHER],AXIS["Geocentric<br>
>>> >> > Z",NORTH],AUTHORITY["EPSG","4896"]] | +proj=geocent +ellps=GRS80<br>
>>> >> > +units=m<br>
>>> >> > +no_de<br>
>>> >> ><br>
>>> >> > The following (from EXPLAIN) proves that the index is being used in<br>
>>> >> > the<br>
>>> >> > working query shown above:<br>
>>> >> > -> Index Scan using tpostalcoordinate_pk on tpostalcoordinate t1<br>
>>> >> > (cost=0.42..8.45 rows=1 width=32)<br>
>>> >> ><br>
>>> >> > What do I need to do to get this in meters without losing my index<br>
>>> >> > access?<br>
>>> >> ><br>
>>> >> > TIA,<br>
>>> >> > Mike<br>
>>> >> ><br>
>>> >> > _______________________________________________<br>
>>> >> > postgis-users mailing list<br>
>>> >> > <a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br>
>>> >> > <a href="http://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" target="_blank">http://lists.osgeo.org/mailman/listinfo/postgis-users</a><br>
>>> >> _______________________________________________<br>
>>> >> postgis-users mailing list<br>
>>> >> <a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br>
>>> >> <a href="http://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" target="_blank">http://lists.osgeo.org/mailman/listinfo/postgis-users</a><br>
>>> ><br>
>>> ><br>
>>> > Here is what we do to create the column on the table:<br>
>>> ><br>
>>> ><br>
>>> ><br>
>>> > 1. SELECT AddGeometryColumn('qsn_app','tpostalcoordinate',<br>
>>> > 'geo_position',4326,'POINT',2);<br>
>>> ><br>
>>> ><br>
>>> ><br>
>>> > 2. UPDATE tpostalcoordinate pc set pc.geo_position =<br>
>>> > ST_SetSRID(ST_Point(pc.longitude, pc.latitude), 4326);<br>
>>> ><br>
>>> ><br>
>>> ><br>
>>> > 3. CREATE INDEX tpostal_geo_position_idx ON tpostalcoordinate USING<br>
>>> > gist(geo_position);<br>
>>> ><br>
>>> ><br>
>>> > I am on EPSG:4896 because I tried about 10 others that were also meters<br>
>>> > and<br>
>>> > that's just the one I ended up on . There is one thing you said that<br>
>>> > leads<br>
>>> > me to believe there is something I don't understand: "...geocentric<br>
>>> > system...". My understanding is that there are two datatypes to choose<br>
>>> > from<br>
>>> > ; geometry and geography. I have chosen geometry because 1. it should<br>
>>> > be<br>
>>> > faster due to simpler calculations and 2. It should be faster because<br>
>>> > it can<br>
>>> > use the index. Beyond that I thought that I would then need to choose<br>
>>> > a<br>
>>> > proper SRID. The two criteria for choosing an SRID for me are 1. the<br>
>>> > units<br>
>>> > are meters and 2. it is based on a datum that is appropriate for the<br>
>>> > USA.<br>
>>> > Now the fact that you said: "geocentric system" makes me think that<br>
>>> > maybe<br>
>>> > there is something else I need to be looking at when choosing a proper<br>
>>> > SRID.<br>
>>> > A little more on the application: The user wants to know all of the<br>
>>> > schools<br>
>>> > that are within a specific distance of his current location, his zip<br>
>>> > code.<br>
>>> > Precision is not important. We get the user's zip code from the screen<br>
>>> > and<br>
>>> > we have a list of schools and their zip codes.<br>
>>> ><br>
>>> ><br>
>>> > Thanks,<br>
>>> ><br>
>>> > Mike<br>
>>> ><br>
>>> ><br>
>>> > _______________________________________________<br>
>>> > postgis-users mailing list<br>
>>> > <a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br>
>>> > <a href="http://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" target="_blank">http://lists.osgeo.org/mailman/listinfo/postgis-users</a><br>
>>> _______________________________________________<br>
>>> postgis-users mailing list<br>
>>> <a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br>
>>> <a href="http://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" target="_blank">http://lists.osgeo.org/mailman/listinfo/postgis-users</a><br>
>><br>
>><br>
>><br>
><br>
><br>
> _______________________________________________<br>
> postgis-users mailing list<br>
> <a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br>
> <a href="http://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" target="_blank">http://lists.osgeo.org/mailman/listinfo/postgis-users</a><br>
_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" target="_blank">http://lists.osgeo.org/mailman/listinfo/postgis-users</a></div></div></blockquote></div><br></div>