3

In QGIS I have two shapefiles: one point layer representing power plants (654 elements) and one polygon layer representing cities and villages connected to gas grid (1348 elements).

For every power plant I need to find all the areas, which are maximum in 10 km distance from the plant. What I can imagine in QGIS I need to find all the areas to whose border the distance from the point is smaller or equal 10 km.

As a result I want to obtain a table with the list of power plants, selected areas and their distances.

How can I make it in QGIS? I am familiar with Python, but never worked with PostGIS or GRASS.

Improve this question
2
  • Are you saying you want to buffer the power plants and use the buffer geometry to select the cities/villages whose borders touch those boundaries? Then you want the distance each power plant is from the boundary of each city/town within 10km? – DPSSpatial Mar 15 '17 at 16:14
  • Actually somehow I did not think about buffering but in this case it can work as a really easy solution to pick up interesting cities. But then I still need to calculate certain distance to the border of selected cities, yes. – sweet.sugar.cola Mar 16 '17 at 7:06
2

Next PyQGIS code does the work.

registry = QgsMapLayerRegistry.instance()

power_plants = registry.mapLayersByName('power_plants')
cities = registry.mapLayersByName('cities')

feats_pp = [ feat for feat in power_plants[0].getFeatures() ]

feats_cities = [ feat for feat in cities[0].getFeatures() ]

selected_feats = []
distances = []
city_id = []

for feat_c in feats_cities:
    for feat_pp in feats_pp:
        d = feat_c.geometry().distance(feat_pp.geometry())
        if d <= 10000:
            selected_feats.append(feat_pp)
            distances.append(d)
            city_id.append(feat_c.id())

epsg = power_plants[0].crs().postgisSrid()

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer&field=distance:double(20,5)&field=city_id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'selected_plants',
                           'memory')

prov = mem_layer.dataProvider()

for i, feat in enumerate(selected_feats):
    feat.setAttributes([i, distances[i], city_id[i]])

prov.addFeatures(selected_feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

I tested it with shapefiles (power_plants and cities) of next image:

enter image description here

After running the code at the Python Console of QGIS, it was obtained a memory layer (named 'selected_plants') whose attributes table has the list of power plants, their distances and selected areas (expressed as 'city_id').

enter image description here

Editing Note:

Creating in each shapefile its respective field name for each feature (in my example 'name_plant' and 'name_city'), next code put in these fields of memory layer the name of attribute that corresponds.

registry = QgsMapLayerRegistry.instance()

power_plants = registry.mapLayersByName('power_plants')
cities = registry.mapLayersByName('cities')

feats_pp = [ feat for feat in power_plants[0].getFeatures() ]

feats_cities = [ feat for feat in cities[0].getFeatures() ]

selected_feats = []
distances = []
plant_names = []
city_names = []

for feat_c in feats_cities:
    for feat_pp in feats_pp:
        d = feat_c.geometry().distance(feat_pp.geometry())
        if d <= 10000:
            selected_feats.append(feat_pp)
            distances.append(d)
            p_name = feat_pp.attribute('plant_name')
            plant_names.append(p_name)
            c_name = feat_c.attribute('city_name')
            city_names.append(c_name)

epsg = power_plants[0].crs().postgisSrid()

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer&field=distance:double(20,5)&field=plant_names:string(30)&field=city_names:string(30)""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'selected_plants',
                           'memory')

prov = mem_layer.dataProvider()

for i, feat in enumerate(selected_feats):
    feat.setAttributes([i, distances[i], plant_names[i], city_names[i]])

prov.addFeatures(selected_feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

After running the above code, I got:

enter image description here

Improve this answer
3
  • Thank you very much, looks like really easy solution for my problem! my question is: how is this code creating an id number of feature? In other words, how can I find out which feature in selected_plants table matches my original feature in plants layer? and the same about cities – sweet.sugar.cola Mar 16 '17 at 8:04
  • I added an 'Editing Note' where I explain how to answer this question. I hope that it helps. By the way, I assumed that you have in attributes table for each shapefile its respective features name (in my case, of course, they are fictional). – xunilk Mar 16 '17 at 8:53
  • Thank you for explanation, naming like this works. But it looks like the code is giving me only parts of the answers. It is duplicating answers, which you can see in my results here: link and giving one answer where there should be more, which you can see here. (the interested point is selected as yellow) link – sweet.sugar.cola Mar 16 '17 at 10:57
2

I have figured out the solution based on @xunilk answer, with the small change of code (in his version in my case the code was overwriting the attributes of power plants). Now it works as it should!

registry = QgsMapLayerRegistry.instance()

power_plants = registry.mapLayersByName('HydroPowerPlants_all')
cities = registry.mapLayersByName('living_areas_switzerland_with_gas_dissolved')

feats_pp = [ feat for feat in power_plants[0].getFeatures() ]

feats_cities = [ feat for feat in cities[0].getFeatures() ]

distances = []
plant_names = []
city_names = []

for feat_c in feats_cities:
    for feat_pp in feats_pp:
        d = feat_c.geometry().distance(feat_pp.geometry())
        if d <= 10000:
            distances.append(d)
            p_name = feat_pp.attribute('ZE-nr')
            plant_names.append(p_name)
            c_name = feat_c.attribute('NAMN1')
            city_names.append(c_name)


epsg = power_plants[0].crs().postgisSrid()

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer&field=distance:double(20,5)&field=plant_names:string(30)&field=city_names:string(30)""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'selected_plants',
                           'memory')

prov = mem_layer.dataProvider()
selected_feats = [QgsFeature() for _ in range(len(distances))]

for i, feat in enumerate(selected_feats):
    feat.setAttributes([i, distances[i], plant_names[i], city_names[i]])

prov.addFeatures(selected_feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.