[postgis-users] Need better strategy for a query

Stephen Woodbridge stephenwoodbridge37 at gmail.com
Thu Nov 21 07:18:58 PST 2019


A follow up on this I added st_simplifyperservetopology to smooth the curves and remove 100x the points. This is very useful after gdal_contour on a raster. 

-Steve 

Sent from my iPhone

> On Nov 19, 2019, at 10:37 PM, Stephen Woodbridge <stephenwoodbridge37 at gmail.com> wrote:
> 
> All hadn't thought of that because it's looking at number of vertices, but playing around with it something like this might work:
> 
> with segments as (
>     select st_subdivide(geom,
>         case when st_length(geom) > 0.5  -- avoid short lines
>         then (st_npoints(geom)/st_length(geom)/2.0)::integer  -- break them into roughly 0.5 deg segments
>         else case when st_npoints(geom) < 8 then 8 else st_npoints(geom) end  end ) as geom     -- don't subdivide short lines and ask for atleast 8 points
>     from c100
>     where gid < 11854891 and st_length(geom)>0    -- ignore zero length lines
> )
> select avg(st_npoints(geom)) as npoints,
>         avg(st_length(geom)) as length,
>         max(st_npoints(geom)) as mpoints,
>         max(st_length(geom)) as mlength
> from segments;
>        npoints       |      length       | mpoints |      mlength
> ---------------------+-------------------+---------+-------------------
>  74.4378379677473623 | 0.297632635835028 |     129 | 0.505588403819749
> 
> Since I'm chopping this up into a new table I can
> 1. copy all the shorter lines <= 0.5
> 2. and then add the chopped lines
> 
> Also since the lines are generated based on gdal_contour, the distance between adjacent vertices is roughly the raster resolution which is why npoint/length works/2 ~= 0.5 deg.
> 
> Thanks Paul and Martin for you ideas this should run much faster. Oh yeah, and st_linemerge is not needed in my case because I'm not using multilinestrings. It cost 200 sec out of 4945 sec to process the first 10,000 segments in my test case. Using st_subdivide the time went to 10 sec!!!!
> 
> Awesome! and many thanks.
> 
> -Steve W
> 
> 
> 
>> On 11/19/2019 4:34 PM, Paul Ramsey wrote:
>> No? http://postgis.net/docs/ST_Subdivide.html
>> 
>> P
>> 
>>>> On Nov 19, 2019, at 1:33 PM, Stephen Woodbridge <stephenwoodbridge37 at gmail.com <mailto:stephenwoodbridge37 at gmail.com>> wrote:
>>> 
>>> Hi,
>>> 
>>> I have a global dataset in SRS EPSG:4326 that are ocean depth contours. The problem is that these tend to be long and have huge bbox so spatial index does not help when trying to render them in mapserver.
>>> 
>>> I plan was to chop these into shorter segments with a command like:
>>> 
>>> bathymetry=# explain select depth,
>>>         feet,
>>>         st_linesubstring(geom, 0.5*n/length,
>>>         case when 0.5*(n+1)<length then 0.5*(n+1)/length else 1 end) as geom
>>>     from (
>>>         select a.depth,
>>>             -round(a.depth/0.3048) as feet,
>>>             st_linemerge(a.geom) as geom,
>>>             st_length(a.geom) as length
>>>         from c100 a
>>>         where st_length(a.geom)>0
>>>     ) as t
>>>     cross join generate_series(0,10000) as n
>>>     where n*0.5/length < 1;
>>>                                                 QUERY PLAN
>>> -----------------------------------------------------------------------------------------------------------
>>>  Nested Loop  (cost=0.00..549466002.55 rows=1345242667 width=1723)
>>>    Join Filter: (((((n.n)::numeric * 0.5))::double precision / st_length(a.geom)) < '1'::double precision)
>>>    ->  Seq Scan on c100 a  (cost=0.00..1279615.77 rows=4035728 width=1719)
>>>          Filter: (st_length(geom) > '0'::double precision)
>>>    ->  Function Scan on generate_series n (cost=0.00..10.00 rows=1000 width=4)
>>> (5 rows)
>>> 
>>> but running a restricted test query on 10,000 lines took 1.5 hours and with about 12M lines in the source table I estimate about 75 days to compute :(
>>> 
>>> The 0.5 in the query is for 1/2 degree (~35 mile) max length. I can adjust that to say 1 deg but I suspect that will only nominally impact the run time (running a test to see).
>>> 
>>> So is there a faster way to do this?
>>> 
>>> Would smoothing the lines first help or would the cost of smooth and then dividing be about the same?
>>> 
>>> -Steve W
>>> _______________________________________________
>>> postgis-users mailing list
>>> postgis-users at lists.osgeo.org <mailto:postgis-users at lists.osgeo.org>
>>> https://lists.osgeo.org/mailman/listinfo/postgis-users
>> 
>> 
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/postgis-users
> 


More information about the postgis-users mailing list