[postgis-users] Need better strategy for a query

[Stephen Woodbridge]

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