[postgis-users] Viewshed Isovist

[diplonics]

Hi,
 Thanks for your replies and advice. For your information, and possibly
critical analysis/opinion, I have written the following script that, so far,
does what I want. Probably not very elegant, generalised or ideal, and could
probably be written much more efficiently, however I don't have time at
present to improve it. Its also certainly predicated on the fact that I know
the query geometry well as I can specify exactly how the unobstructed
viewshed should be deconstructed. Basically it works by breaking the
viewshed into loads of lines from the base viewing line to the viewshed
extent line, from left to right. It then calculates the first point on each
of these lines where they intersect with the buildings model. All the points
are stored and a polygon built from them. Its accuracy is dependent on how
finely spaced the intersection lines are; the more there are the longer its
running time. So far its worked fine on my test cases but it'll be a few
hours yet to prove its generally fine as its currently processing 15,000+ of
these things.

<*****CODE START*****>
DECLARE @I, @J, @K, @LoopCrtl1, @TempVal, @TempVal2, @TempTable1,
@TempTable2, @TempTable3, @TempTable4;

/*START: - Create Table to hold all polygons that intersect with the
buildings model.*/
SET @TempTable1 = DROP TABLE IF EXISTS "Temp_Store5";
SET @TempTable1 = CREATE TABLE "Temp_Store5" (sv_id integer) WITH oids;
SET @TempTable1 = SELECT AddGeometryColumn( 'Temp_Store5', 'polygon_geom',
4326, 'POLYGON', 2 );
SET @TempTable1 = INSERT INTO "Temp_Store5" (sv_id, polygon_geom)
                              (SELECT DISTINCT(sv_id), polygon_geom FROM
svindex, model WHERE
                               ST_Intersects(svindex.polygon_geom,
model.the_geom) ORDER BY sv_id);
SET @TempTable1 = VACUUM "Temp_Store5";
/*END: - Create Table to hold all polygons that intersect with the buildings
model.*/

SET @TempTable1 = SELECT sv_id FROM "Temp_Store5";
SET @LoopCtrl1 = LINES(@TempTable1);
Print @LoopCtrl1;

/*START: - Create table to hold all the Viewshed analysed polygons*/
SET @TempTable4 = DROP TABLE IF EXISTS "viewshedpolys";
SET @TempTable4 = CREATE TABLE "viewshedpolys" (sv_id integer) WITH oids;
SET @TempTable4 = SELECT AddGeometryColumn( 'viewshedpolys', 'polygon_geom',
4326, 'POLYGON', 2 );
/*END: - Create table to hold all the Viewshed analysed polygons*/


/*While there are polygons to Viewshed analyse*/
SET @I = 0;
WHILE(@I < @LoopCtrl1)
BEGIN
   SET @TempVal = @TempTable1[@I][0];
   
   SET @TempTable2 = DROP TABLE IF EXISTS "Temp_Store6";
   SET @TempTable2 = CREATE TABLE "Temp_Store6" (id integer) WITH oids;
   SET @TempTable2 = SELECT AddGeometryColumn( 'Temp_Store6', 'line_geom',
4326, 'LINESTRING', 2 );
   SET @TempTable2 = INSERT INTO "Temp_Store6" (id, line_geom) VALUES 
                              (1, (SELECT
ST_MakeLine(ST_PointN(ST_ExteriorRing(polygon_geom), 1),
                                                     
ST_PointN(ST_ExteriorRing(polygon_geom), 2))FROM "Temp_Store5" WHERE sv_id =
'@TempVal')),
                              (2, (SELECT
ST_MakeLine(ST_PointN(ST_ExteriorRing(polygon_geom), 4),
                                                     
ST_PointN(ST_ExteriorRing(polygon_geom), 3))FROM "Temp_Store5" WHERE sv_id =
'@TempVal')),
                              (3, (SELECT
ST_MakeLine(ST_PointN(ST_ExteriorRing(polygon_geom), 1),
                                                     
ST_PointN(ST_ExteriorRing(polygon_geom), 4))FROM "Temp_Store5" WHERE sv_id =
'@TempVal'));
   SET @TempTable2 = VACUUM "Temp_Store6";

   SET @TempTable3 = DROP TABLE IF EXISTS "Temp_Store7";
   SET @TempTable3 = CREATE TABLE "Temp_Store7" (id integer) WITH oids;
   SET @TempTable3 = SELECT AddGeometryColumn( 'Temp_Store7', 'point_geom',
4326, 'POINT', 2 );
   
   /*Calculate Viewshed based on generating line 3  between lines 1 and 2 in
K increments i.e. 1%*/
   SET @J = 1;
   SET @K = 0.01;
   WHILE(@K <= 1.01)
   BEGIN   
       SET @TempTable3 = SELECT ST_PointN(ST_Intersection(model.the_geom,
temp.line_geom),1)
                         FROM "Temp_Store6" AS temp, model WHERE temp.id = 3
