Monday, 27 February 2017

qgis - Number of buildings between two points


I am using QGIS 2.18.13.



I have a polygon layer with building footprints, and a point layer with a given number of points. I would like to draw lines (randomly) between two points and see how many buildings the line crosses. How I can do that with QGIS?



Answer



One approach is through Virtual Layers / SpatiaLite syntax.


An example with 12 building and three points.


enter image description here


For illustration purpose add connecting lines between points with the Virtual Layer SQL:


select p1.id as id1, p2.id as id2, MakeLine(p1.geometry, p2.geometry) as geometryline
from point p1
cross join point p2
where p1.id <> p2.id


Use the Add / Edit Virtual Layer button:


enter image description here


In the QGIS Create a virtual layer dialog it looks like:


enter image description here


It will give you:


enter image description here


Add another virtual layer doing the complete query. Remember to rename the layer so you won't overwrite the first virtual layer:


with 
crossinglines(geometry) as (

select distinct MakeLine(p1.geometry, p2.geometry) as geometry
from point p1
cross join point p2
)
select distinct b.geometry
from building b
inner join crossinglines cl on st_crosses(cl.geometry, b.geometry)

In the layer control context menu Show feature count for the last virtual layer - showing the number of buildings with a line crossing.


enter image description here



The lines of the first SQL are doubled. This is handled in the second SQL.


Notice that the virtual layers are dynamic to their base layers. Adding a point will recalculate the virtual layers when they are refreshed with the QGIS refresh button.


enter image description here


enter image description here


UPDATE:


Now for adding the building id to the building being crossed add a new virtual layer with the following SQL.


select v2.geometry, b.osm_id
from virtual_layer2 v2
inner join building b on st_equals(b.geometry, v2.geometry)


Where virtual_layer2 is you joined building layer that crosses the lines. Change the column name osm_id to your building id.


enter image description here


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...