[postgis-users] Splitting/merging linework into equal intervallinestrings

[Paragon Corporation]

 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: postgis-users-bounces at postgis.refractions.net
[mailto:postgis-users-bounces at postgis.refractions.net] On Behalf Of Dane
Blakely Springmeyer
Sent: Wednesday, July 09, 2008 11:25 PM
To: postgis-users at 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 at postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users