AND ST_Crosses(model.the_geom, temp.line_geom);
       SET @TempVal2 = LINES(@TempTable3);

       IF @TempVal2 > 0
       BEGIN
          PRINT 'INTERSECTION';
          SET @TempVal2 = @TempTable3[0][0];
          SET @TempTable3 = INSERT INTO "Temp_Store7" (id, point_geom)
VALUES (@J, '@TempVal2');          
       END
       ELSE
       BEGIN
          PRINT 'NO INTERSECTION';
          SET @TempTable3 = INSERT INTO "Temp_Store7" (id, point_geom)
VALUES 
                              (@J, (SELECT ST_PointN(line_geom, 2) FROM
"Temp_Store6" WHERE id = 3));
       END                                                                                          

       IF @K <= 1
          SET @TempTable2 = UPDATE "Temp_Store6" SET line_geom = (SELECT
ST_MakeLine(ST_Line_Interpolate_Point((SELECT line_geom FROM "Temp_Store6"
WHERE id = 1), @K),
                                                                 
ST_Line_Interpolate_Point((SELECT line_geom FROM "Temp_Store6" WHERE id =
2), @K))) WHERE id = 3;

       SET @J = @J + 1;
       SET @K = @K + 0.01;
   END
   SET @TempTable3 = INSERT INTO "Temp_Store7" (id, point_geom) VALUES 
                              (@J, (SELECT ST_PointN(line_geom, 2) FROM
"Temp_Store6" WHERE id = 1));
   SET @J = @J + 1;
   SET @TempTable3 = INSERT INTO "Temp_Store7" (id, point_geom) VALUES 
                              (@J, (SELECT ST_PointN(line_geom, 1) FROM
"Temp_Store6" WHERE id = 1));
   SET @J = @J + 1;                           
   SET @TempTable3 = INSERT INTO "Temp_Store7" (id, point_geom) VALUES 
                              (@J, (SELECT point_geom FROM "Temp_Store7"
WHERE id = 1));                           
   SET @TempTable3 = VACUUM "Temp_Store7";

   SET @TempTable4 = INSERT INTO "viewshedpolys" (sv_id, polygon_geom)
VALUES
                              (@TempVal, (SELECT
ST_MakePolygon(ST_MakeLine(point_geom)) FROM "Temp_Store7"));
                               
   SET @I = @I + 1;
END

SET @TempTable4 = VACUUM "viewshedpolys"; 
PRINT 'ALL DONE';
<*****CODE END*****>

Also for your information there is another source of viewshed analysis at
http://www.vr.ucl.ac.uk/depthmap/, which I may look at in the future,
however I don't think it'll work with a PostGIS DB easily.
-- 
View this message in context: http://www.nabble.com/Viewshed-Isovist-tp23546020p23636415.html
Sent from the PostGIS - User mailing list archive at Nabble.com.