Geoadat-csomag előkészítése

Töltsük be a következő rétegeket:

  • Települések
  • Járások
  • Úthálózat (OSM)

Oldjuk meg a következő feladatokat programozás révén:

  1. Iteráljunk végig az adminisztratív egységeken, miközben:

    1. Készítsünk mappát minden település nevével
    2. Mentsük bele az egyes elemeket - adminisztratív egységeket - a megfelelő mappába
    3. Vágjuk ki az adott adminsztratív egységre eső úthálózatot

Szelektáljuk ki azokat a településeket, amelyek a Tapolcai járás területébe esnek:

Vektor / Kutató eszközök / Kiválasztás pozíció alapján (Benne vannak)

1. ábra: Tapolcai járás településeinek szelekciója.
1. ábra: Tapolcai járás településeinek szelekciója.

Mentsük a szelektált elemeket egy külön rétegre.

Ezt követően kezdjük el írni a kódot. Töltsük be a településeket, illetve az OpenStreetMap úthálózatát:

path = "d:/indikatrix-org/geoadat-csomag-elokeszitese/" 

vector_filename1 = "tapolcai_telepulesek_wgs84"
vector_filename2 = "gis_osm_roads_free_1"

vectorlayer1 = QgsVectorLayer(path + vector_filename1 + ".shp", vector_filename1, "ogr")
vectorlayer2 = QgsVectorLayer(path + vector_filename2 + ".shp", vector_filename2, "ogr")

Bővítsük a kódot; gyűjtsük egy listába a településeket - ezeken fogunk végigmenni. Tesztként nyomtassuk ki a települések neveit:

idx = vectorlayer1.fields().lookupField('name')    
settlements = vectorlayer1.getFeatures()
for settlement in settlements:
    print(settlement.attributes()[idx])

Hozzunk létre mappákat a munkakönyvtárunkba minden településnek:

import os

path = "d:/indikatrix-org/geoadat-csomag-elokeszitese/" 

vector_filename1 = "tapolcai_telepulesek_wgs84"
vector_filename2 = "gis_osm_roads_free_1"

vectorlayer1 = QgsVectorLayer(path + vector_filename1 + ".shp", vector_filename1, "ogr")
vectorlayer2 = QgsVectorLayer(path + vector_filename2 + ".shp", vector_filename2, "ogr")

idx = vectorlayer1.fields().lookupField('name')    
settlements = vectorlayer1.getFeatures()
for settlement in settlements:
    settlement_name = settlement.attributes()[idx]
    print(settlement_name)
    dirname = path + settlement name
    try:  
      os.mkdir(dirname)
    except OSError:  
      print("Creation of the directory %s failed" % dirname)
    else:  
      print("Successfully created the directory %s " % dirname)

Megjegyzés: A településnevek furcsa karaktereit úgy javíthatjuk ki, ha előtte a réteg kódolását UTF-8-ra állítjuk, majd átmentjük egy új rétegbe, szintén UTF-8 kódolással.

Bővítsük tovább a kódot, mentsük ki a településeket külön rétegekre, illetve töltsük be őket a grafikus felületre:

import os

path = "d:/indikatrix-org/geoadat-csomag-elokeszitese/" 

vector_filename1 = "tapolcai_telepulesek_wgs84"
vector_filename2 = "gis_osm_roads_free_1"

vectorlayer1 = QgsVectorLayer(path + vector_filename1 + ".shp", vector_filename1, "ogr")
vectorlayer2 = QgsVectorLayer(path + vector_filename2 + ".shp", vector_filename2, "ogr")

crs = vectorlayer1.crs()

idx = vectorlayer1.fields().lookupField('NAME')    
settlements = vectorlayer1.getFeatures()
for settlement in settlements:
    id = [settlement.id()]
    vectorlayer1.select(id)
    settlement_name = settlement.attributes()[idx]
    print(settlement_name)
    dirname = path + settlement_name
    try:  
      os.mkdir(dirname)
    except OSError:  
      print("Creation of the directory %s failed" % dirname)
    else:  
      print("Successfully created the directory %s " % dirname)
    QgsVectorFileWriter.writeAsVectorFormat( vectorlayer1, dirname + "/adm_area.shp", "system", crs, "ESRI Shapefile", 1)
    vectorlayer1.removeSelection()
    templayer = QgsVectorLayer(dirname + "/adm_area.shp", settlement_name, "ogr")
    QgsProject.instance().addMapLayer(templayer)

2. ábra: Tapolcai járás települései külön rétegeken, véletlen színekkel.
2. ábra: Tapolcai járás települései külön rétegeken, véletlen színekkel.

Végül bővítsük úgy a kódot, hogy a ciklus során készüljön el az úthálozatnak az adott településre vonatkozó kivágata:

import os
import processing

path = "d:/indikatrix-org/geoadat-csomag-elokeszitese/" 

vector_filename1 = "tapolcai_telepulesek_wgs84"
vector_filename2 = "gis_osm_roads_free_1"

vectorlayer1 = QgsVectorLayer(path + vector_filename1 + ".shp", vector_filename1, "ogr")
vectorlayer2 = QgsVectorLayer(path + vector_filename2 + ".shp", vector_filename2, "ogr")

crs = vectorlayer1.crs()

idx = vectorlayer1.fields().lookupField('NAME')    
settlements = vectorlayer1.getFeatures()
for settlement in settlements:
    id = [settlement.id()]
    vectorlayer1.select(id)
    settlement_name = settlement.attributes()[idx]
    print(settlement_name)
    dirname = path + settlement_name
    try:  
      os.mkdir(dirname)
    except OSError:  
      print("Creation of the directory %s failed" % dirname)
    else:  
      print("Successfully created the directory %s " % dirname)
    QgsVectorFileWriter.writeAsVectorFormat( vectorlayer1, dirname + '/' + "adm_area" + ".shp", "system", crs, "ESRI Shapefile", 1)
    vectorlayer1.removeSelection()
    templayer = QgsVectorLayer(dirname + "/adm_area.shp", settlement_name, "ogr")
    QgsProject.instance().addMapLayer(templayer)
    #processing.algorithmHelp("native:clip")
    parameters = {'INPUT': vectorlayer2,
              'OVERLAY': dirname + "/adm_area.shp",
              'OUTPUT': dirname + "/road_clip.shp"
              }
    processing.run("native:clip", parameters)

3. ábra: Tapolcai járás települései, illetve az első három település úthálózata: Ábrahámhegy, Badacsonytomaj, Badacsonytördemic.
3. ábra: Tapolcai járás települései, illetve az első három település úthálózata: Ábrahámhegy, Badacsonytomaj, Badacsonytördemic

A kód továbbfejlesztését a kedves olvasóra bízom, elkészíthetők még pl.:

  • Tematikus kivágatok készítése:

    expr = QgsExpression( "fclass='primary'" )
    primary_roads = vectorlayer2.getFeatures(QgsFeatureRequest(expr))
  • A rétegek Egységes Országos Vetületbe (EOV) történő átszámítása:

    crs = QgsCoordinateReferenceSystem(23700, 2)

© Dr. Wirth Ervin