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
>
_______________________________________________
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
