[postgis-users] Need help on basic concepts, I just need really simple calculations

Sun Mar 22 21:20:55 PDT 2015

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....
There is also ST_Contains(), ST_Fully_Within(), ST_Distance_Sphere(), ST_Distance_Spheroid(), aqnd others that you could use.
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.

you can also use the ST_Distance_Sphere() or distance ST_Distance_Spheroid() functions as shown in my docs below. Or see:http://postgis.net/docs/manual-2.1/ST_Distance_Sphere.htmlhttp://postgis.net/docs/manual-2.1/ST_Distance_Spheroid.html

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.You can use the spheroid definition in the example above, or in mine below.

Generally get the command working first, & optimise for speed if you need/want to. Getting it working reiably is the main objective.

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().

Without the error message I can't see anything obvious wrong with your create index command...

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:
....
Querying for area measurements.
Given the basic understanding of projections covered above, it should be apparent that the way to
retrieve an accurate measurement depends on doing the measuring in a suitable projection. Within
PostGIS you can create custom or non-standard projections. What is needed is a unique number to
identify the projection, and the parameters needed for the projection transformation. PostGIS, like
several other applications dealing with spatial data, uses the Proj.4 application/library to do any
reprojection required. So, as well as this unique id, PostGIS also stores the parameter Proj.4 needs.
EPSG (the European Petroleum Survey Group) established a set of codes and standard projections,
called EPSG codes. These are all supplied with PostGIS, and are stored in the table spatial_ref_sys, any
new projections needed can be added to this table.
We need to create an equal area projection suitable for use around NZ. To do this, we'll add a custom
equal area projection to the standard list.
insert into spatial_ref_sys
values ( 27201,
'WOODB, NIWA',
27201,
null,
'+proj=aea +lat_1=-30 +lat_2=-50 +lat=-40 +lon_0=175 +x_0=0 +y_0=0 +ellps=WGS84
+datum=WGS84 +units=m +no_defs');

The SRID for NZMG is 27200, so I arbitrarily chose 27201 for this projection (it was not being used). The last string entered contains the parameters required by the Proj.4 reprojection library used by PostGIS.
Having gone through this exercise, which shows the preferred way to measure areas, PostGIS has a
couple of options which make life easier, at least for linear measurements. The functions
distance_sphere, distance_spheroid and length_spheroid always return values in meters. They use a
sphere (or spheroid, respectively) instead of the native spheroid/datum, so are generally not as accurate,
but will often meet all that is needed for a reasonably accurate distance measurement. Note that the
*_spheroid functions require the spheroid radius & offset to be given. There are no equivalents for areal
measurements (yet), so these will always require reprojection.
To see the difference between spherical, spheroidal & projected distance measurements, try this
modified version of the projection example SQL above:

select trip_code,
station_no as stn,
ST_length(trackline)::decimal(8,6) as native,
(ST_length(trackline)*(60*1.852))::decimal(6,3) as native_km,
(ST_length(ST_transform(trackline,27200))/1000)::decimal(6,3) as NZMG_km,
(ST_length(ST_transform(trackline,27201))/1000)::decimal(6,3) as AEA_km,
(ST_distance_sphere(startp, endp)/1000)::decimal(6,3) as sphere,
(ST_length_spheroid(trackline,'SPHEROID[\"WGS84\",6378137,298.257223563]')/1000)::decimal(6,3) as spheroid
from station
order by trip_code,
station_no;
This returns six distances, cartesian degrees, cartesian degrees converted to km (1.852km/minute), km
projected to NZMG, km projected to the custom AEA projections, km projected to a sphere and km
projected to the WGS84 spheroid. Most of these give similar answers, all are correct answers, as far as
they go. Which is most correct is debatable. NZMG has increasing distortion the further from land we
go; the AEA is strictly an equal area projection, so has small inherent issues with distances, partly
depending on the direction between the points and distance from the centre; spheroidal is likely to give
a good estimate anywhere, but will vary depending on the spheroid used; spherical is quick but subject
to some error as the world is not a sphere. The native cartesian degree value is highly dependent on the
direction being measured and the distance from the equator, and is pretty definitely the least accurate.
When it comes to areas, the AEA projection will be reasonably reliable & consistent around New
Zealand. To get the specified areas for each stratum, and the one calculated by PostGIS, using the new
AEA projection and NZMG, try this sql.
select t.trip_code,
s.stratum,
s.area_km2,
(area(transform(m.geom,27201))/1000000)::decimal(7,2) as AEA_area,
(area(transform(m.geom,27200))/1000000)::decimal(7,2) as NZMG_area,
(area(transform(m.geom,2193))/1000000)::decimal(7,2) as NZTM_area
from stratum s,
stratum_def m,
trip t,
stratum_trip st
where st.stratum_key = m.stratum_key
and s.stratum=st.stratum
and st.trip_code = t.trip_code
and s.trip_code=t.trip_code
order by trip_code,
stratum;

     From: Aaron Lewis <the.warl0ck.1989 at gmail.com>
 To: Brent Wood <pcreso at pcreso.com> 
Cc: PostGIS Users Discussion <postgis-users at lists.osgeo.org> 
 Sent: Monday, March 23, 2015 2:19 PM
 Subject: Re: [postgis-users] Need help on basic concepts, I just need really simple calculations

Thanks Brent, I think I understand GIS a little bit now ...

Just two more questions,
1. When looking for points within a "circle", is this the recommended way?
    Which is,
    Storing data in `location geometry(point, 4326)` , convert it to
geography then do the calculation

2. I tried to create a GIST index and it's not working ... The
document says ST_DWithin uses index anyway.

gis=# create index user_loc_gist_idx on users using gist (location);
gis=# VACUUM ANALYSE;

gis=# explain SELECT * FROM users WHERE ST_DWithin (users.location,
st_setsrid(st_makepoint (45.3, 35.2), 4326), 100);

            QUERY PLAN
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 Seq Scan on users  (cost=0.00..1.80 rows=1 width=146)
  Filter: ((location &&
'0103000020E610000001000000050000009A99999999594BC033333333333350C09A99999999594BC06666666666E660409A999999992962406666666666E660409A9999999929624033333333333350C09A99999999594BC033333333333350C0'::geometry)
AND ('0101000020E61000006666666666A646409A99999999994140'::geometry &&
st_expand(location, 100::double precision)) AND _st_dwithin(location,
'0101000020E61000006666666666A646409A99999999994140'::geometry,
100::double precision))
(2 rows)

