3

I have a bunch of UK postcodes, which are basically lat/lng points. I currently map these to user-defined "areas" as follows:

  1. Create a very small box around the point.
  2. Search for areas which intersect that.
  3. If I find any, choose the one with the smallest area.
  4. If I don't, step up the box and try again.

Some notes:

  • Point 3 is crucial. One area might be "England"; another might be "Camden". I want to map to Camden, not England.
  • The areas have various geometries, though I could make them all polygons if need be.
  • The areas may overlap. There are gaps between areas.
  • A point may therefore be within one or more areas, in which case I want the smallest, or not within any, in which case I want the closest.
  • I have min and max bounds I want to apply to the size of the areas, but if need be I could do that outside the scope of this question by pre-filtering the data.

What I have currently, using MySQL and spatial indexes, is like this (edited for brevity):

    $swlat = round($loc['lat'], 2) - self::LOCATION_STEP;
    $swlng = round($loc['lng'], 2) - self::LOCATION_STEP;
    $nelat = round($loc['lat'], 2) + self::LOCATION_STEP;
    $nelng = round($loc['lng'], 2) + self::LOCATION_STEP;

    $count = 0;

    do {
        $poly = "POLYGON(($swlng $swlat, $swlng $nelat, $nelng $nelat, $nelng $swlat, $swlng $swlat))";

        # We want to find the smallest nearby or containing area.
        #
        # Exclude locations which are very large, e.g. Greater London or too small (probably just a
        # single building.
        $sql = "SELECT locations.id, locations.name, locations.maxdimension, ST_AsText(locations_spatial.geometry) AS geom,
CASE WHEN ST_Intersects(locations_spatial.geometry, ST_GeomFromText('POINT({$loc['lng']} {$loc['lat']})', {$this->dbhr->SRID()})) THEN 0
ELSE ST_Distance(locations_spatial.geometry, ST_GeomFromText('POINT({$loc['lng']} {$loc['lat']})', {$this->dbhr->SRID()})) END AS dist
FROM locations_spatial INNER JOIN `locations` ON locations.id = locations_spatial.locationid 
LEFT OUTER JOIN locations_excluded ON locations_excluded.locationid = locations.id 
WHERE 
  ST_Intersects(locations_spatial.geometry, ST_GeomFromText('$poly', {$this->dbhr->SRID()})) AND type != 'Postcode' 
  AND ST_Dimension(locations_spatial.geometry) = 2 AND locations_excluded.locationid IS NULL 
  AND locations_spatial.locationid != $id AND maxdimension < " . self::TOO_LARGE . " AND maxdimension > " . self::TOO_SMALL . " 
ORDER BY maxdimension ASC, dist ASC LIMIT 1;";
        $nearbyes = $this->dbhr->preQuery($sql);

        if (count($nearbyes) === 0) {
            $swlat -= self::LOCATION_STEP;
            $swlng -= self::LOCATION_STEP;
            $nelat += self::LOCATION_STEP;
            $nelng += self::LOCATION_STEP;
            $count++;
        }

    } while (count($nearbyes) == 0 && $count < self::LOCATION_MAX);

...but it is slow and I would like to speed it up.

What is an efficient way of doing this nowadays? More specifically, is there a way to avoid this loop approach (either using PostGIS or other tech)?

I'm familiar with the KNN problem, but I don't think this is exactly the same. I have dabbled with PostgreSQL/PostGIS, but I'm not familiar in depth with what is/isn't efficient there.

6
  • postgis will be faster than mysql here is an example of ST_Envelope gis.stackexchange.com/a/160335/276 reference postgis.net/docs/ST_Envelope.html
    – Mapperz
    Dec 14 '21 at 17:59
  • @Mapperz I've prototyped using postgis, and it is faster (as you'd expect because the spatial indexing is better). But I would still be using the loop approach, and I'm wondering whether there is a way to do this in a single query, either in PostGIS or using other tech. Edited the question to clarify. Dec 14 '21 at 18:02
  • It's a fairly big question, so here is a suggestion as a comment only: KNN is indeed the solution to remove the loop: start by selecting the nearest "region" and compute the distance between it and the point. Select any region within this distance (use st_dwithin), optionally increased by your "step size". Then identify the smallest/largest found region and proceed with the proper one.
    – JGH
    Dec 14 '21 at 18:28
  • @JGH Thanks for your reply. The difficulty I see with that is that a point might be within one or more large areas, and close to (but not within) one or more much smaller areas. In that case it is the small close areas that I want to find, not a large containing one. The loop approach tends to find the small close areas first. Do you see the problem? I can't quite work out how to express the balance between closeness and area size in a single query (and then make it well-indexed). Dec 14 '21 at 18:54
  • 1
    You need to come up with a well-specified metric calculation that takes into account both closeness and within-area. Then you can compute this for all polygons within the maximum distance, and pick the smallest.
    – dr_jts
    Dec 15 '21 at 6:17
