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

[Ofer Zilberstein]

Hi,

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 

Thanks...

Ofer Zilberstein

[Mark Wynter]

Grass r.viewshed is very good and fast, and plays nicely with postgis et al. We are talking about seconds per viewshed.

You keep your DEM in grass. Then set the g.extent of the DEM raster to define your area of interest. This is very important to get fast calculation times.

Once you've called r.viewshed, you need to convert to a vector using r.to.vect. Then export to Postgis using v.out.ogr

This will give you a viewshed in postgis which you can then query in the usual way.

If you're looking to parallelise, then g.extent doesn't - each map set can only ever have 1 extent - there are ways around this.

I've run this r.viewshed to PostGIS process at scale for mobile telco operators with 20,000 plus towers x N scenarios per tower. Works well.

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

[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

[Paolo Cavallini]

Il 09/03/2015 04:16, Mark Wynter ha scritto:
> 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

Hi all.
Please note that in GRASS7, recently released, r.los has been replaced
by the (much more efficient) r.viewshed.
All the best.

--
Paolo Cavallini - www.faunalia.eu
QGIS & PostGIS courses: http://www.faunalia.eu/training.html

[Rémi Cura]

Very cool solution !

Could you give a better estimate of run time?

(example of DEM resolution, time?)

Cheers,

Rémi-C

[Mark Wynter]

Remi - the project was almost 4 years ago, we were running Grass 6.4, and at the time r.viewshed was a module addin.
The DEM was for whole of Australia, 1 second (about 30 meter intervals).
Run time for batches of 7000 towers as a single process was about 7 hours max.
1000 towers per hr, with a max distance of up to 1km (average would have been about 500m), was what we were consistently achieving. This is where g.region becomes very important to set the extent - and to get results very quickly.
So we’re talking no more than several seconds per view shed.
I was surprised by how fast it was.

I haven’t tested in Grass7 as I haven’t loaded the DEM - at some point I might get around to it - not a priority ATM.

Would be keen for someone else to have a go, and take the original PL/R code pattern, and up date it (e.g. so the function returns WKT). And we share it back to the PostGIS community.

Its a neat little example of how you can embed both R and grass directly in a PostGIS workflow. So rather than the temptation of trying to do everything in PostGIS, you now have R and grass at your fingertips, still within the comforts of PostGIS.

WDYT?

> Very cool solution !
>
> Could you give a better estimate of run time?
> (example of DEM resolution, time?)

[Mark Wynter]

Hi everyone

I’ve been talking about sharing this code for some time - well here it is. You can find the code on GitHub.

Here’s the write up

http://dimensionaledge.com/quadgrid-generation-using-recursive-sql-function/

The blog post includes a GIF animation that shows the recursive quadgrid in action, plus a 3D representation of Sydney’s population density.

Quadgrids are not limited to underlying points. You can use any value criteria, such as the length of line strings, if for say quadgridding a set of watercourses or roads.

Best regards
Mark

[Andre Mano]

Thanks for the sharing Mark - very helpful stuff!

--

..................................
André Mano
http://opussig.blogspot.com/

[Rémi Cura]

Most interesting !

I really like the LATERAL, even at the price of postgres >=9.3, it is well worth it.

I spent a lot of time to efficiently build quad tree over point cloud (plpgsql, then python),

and I must say 20 sec for 1 million points is good.

Cheers,

Rémi-C