On Mon, Mar 23, 2015 at 4:14 AM, Brent Wood <pcreso at pcreso.com> wrote:
> Hi Aaron,
>
> Hopefully this (simplistic) description helps.
>
> Setting a SRID value for a feature tells Postgis what coordinate reference
> system (CRS) the coordinates are in. This includes the unit, which can be
> any linear unit of measurement, such as degrees (measured at the surface of
> the earth), meters, feet, etc. Obviously getting this right is important to
> calculate distances correctly.
>
> Different CRS's generally apply to different parts of the globe, and include
> the projection parameters. Note that every projection applies some
> distortion - trying to map (project) a part of a spherical surface to a flat
> plane has this problem. There are three main types of distortion - angular
> (changes shapes), area and distance. Normally to measure distance using a
> projected CRS, you'd use an equidistant projection (which minimises distance
> distortions) centered near the location you are trying to measure.
>
> An alternative approach is to measure it on a spheroid, in 3D space, instead
> of on a projected 2D space. This is basically what a Postgis geography
> allows you to do. But the coordinate units in a geography are degree
> coordinates, so you need to specify that your coordinates are unprojected
> lon/lat values when you use them in a geography. The SRID (Spatial Reference
> ID) for such a CRS is 4326.
>
> In your case, try:
> SELECT * FROM users
> WHERE ST_DWithin (users.location::geography, st_setsrid(st_makepoint (146.0,
> 138.19), 4326)::geography, 100);
>
> This assumes your location is a geometry with a SRID of 4326 (ie: the
> coordinate values are unprojected lon/lat degrees). It then converts this
> geometry to a geography datatype, which it tests against a point geometry in
> SRID 4326, which is also converted to a geography for the test, to see if it
> is within 100m. So this SQL tests geography against geography datatype.
>
> If your location feature is not SRID 4326, you'll need to reproject
> (transform) it to 4326 to for this to work:
>
> SELECT * FROM users
> WHERE ST_DWithin (ST_Transform(users.location,4326)::geography,
> st_setsrid(st_makepoint (146.0,
> 138.19), 4326)::geography, 100);
>
> (I haven't tested it, but I think this should work)
>
> Cheers
>
> Bent Wood
>
> ________________________________
> From: Aaron Lewis <the.warl0ck.1989 at gmail.com>
> To: postgis-users at lists.osgeo.org
> Sent: Monday, March 23, 2015 12:09 AM
> Subject: [postgis-users] Need help on basic concepts, I just need really
> simple calculations
>
> Hi,
>
> I've been searching online for days. Trying to understand why SRID is
> required. So I picked some random value.
>
> Now I'm need to retrieve POINTs within a circle range, e.g a circle at
> (146.0, 138.19) with radius of 100 meters:
>
> SELECT * FROM users
> WHERE ST_DWithin (users.location, st_setsrid(st_makepoint (146.0,
> 138.19), 2600), 100);
>
> It's very simple, but the result seems wrong. I have a record contains
> a POINT(55 43) that matches this query.
>
> Anyone know what's wrong with it?
>
> --
> Best Regards,
> Aaron Lewis - PGP: 0x13714D33 - http://pgp.mit.edu/
> Finger Print:  9F67 391B B770 8FF6 99DC  D92D 87F6 2602 1371 4D33
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

>
>

-- 
Best Regards,
Aaron Lewis - PGP: 0x13714D33 - http://pgp.mit.edu/
Finger Print:  9F67 391B B770 8FF6 99DC  D92D 87F6 2602 1371 4D33

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20150323/8776539a/attachment.html>