[postgis-users] Line of sight (LOS) Calculation

[Mark Wynter]

An example of how you can generate Viewsheds in PostGIS via PL/R and Grass using a single SQL query statement.

CREATE OR REPLACE FUNCTION generate_viewshed(coordxy text) RETURNS text as 
$$
library(spgrass6)
initGRASS(gisBase = "/usr/lib/grass64/", home = tempdir(), gisDbase = "/data/grassdata/", location = "aust_dem_1sec_wgs84", mapset = "postgres", SG="aust_dem_1sec", override = TRUE)
execGRASS("g.region", parameters = list(n="-34.9759927736275", s="-34.9985728860995", e="138.641896696527", w="138.619947232807", align="aust_dem_1sec"))
execGRASS("r.viewshed", parameters = list(input = "aust_dem_1sec", output = "viewshed_raster", coordinate = coordxy, obs_elev = 25, max_dist = 1000), flags = c("b", "overwrite")) 
execGRASS("r.to.vect", parameters = list(input = "viewshed_raster", output = "viewshed_vector", feature = "area"), flags = c("overwrite")) 
execGRASS("v.out.ogr", parameters = list(input = "viewshed_vector", dsn = "PG:host=localhost dbname=mydb user=postgres password=password", olayer = "viewshed_vector", format = "PostgreSQL", type = "area"), flags = c("c")) 
print('done');
$$
LANGUAGE 'plr';

You can pass more arguments into the PL/R function beside the x,y (tower coordinate) - e.g. the map set extents.

Instead of the last two lines, namely v.out.ogr and print’done', you can use the spgrass6 function “readVECT6()” which converts the grass vector to as a spatial data frame.
Then write the spatial data frame to WKT using writeWKT() in the rgeos package.   Then have the PL/R function return the WKT object - which you can then insert into your postgis table - like this….

INSERT INTO viewsheds SELECT GeomFromText(generate_viewshed('138.630922,-34.987283’))

This small modification means you can insert multiple view shed results into a single table -  all in memory and as part of a pl/pgsql function.  Whereas when I wrote the above code, v.out.ogr (certainly with spgrass6) didn’t play nicely with pl/pgsql as I couldn’t use 'v.out.ogr —overwrite' within a FOR EACH tower LOOP.  Because the table actually hadn’t been COMMITTED to disk.  As Joe Conway suggested to me at the time, I could have used DBLink to force the commits.

Its now been a few years on, so I think a PL/R function which returns WKT is a much cleaner solution.

If you’re running grass7 then you’ll need the spgrass7 library in R.

Hope this helps
Mark



> 
> I would like to calculate LOS (visibility)  from a point to an area fast as
> possible.
> For example looking from one observation point to 1000 targets that
> represents the area
> The option is to use terrain existing tools such as ARC Toolbox or using
> function in POSTGIS
> I wander if there is an experience with such task and what will be the
> execution time for a single query
>