<html><body><div style="color:#000; background-color:#fff; font-family:verdana, helvetica, sans-serif;font-size:13px"><div dir="ltr" id="yui_3_16_0_1_1426561486601_1083603"><span id="yui_3_16_0_1_1426561486601_1083606">In terms of the efficiencies underlying the query - I don't know enough to advise you. there are a few ways to "phrase" such a query, & I don't know the advantages/disadvantages. My guess is that using 2D geometrires will be faster than 3D geographies - but the casting of lon/lat coord geometries to geographies to get a measure in meters may still be faster than a transformation to the appropriate CRS (assuming you even you know one). All will work, but the answers will vary slightly....</span></div><div style="" class="" dir="ltr" id="yui_3_16_0_1_1426561486601_1085278"><br></div><div style="" class="" dir="ltr">There is also ST_Contains(), ST_Fully_Within(), ST_Distance_Sphere(), ST_Distance_Spheroid(), aqnd others that you could use.</div><div dir="ltr" id="yui_3_16_0_1_1426561486601_1085278"><br></div><div id="yui_3_16_0_1_1426561486601_1098687" dir="ltr">when you say "circle" - is this a predefined circle (polygon feature), or one created on the fly by buffering a point with a distance? If location is a point feature, forget about circles and ST_Within() and just use a Distance function.<br></div><br><div id="yui_3_16_0_1_1426561486601_1093441" dir="ltr">you can also use the ST_Distance_Sphere() or distance ST_Distance_Spheroid() functions as shown in my docs below. Or see:</div><div id="yui_3_16_0_1_1426561486601_1096650" dir="ltr"><a id="yui_3_16_0_1_1426561486601_1096653" href="http://postgis.net/docs/manual-2.1/ST_Distance_Sphere.html">http://postgis.net/docs/manual-2.1/ST_Distance_Sphere.html</a></div><div id="yui_3_16_0_1_1426561486601_1096683" dir="ltr"><a id="yui_3_16_0_1_1426561486601_1096822" href="http://postgis.net/docs/manual-2.1/ST_Distance_Spheroid.html">http://postgis.net/docs/manual-2.1/ST_Distance_Spheroid.html</a><br></div><div id="yui_3_16_0_1_1426561486601_1093440"><br></div><div dir="ltr" id="yui_3_16_0_1_1426561486601_1093163">I recommend you use ST_Distance_Spheroid(), as a small distance like 100m with lat & longs in degrees needs the extra precision of a spheroid rather than a sphere to work reliably.</div><div id="yui_3_16_0_1_1426561486601_1097704" dir="ltr">You can use the spheroid definition in the example above, or in mine below.<br></div><div id="yui_3_16_0_1_1426561486601_1097137"><br></div><div id="yui_3_16_0_1_1426561486601_1097558"><br></div><div id="yui_3_16_0_1_1426561486601_1086523" dir="ltr">Generally get the command working first, & optimise for speed if you need/want to. Getting it working reiably is the main objective.<br></div><br><div dir="ltr" id="yui_3_16_0_1_1426561486601_1086521">Your SQL below won't work. You are telling QGIS the data unit is degrees (SRID=4326) and then giving it a distance of 100, which means 100 degrees to measure against - you want 100m, not 100 degrees do you not? You need to convert your data to something with a unit of meters for distance measurement - either reproject away from 4326 to an appropriate projection or convert to geography - or use ST_Distance_Spheroid().<br></div><div id="yui_3_16_0_1_1426561486601_1086522"><br></div><div id="yui_3_16_0_1_1426561486601_1085276" dir="ltr">Without the error message I can't see anything obvious wrong with your create index command...</div><div id="yui_3_16_0_1_1426561486601_1091097" dir="ltr"><br></div><div id="yui_3_16_0_1_1426561486601_1091098" dir="ltr"><br></div><div id="yui_3_16_0_1_1426561486601_1091099" dir="ltr">I have included the relevant text from a Postgis tutorial I wrote some years ago - before geographies were implemented. You won't have the dafta to actually try these, but the may help demonstrate the principle:</div><div id="yui_3_16_0_1_1426561486601_1091100" dir="ltr"><br></div><div id="yui_3_16_0_1_1426561486601_1092467" dir="ltr">....<br></div><div id="yui_3_16_0_1_1426561486601_1091101" dir="ltr"><b id="yui_3_16_0_1_1426561486601_1091621">Querying for area measurements.</b></div><div id="yui_3_16_0_1_1426561486601_1091620" dir="ltr"><br style="" class="">Given the basic understanding of projections covered above, it should be apparent that the way to<br style="" class="">retrieve an accurate measurement depends on doing the measuring in a suitable projection. Within<br style="" class="">PostGIS you can create custom or non-standard projections. What is needed is a unique number to<br style="" class="">identify the projection, and the parameters needed for the projection transformation. PostGIS, like<br style="" class="">several other applications dealing with spatial data, uses the Proj.4 application/library to do any<br style="" class="">reprojection required. So, as well as this unique id, PostGIS also stores the parameter Proj.4 needs.<br style="" class="">EPSG (the European Petroleum Survey Group) established a set of codes and standard projections,<br style="" class="">called EPSG codes. These are all supplied with PostGIS, and are stored in the table spatial_ref_sys, any<br style="" class="">new projections needed can be added to this table.<br style="" class="">We need to create an equal area projection suitable for use around NZ. To do this, we'll add a custom<br style="" class="">equal area projection to the standard list.</div><div id="yui_3_16_0_1_1426561486601_1091477" dir="ltr"><br style="" class=""><i id="yui_3_16_0_1_1426561486601_1091476">insert into spatial_ref_sys<br style="" class="">values ( 27201,<br style="" class="">'WOODB, NIWA',<br style="" class="">27201,<br style="" class="">null,<br style="" class="">'+proj=aea +lat_1=-30 +lat_2=-50 +lat=-40 +lon_0=175 +x_0=0 +y_0=0 +ellps=WGS84<br style="" class="">+datum=WGS84 +units=m +no_defs');</i><br style="" class=""></div><div id="yui_3_16_0_1_1426561486601_1091912" dir="ltr"><br></div><div id="yui_3_16_0_1_1426561486601_1091622" dir="ltr">The SRID for NZMG is 27200, so I arbitrarily chose 27201 for this projection (it was not being used). The</div><div id="yui_3_16_0_1_1426561486601_1092746" dir="ltr"> last string entered contains the parameters required by the Proj.4 reprojection library used by PostGIS.</div><div id="yui_3_16_0_1_1426561486601_1092747" dir="ltr"><br style="" class="">Having gone through this exercise, which shows the preferred way to measure areas, PostGIS has a<br style="" class="">couple of options which make life easier, at least for linear measurements. The functions<br style="" class="">distance_sphere, distance_spheroid and length_spheroid always return values in meters. They use a<br style="" class="">sphere (or spheroid, respectively) instead of the native spheroid/datum, so are generally not as accurate,<br style="" class="">but will often meet all that is needed for a reasonably accurate distance measurement. Note that the<br style="" class="">*_spheroid functions require the spheroid radius & offset to be given. There are no equivalents for areal<br style="" class="">measurements (yet), so these will always require reprojection.<br style="" class="">To see the difference between spherical, spheroidal & projected distance measurements, try this<br style="" class="">modified version of the projection example SQL above:<br style="" class=""><br style="" class=""><i id="yui_3_16_0_1_1426561486601_1091761">select trip_code,<br style="" class="">station_no as stn,<br style="" class="">ST_length(trackline)::decimal(8,6) as native,<br style="" class="">(ST_length(trackline)*(60*1.852))::decimal(6,3) as native_km,<br style="" class="">(ST_length(ST_transform(trackline,27200))/1000)::decimal(6,3) as NZMG_km,<br style="" class="">(ST_length(ST_transform(trackline,27201))/1000)::decimal(6,3) as AEA_km,<br style="" class="">(ST_distance_sphere(startp, endp)/1000)::decimal(6,3) as sphere,<br style="" class="">(ST_length_spheroid(trackline,'SPHEROID[\"WGS84\",6378137,298.257223563]')/1000)::decimal(6,3) as spheroid<br style="" class="">from station<br style="" class="">order by trip_code,<br style="" class="">station_no;</i></div><div id="yui_3_16_0_1_1426561486601_1091471" dir="ltr"><br style="" class="">This returns six distances, cartesian degrees, cartesian degrees converted to km (1.852km/minute), km<br style="" class="">projected to NZMG, km projected to the custom AEA projections, km projected to a sphere and km<br style="" class="">projected to the WGS84 spheroid. Most of these give similar answers, all are correct answers, as far as<br style="" class="">they go. Which is most correct is debatable. NZMG has increasing distortion the further from land we<br style="" class="">go; the AEA is strictly an equal area projection, so has small inherent issues with distances, partly<br style="" class="">depending on the direction between the points and distance from the centre; spheroidal is likely to give<br style="" class="">a good estimate anywhere, but will vary depending on the spheroid used; spherical is quick but subject<br style="" class="">to some error as the world is not a sphere. The native cartesian degree value is highly dependent on the<br style="" class="">direction being measured and the distance from the equator, and is pretty definitely the least accurate.<br style="" class="">When it comes to areas, the AEA projection will be reasonably reliable & consistent around New<br style="" class="">Zealand. To get the specified areas for each stratum, and the one calculated by PostGIS, using the new<br style="" class="">AEA projection and NZMG, try this sql.</div><div id="yui_3_16_0_1_1426561486601_1091911" dir="ltr"><br style="" class=""><i id="yui_3_16_0_1_1426561486601_1091910">select t.trip_code,<br style="" class="">s.stratum,<br style="" class="">s.area_km2,<br style="" class="">(area(transform(m.geom,27201))/1000000)::decimal(7,2) as AEA_area,<br style="" class="">(area(transform(m.geom,27200))/1000000)::decimal(7,2) as NZMG_area,<br style="" class="">(area(transform(m.geom,2193))/1000000)::decimal(7,2) as NZTM_area<br style="" class="">from stratum s,<br style="" class="">stratum_def m,<br style="" class="">trip t,<br style="" class="">stratum_trip st<br style="" class="">where st.stratum_key = m.stratum_key<br style="" class="">and s.stratum=st.stratum<br style="" class="">and st.trip_code = t.trip_code<br style="" class="">and s.trip_code=t.trip_code<br style="" class="">order by trip_code,<br style="" class="">stratum;<br style="" class=""></i><br></div><div dir="ltr"><br></div><div dir="ltr"><br></div><div id="yui_3_16_0_1_1426561486601_1084855" style="font-family: verdana, helvetica, sans-serif; font-size: 13px;"> <div id="yui_3_16_0_1_1426561486601_1084854" style="font-family: HelveticaNeue, Helvetica Neue, Helvetica, Arial, Lucida Grande, sans-serif; font-size: 16px;"> <div id="yui_3_16_0_1_1426561486601_1084853" dir="ltr"> <hr id="yui_3_16_0_1_1426561486601_1084865" size="1"> <font id="yui_3_16_0_1_1426561486601_1084858" face="Arial" size="2"> <b><span style="font-weight:bold;">From:</span></b> Aaron Lewis <the.warl0ck.1989@gmail.com><br> <b id="yui_3_16_0_1_1426561486601_1084857"><span id="yui_3_16_0_1_1426561486601_1084856" style="font-weight: bold;">To:</span></b> Brent Wood <pcreso@pcreso.com> <br><b><span style="font-weight: bold;">Cc:</span></b> PostGIS Users Discussion <postgis-users@lists.osgeo.org> <br> <b><span style="font-weight: bold;">Sent:</span></b> Monday, March 23, 2015 2:19 PM<br> <b id="yui_3_16_0_1_1426561486601_1084860"><span id="yui_3_16_0_1_1426561486601_1084859" style="font-weight: bold;">Subject:</span></b> Re: [postgis-users] Need help on basic concepts, I just need really simple calculations<br> </font> </div> <div id="yui_3_16_0_1_1426561486601_1084861" class="y_msg_container"><br>Thanks Brent, I think I understand GIS a little bit now ...<br clear="none"><br clear="none">Just two more questions,<br clear="none">1. When looking for points within a "circle", is this the recommended way?<br clear="none"> Which is,<br clear="none"> Storing data in `location geometry(point, 4326)` , convert it to<br clear="none">geography then do the calculation<br clear="none"><br clear="none">2. I tried to create a GIST index and it's not working ... The<br clear="none">document says ST_DWithin uses index anyway.<br clear="none"><br clear="none">gis=# create index user_loc_gist_idx on users using gist (location);<br clear="none">gis=# VACUUM ANALYSE;<br clear="none"><br clear="none">gis=# explain SELECT * FROM users WHERE ST_DWithin (users.location,<br clear="none">st_setsrid(st_makepoint (45.3, 35.2), 4326), 100);<br clear="none"><br clear="none"><br clear="none"><br clear="none"> QUERY PLAN<br clear="none">------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------<br clear="none"> Seq Scan on users (cost=0.00..1.80 rows=1 width=146)<br clear="none"> Filter: ((location &&<br clear="none">'0103000020E610000001000000050000009A99999999594BC033333333333350C09A99999999594BC06666666666E660409A999999992962406666666666E660409A9999999929624033333333333350C09A99999999594BC033333333333350C0'::geometry)<br clear="none">AND ('0101000020E61000006666666666A646409A99999999994140'::geometry &&<br clear="none">st_expand(location, 100::double precision)) AND _st_dwithin(location,<br clear="none">'0101000020E61000006666666666A646409A99999999994140'::geometry,<br clear="none">100::double precision))<br clear="none">(2 rows)<br clear="none"><br clear="none"><br clear="none"><br clear="none"><br clear="none"><br clear="none"><br clear="none"><br clear="none">On Mon, Mar 23, 2015 at 4:14 AM, Brent Wood <<a shape="rect" ymailto="mailto:pcreso@pcreso.com" href="mailto:pcreso@pcreso.com">pcreso@pcreso.com</a>> wrote:<br clear="none">> Hi Aaron,<br clear="none">><br clear="none">> Hopefully this (simplistic) description helps.<br clear="none">><br clear="none">> Setting a SRID value for a feature tells Postgis what coordinate reference<br clear="none">> system (CRS) the coordinates are in. This includes the unit, which can be<br clear="none">> any linear unit of measurement, such as degrees (measured at the surface of<br clear="none">> the earth), meters, feet, etc. Obviously getting this right is important to<br clear="none">> calculate distances correctly.<br clear="none">><br clear="none">> Different CRS's generally apply to different parts of the globe, and include<br clear="none">> the projection parameters. Note that every projection applies some<br clear="none">> distortion - trying to map (project) a part of a spherical surface to a flat<br clear="none">> plane has this problem. There are three main types of distortion - angular<br clear="none">> (changes shapes), area and distance. Normally to measure distance using a<br clear="none">> projected CRS, you'd use an equidistant projection (which minimises distance<br clear="none">> distortions) centered near the location you are trying to measure.<br clear="none">><br clear="none">> An alternative approach is to measure it on a spheroid, in 3D space, instead<br clear="none">> of on a projected 2D space. This is basically what a Postgis geography<br clear="none">> allows you to do. But the coordinate units in a geography are degree<br clear="none">> coordinates, so you need to specify that your coordinates are unprojected<br clear="none">> lon/lat values when you use them in a geography. The SRID (Spatial Reference<br clear="none">> ID) for such a CRS is 4326.<br clear="none">><br clear="none">> In your case, try:<br clear="none">> SELECT * FROM users<br clear="none">> WHERE ST_DWithin (users.location::geography, st_setsrid(st_makepoint (146.0,<br clear="none">> 138.19), 4326)::geography, 100);<br clear="none">><br clear="none">> This assumes your location is a geometry with a SRID of 4326 (ie: the<br clear="none">> coordinate values are unprojected lon/lat degrees). It then converts this<br clear="none">> geometry to a geography datatype, which it tests against a point geometry in<br clear="none">> SRID 4326, which is also converted to a geography for the test, to see if it<br clear="none">> is within 100m. So this SQL tests geography against geography datatype.<br clear="none">><br clear="none">> If your location feature is not SRID 4326, you'll need to reproject<br clear="none">> (transform) it to 4326 to for this to work:<br clear="none">><br clear="none">> SELECT * FROM users<br clear="none">> WHERE ST_DWithin (ST_Transform(users.location,4326)::geography,<br clear="none">> st_setsrid(st_makepoint (146.0,<br clear="none">> 138.19), 4326)::geography, 100);<br clear="none">><br clear="none">> (I haven't tested it, but I think this should work)<br clear="none">><br clear="none">> Cheers<br clear="none">><br clear="none">> Bent Wood<br clear="none">><br clear="none">> ________________________________<br clear="none">> From: Aaron Lewis <<a shape="rect" ymailto="mailto:the.warl0ck.1989@gmail.com" href="mailto:the.warl0ck.1989@gmail.com">the.warl0ck.1989@gmail.com</a>><br clear="none">> To: <a shape="rect" ymailto="mailto:postgis-users@lists.osgeo.org" href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br clear="none">> Sent: Monday, March 23, 2015 12:09 AM<br clear="none">> Subject: [postgis-users] Need help on basic concepts, I just need really<br clear="none">> simple calculations<br clear="none">><br clear="none">> Hi,<br clear="none">><br clear="none">> I've been searching online for days. Trying to understand why SRID is<br clear="none">> required. So I picked some random value.<br clear="none">><br clear="none">> Now I'm need to retrieve POINTs within a circle range, e.g a circle at<br clear="none">> (146.0, 138.19) with radius of 100 meters:<br clear="none">><br clear="none">> SELECT * FROM users<br clear="none">> WHERE ST_DWithin (users.location, st_setsrid(st_makepoint (146.0,<br clear="none">> 138.19), 2600), 100);<br clear="none">><br clear="none">> It's very simple, but the result seems wrong. I have a record contains<br clear="none">> a POINT(55 43) that matches this query.<br clear="none">><br clear="none">> Anyone know what's wrong with it?<br clear="none">><br clear="none">> --<br clear="none">> Best Regards,<br clear="none">> Aaron Lewis - PGP: 0x13714D33 - <a shape="rect" href="http://pgp.mit.edu/" target="_blank">http://pgp.mit.edu/</a><br clear="none">> Finger Print: 9F67 391B B770 8FF6 99DC D92D 87F6 2602 1371 4D33<br clear="none">> _______________________________________________<br clear="none">> postgis-users mailing list<br clear="none">> <a shape="rect" ymailto="mailto:postgis-users@lists.osgeo.org" href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br clear="none">> <a shape="rect" href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><div class="qtdSeparateBR"><br><br></div><div class="yqt0975809398" id="yqtfd73100"><br clear="none">><br clear="none">><br clear="none"><br clear="none"><br clear="none"><br clear="none">-- <br clear="none">Best Regards,<br clear="none">Aaron Lewis - PGP: 0x13714D33 - <a shape="rect" href="http://pgp.mit.edu/" target="_blank">http://pgp.mit.edu/</a><br clear="none">Finger Print: 9F67 391B B770 8FF6 99DC D92D 87F6 2602 1371 4D33</div><br><br></div> </div> </div> </div></body></html>