<html><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">Regina and Brent:<div><br></div><div>Wow. Thanks so much for the pointers.</div><div><br></div><div>I've now got a single query that is working excellently so far. Without indexes it runs on the entire dataset in about 5 minutes when splitting lines to ~ 150 m, which is just fine.</div><div><br></div><div>I was able to use the ST_Reverse function to make sure that any fragments at the end of the linestrings would occur at the upstream end, which I think handles my minimum requirement (for now).</div><div><br></div><div>Below is what just a slight reworking of Regina's second email example looks like.</div><div><br></div><div>Thanks!</div><div><br></div><div>Dane</div><div><br></div><div><br></div><div>------</div><div><br></div><div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(103, 56, 32); ">drop<span style="color: #23230e"> </span>table<span style="color: #23230e"> </span>if<span style="color: #23230e"> </span>exists<span style="color: #23230e"> test2;</span></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "><span style="color: #673820">create</span> <span style="color: #673820">table</span> <span style="color: #003369">test2 </span>(<span style="color: #673820">oid</span> serial,strahler <span style="color: #881350">varchar</span>, geometry geometry);</div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(103, 56, 32); ">INSERT</div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "><span style="color: #673820">INTO</span> <span style="color: #003369">test2 </span>(strahler, geometry)</div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "><span style="color: #673820">SELECT</span> strahler, </div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 110, 37); "><span style="color: #23230e"> </span>-- take a substring if the length is greater than 152.4 meters otherwise take the remainder</div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "> <span style="color: #003369">ST_Line_Substring</span>(geometry, <span style="color: #0000ff">152.4</span>*n/<span style="color: #440088">length</span>,</div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(103, 56, 32); "><span style="color: #23230e"> </span>CASE</div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "> <span style="color: #673820">WHEN</span> <span style="color: #0000ff">152.4</span>*(n+<span style="color: #0000ff">1</span>) < <span style="color: #440088">length</span> <span style="color: #673820">THEN</span> <span style="color: #0000ff">152.4</span>*(n+<span style="color: #0000ff">1</span>)/<span style="color: #440088">length</span></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "> <span style="color: #673820">ELSE</span> <span style="color: #0000ff">1</span></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "> <span style="color: #673820">END</span>) <span style="color: #673820">As</span> geometry</div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(103, 56, 32); ">FROM<span style="color: #23230e"> </span></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "> (<span style="color: #673820">SELECT</span> strahler, </div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "><span class="Apple-style-span" style="color: rgb(35, 110, 37); "><span style="color: rgb(35, 35, 14); "> </span>-- reverse the vertex order so that the upstream end of each linestring will be the remainder</span></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "> <span style="color: #003369">ST_LineMerge</span>(<span style="color: #003369">ST_Reverse</span>(ts.geometry)) <span style="color: #673820">AS</span> geometry,</div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "> <span style="color: #003369">ST_Length</span>(ts.geometry) <span style="color: #673820">As</span> <span style="color: #440088">length</span> </div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "> <span style="color: #673820">FROM</span> sshiap_lines ts</div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "> ) t </div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(0, 51, 105); "><span style="color: #673820">CROSS</span><span style="color: #23230e"> </span><span style="color: #673820">JOIN</span><span style="color: #23230e"> </span>generate_series<span style="color: #23230e">(</span><span style="color: #0000ff">0</span><span style="color: #23230e">,</span><span style="color: #0000ff">10000</span><span style="color: #23230e">) n </span></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "><span style="color: #673820">WHERE</span> n*<span style="color: #0000ff">152.4</span>/<span style="color: #440088">length</span> < <span style="color: #0000ff">1</span>;</div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "><br></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "><br></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 12px/normal 'Lucida Grande'; color: rgb(35, 35, 14); "><br></div><div><div>On Jul 9, 2008, at 10:21 PM, Paragon Corporation wrote:</div><br class="Apple-interchange-newline"><blockquote type="cite"><div> Dane,<br><br>I recall doing something along 2 and 3 before and I did it in python<br>(although I was interested in getting each measure at x distance along the<br>line so I was collecting<br>Points every x distance meters along the way rather than breaking up the<br>lines. So my translation of what I did to what you want may have a lot of<br>logical errors)<br><br>For 2 - I think your logic would look something like<br><br>2) SELECT gid, length*frac As dist, ST_Line_Substring(the_geom,startfrac,<br>endfrac) As the_geom<br>FROM (SELECT gid, 500*n/length As startfrac, CASE WHEN 500*(n+1) < length<br>THEN 500*(n+1)/length ELSE 1 END As endfrac, length, the_geom<br><span class="Apple-tab-span" style="white-space:pre"> </span>FROM (SELECT ts.gid, ST_LineMerge(ts.the_geom) As the_geom,<br>ST_Length(ts.the_geom) As length FROM line_segments ts<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>) t CROSS JOIN generate_series(0,10000) n <br><span class="Apple-tab-span" style="white-space:pre"> </span>WHERE n*500/length < 1) As segments_500<br><br><br><br>The above basically assumes your data is in some meter projection and you<br>have no line segment that has got more than 10000 (500 ft segments)<br>If you know for sure there is no way you have a line segment greater than<br>500*10000 then you can safely reduce that generate_series number and get<br>better performance.<br><br>This I did do in python by the way - and my logic is much more complicated<br>than above but hopefully I cut out enough for you to get the basic idea<br><br>3) For this I used python and the gdal and psycopg libraries. Basic logic<br>looks something like Below<br>open up a tile (my above query (maybe taking the centroid of each point I<br>actually limited to only those linesegments that intersected the tile I was<br>loading so I processed one tile at a time) - so my above query had an extra<br>join with the tile extent<br><br>Looks something like below where pt is an array of the centroid points of<br>the above query<br><br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>dataset =<br>osgeo.gdal.Open('%(atilefolder)s/%(atilefile)s'%{'atilefolder':atilefolder,'<br>atilefile':atilefile}, GA_ReadOnly)<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>geot = dataset.GetGeoTransform()<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>#0 0 row col is at top left corner of the tif so we<br>need min_x, max_y that corresponds to pos (0,0) <br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>min_x = geot[0]<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>max_y = geot[3]<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>res_x = abs(geot[1])<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>res_y = abs(geot[5])<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>print 'Origin of tile ', atilefile , '= (',min_x,<br>',',max_y,')'<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>print 'Pixel Size = (',res_x, ',',res_y,')'<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>print 'Size is<br>',dataset.RasterXSize,'x',dataset.RasterYSize, 'x',dataset.RasterCount<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>for pt in r:<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>#for each point figure out the row col<br>position in the tiff<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>gid = pt[0]<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>x = pt[2]<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>y = pt[3]<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>#print max_y - y<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>row = int(math.ceil((max_y - y)/res_y))<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>col = int(math.ceil((x - min_x)/res_x))<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>#If point falls in grid proceed to determine<br>elevation from tiff<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>#Some line segments may have measures that<br>span multiple tiles so we need to skip the measurement points that don't<br>fall in the tile<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if row > 0 and col > 0 and<br>dataset.RasterYSize > row and dataset.RasterXSize > col:<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>numpoints = numpoints + 1<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>pt_dist = pt[1]<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>#grab 1 pixel at row col position<br>in tif and return the band1 - which is the elevation<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>rd_scan = dataset.ReadRaster(col -<br>1, row - 1, 1, 1)<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>rd_area = struct.unpack('f' * 1,<br>rd_scan)<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>elev = int(rd_area[0])<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>#your update statement goes here<br><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><br><br> Hope that helps,<br>Regina<br><br><br><br>-----Original Message-----<br>From: <a href="mailto:postgis-users-bounces@postgis.refractions.net">postgis-users-bounces@postgis.refractions.net</a><br>[<a href="mailto:postgis-users-bounces@postgis.refractions.net">mailto:postgis-users-bounces@postgis.refractions.net</a>] On Behalf Of Dane<br>Blakely Springmeyer<br>Sent: Wednesday, July 09, 2008 11:25 PM<br>To: <a href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a><br>Subject: [postgis-users] Splitting/merging linework into equal<br>intervallinestrings<br><br>PostGIS users,<br><br>I have 1:24k hydrographic linework that I need process in several successive<br>steps using a postgis workflow.<br><br>I'm stuck at step 2 of this overall workflow:<br><br>1) separate all linestrings into distinct Strahler Order groups (done)<br>2) split all linework into linestrings of 500 ft (while avoiding any<br>linestring fragments smaller than 500 ft)<br>3) perform an elevation lookup to raster data to calculate the gradient of<br>each 500 ft segment and..<br>4) combine all adjacent linestrings which fall into similar gradient<br>classes.<br><br>For step 2 I'm investigating using ST_Segmentize() and/or<br>ST_Line_substring(), but I'd really appreciate some help with how to best<br>approach the problem.<br><br>ST_Segmentize seems to only insert more nodes/points into already very high<br>resolution linework (ie nodes already exist at intervals much more frequent<br>than 500ft), thus extracting linestrings from the result of ST_Segmentize<br>clearly isn't as simple as using MakeLine() between each segment.<br><br>ST_Line_Substring works on percentages which is smart, but I have yet to<br>wrap my mind around how to measure each individual line such that I can<br>translate percentage distance along a linestring into equal distance<br>intervals.<br><br>Anyone have examples or guidance?<br><br>Thanks,<br><br>Dane<br>_______________________________________________<br>postgis-users mailing list<br><a href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a><br>http://postgis.refractions.net/mailman/listinfo/postgis-users<br><br><br><br>_______________________________________________<br>postgis-users mailing list<br>postgis-users@postgis.refractions.net<br>http://postgis.refractions.net/mailman/listinfo/postgis-users<br></div></blockquote></div><br></div></body></html>