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