4

Given your note

A point may therefore be within one or more areas, in which case I want the smallest, or not within any, in which case I want the closest.

things are simple: rank the ordered distance to nearby (or containing) areas, and find the smallest one in that set:

SELECT id,
       geom
FROM   (
  SELECT id,
         geom,
         RANK() OVER(ORDER BY geom <-> '<POI_WKT>'::GEOMETRY) AS _rnk
  FROM   <areas>
  -- WHERE ST_Area(geom) BETWEEN <min> AND <max>
  LIMIT  <areas_row_count/3>
) q
WHERE  _rnk = 1
ORDER BY
       ST_Area(geom)
LIMIT  1
;

Here

  • all <areas>.geom having an equal distance to the given poi will get the same RANK (this includes 0.0 for those that contain the poi), thus
    • if a poi is within one or multiple <areas>.geom, or if more than one <areas>.geom are equally close to the poi, the query returns the one with the smallest ST_Area, otherwise it returns the single closest <areas>.geom to the given poi (no matter the ST_Area)
  • because of the ORDER BY in the RANK Window function and the inner LIMIT, the operation is a fully index driven (K)NN search; the actual LIMIT value is not important as long as it is less than one third of the <areas> table row count

If that rule set is not sufficient and e.g. in the case where a poi does not lie within an <areas>.geom and you want to consider the ST_Area of the k nearest neighbors, too, to to decide which one to return, you'd need to RANK over a factor between ST_Distance & ST_Area (as mentioned by @dr_jts) - which, because of the dynamic ST_Distance partner geometry, cannot get indexed.

Update:

You could use an augmented RANK over ST_Area and ORDER BY the sum of both ranked values - using e.g. its SQRT will have a strong bias towards smaller areas; other rank amplifiers are possible, too:

SELECT  id,
        geom
FROM    (
    SELECT  id,
            area_rank,
            geom,
            SQRT(RANK() OVER(ORDER BY ST_Area(geom))) AS _arnk,
            RANK() OVER(ORDER BY geom <-> '<POI_WKT>'::GEOMETRY) AS _drnk
    FROM    <areas>
    -- WHERE ST_Area(geom) BETWEEN <min> AND <max>
    LIMIT   <areas_row_count/3>
) q
ORDER BY
        _arnk + _drnk
LIMIT   1
;
5
  • Note that I am not entirely certain what your areas are, and what your points. This query assumes a running table having Polygons of some sort.
    – geozelot
    Dec 15 '21 at 15:08
  • Thank you for introducing RANK(), which I have never fully got to grips with, and clearly should. Your first example seems very promising and is fast. The augmented version doesn't seem to be indexed. I'll do some more testing and report back on whether this does the trick. Dec 15 '21 at 18:49
  • Nice solution. Why is the one-third factor required?
    – dr_jts
    Dec 16 '21 at 17:43
  • @dr_jts that my own ballpark estimation to have the planner not seq scan the table in favor of an index lookup.
    – geozelot
    Dec 16 '21 at 20:22
  • As @geozelot says, it's hard to get this indexed if you want to consider the area. But for my use-case, I can use the indexed distance query to obtain a reasonable number of candidate locations in a subquery, and then apply some other metrics to pick one from that set. Dec 17 '21 at 10:25
1

Building on @geozelot's answer, it is possible to simulate the approach of finding locations which appear within a gradually increasing box by adding in a CASE statement with several ST_Intersects over an increasing radius. This remains nicely indexed and is therefore fast.

WITH ourpoint AS
(
 SELECT ST_MakePoint(?, ?) as p
)
SELECT
   locationid,
   name,
   ST_Area(location) AS area,
   dist,
   CASE
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.00015625), 3857)) THEN 1
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.0003125), 3857)) THEN 2
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.000625), 3857)) THEN 3
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.00125), 3857)) THEN 4
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.0025), 3857)) THEN 5
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.005), 3857)) THEN 6
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.01), 3857)) THEN 7
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.02), 3857)) THEN 8
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.04), 3857)) THEN 9
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.08), 3857)) THEN 10
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.16), 3857)) THEN 11
       WHEN ST_Intersects(location, ST_SetSRID(ST_Buffer((SELECT p FROM ourpoint),0.32), 3857)) THEN 12
   END AS intersects
  FROM (
    SELECT   locationid,
             name,
             location,
             location <-> ST_SetSRID((SELECT p FROM ourpoint), 3857) AS dist
    FROM     locations_tmp 
    WHERE    ST_Area(location) BETWEEN 0.00001 AND 0.15
    ORDER BY location <-> ST_SetSRID((SELECT p FROM ourpoint), 3857)
    LIMIT 10
) q
WHERE ST_Dimension(location) = 2
ORDER BY intersects ASC, area ASC LIMIT 1;

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.