[postgis-devel] Faster ST_Intersection

Tue Oct 23 12:59:26 PDT 2018

вт, 23 окт. 2018 г. в 22:27, Martin Davis <mtnclimb at gmail.com>:

> That's an interesting use case.  Do you want to split the buildings across
> isochrones, or can each building just be assigned to a single isochrone?
> (I haven't actually looked at the data, so apologies if this is obvious).
>

I want to split, as this is example benchmarking data. I believe for some
other cases population can have a different source dataset, and "complex
polygon dataset" that I use isochrones for is user-supplied.

There are options of converting population to Raster, but it looks like
things like (RT_)ST_Intersection vectorize the polygon back and go to GEOS,
in their proportional/vector form. Or it there a recipe I missed?

>
> There are a zillion polygon overlay papers, and some of them even have
> implementations  8^)   JTS has used monotone chains in overlay since
> forever... sigh, shoulda published first.   I would guess that there is no
> significant further breakthrough in that paper.  Franklin's idea of
> avoiding building topology seems like the big potential win.  Shouldn't be
> too hard to code up (although needs some indexing and attention to
> robustness).
>
>
>
> On Tue, Oct 23, 2018 at 12:06 PM Darafei "Komяpa" Praliaskouski <
> me at komzpa.net> wrote:
>
>> Martin,
>>
>> The case for intersections only is the same as 25 years ago:
>> interactively summarize demographics by some area.
>>
>> In my benchmarking case I have 33575 buildings with population in
>> metadata and 233597 points in them on one side, and one minute walking
>> isochrones for three hours subdivided into 14642 parts with 5315008 points.
>> I want to have the population proportionally summarized for each walking
>> minute. No spatial viz needed in that case, just numbers.
>>
>> Raw sum(ST_Area(geom)) for both layers returns in below 100ms. Something
>> like 500ms for all the operation would be interactive enough. Today I went
>> from 280s to 17s on 8-thread i5.
>>
>> If you want to play with the data, here it is temporarily:
>> https://www.dropbox.com/s/ixg0ap8m8vpzzux/hackerspace_eta.sqld.gz?dl=0
>>  (614mb)
>> (all OdBL / from OpenStreetMap, isochrones lead to two spots where Minsk
>> Hackerspace resides).
>>
>> Other things I've found today, god bless sci-hub:
>> http://sci-hub.tw/https://ieeexplore.ieee.org/document/5380680
>> monotone-chain polygon intersection algorithm from 2009.
>>
>> Also, here: http://codeforces.com/blog/entry/2514 ("Darts") they manage
>> to calculate area of union without rebuilding the linework, interesting if
>> it can be adapted to intersection with just one sort of edges.
>>
>>
>> вт, 23 окт. 2018 г. в 21:24, Martin Davis <mtnclimb at gmail.com>:
>>
>>> Darafei, what is your use case for needing intersection areas only?
>>>
>>> It seems like this should be a common use case in spatial analysis - as
>>> long as no spatial visualization of the results are required.  (And if they
>>> are, perhaps simply computing an interior point for each resultant would
>>> suffice for visualization purposes?  Which should also be amenable to a
>>> faster algorithm - perhaps called InteriorPointClipped)
>>>
>>> On Tue, Oct 23, 2018 at 10:02 AM Darafei "Komяpa" Praliaskouski <
>>> me at komzpa.net> wrote:
>>>
>>>> Hi,
>>>>
>>>> I have two polygonal layers to intersect and need to measure area of
>>>> intersections.
>>>>
>>>> ST_Area(ST_Intersection(a.geom, b.geom)) works fine, but is rather
>>>> slow. Typical polygon on first side has <10 edges, on second 200..1000
>>>> edges.
>>>>
>>>> What are my options to make it faster? Any ideas, links, experiences to
>>>> share are welcome.
>>>>
>>>> For starters, this seems to work faster than vanilla ST_Intersection in
>>>> my case:
>>>>
>>>> create or replace function clipped_st_intersection(geom1 geometry, geom2 geometry)
>>>>   returns geometry
>>>> as $$select ST_Intersection(
>>>>               case
>>>>                 when ST_MemSize(geom1) > 160 then ST_ClipByBox2D(geom1, geom2)
>>>>                 else geom1 end, case
>>>>                 when ST_MemSize(geom2) > 160 then ST_ClipByBox2D(geom2, geom1)
>>>>                 else geom2 end)$$
>>>> language sql
>>>> immutable
>>>> strict
>>>> parallel safe;
>>>>
>>>> --
>>>> Darafei Praliaskouski
>>>> Support me: http://patreon.com/komzpa
>>>>
>>> _______________________________________________
>>>> postgis-devel mailing list
>>>> postgis-devel at lists.osgeo.org
>>>> https://lists.osgeo.org/mailman/listinfo/postgis-devel
>>>
>>> _______________________________________________
>>> postgis-devel mailing list
>>> postgis-devel at lists.osgeo.org
>>> https://lists.osgeo.org/mailman/listinfo/postgis-devel
>>
>> --
>> Darafei Praliaskouski
>> Support me: http://patreon.com/komzpa
>> _______________________________________________
>> postgis-devel mailing list
>> postgis-devel at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/postgis-devel
>
> _______________________________________________
> postgis-devel mailing list
> postgis-devel at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-devel

-- 
Darafei Praliaskouski
Support me: http://patreon.com/komzpa
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-devel/attachments/20181023/a3fbd1b5/attachment.html>