Menemukan Polygon Neighbor atau Tetangga pada Layer

Ada beberapa kasus penggunaan dimana anda ingin menemukan semua polygon yang posisinya dekat atau tetangga pada sebuah layer. Dengan sebuah skrip phyton pendek. Kita dapat melakukannya dan banyak hal lagi di QGIS. Berikut sebuat contoh skript yang anda dapat gunakan untuk menemukan semua poligon yang mengshare batasnya dengan poligon lain pada sebuah layar dan juga menambah nama mereka pada tabel attribut. Sebagai bonus tambahan, skript ini juga menjumlahkan attribut pilihan anda dari semua poligon tetangga tadi.

Tinjauan Tugas

Untuk mendemonstrasikan bagaimana skript bekerja, kita akan menggunakan sebuah layer dari poligon negara dan menemukan negara yang berbagi batas negara. Kita juga ingin menghitung total puplasi dari negara tetangga.

Mendapatkan data

Kita akan menggunakan dataset Admin 0 - Countries dari Natural Earth.

Unduh Admin 0 - countries shapefile..

Sumber data [NATURALEARTH]

Mendpatkan skript

Unduh neighbors.py script dan simpan pada disk anda.

Prosedur

  1. Buka layer ne_10m_admin_0_countries dengan mengakses Layer ‣ Add Vector Layer.

../_images/140.png
  1. Skript menggunakan 2 field untuk melakukan pekerjaanya. Sebuah field nama dan sebuah field yang ingin anda jumlahkan. Gunakan tool Identify untuk mengklik fitur mana saja dan memeriksa attributnya. Pada kasus ini, field nama adalah NAME dan kita ingin menjumlahkan perkiraan populasi dari field POP_EST.

../_images/226.png
  1. Akses Plugins ‣ Python Console.

../_images/317.png
  1. Pada jendela Python Console , klik tombol Show Editor

../_images/412.png
  1. Pada panel Editor , klik tombol Open file dan jelajah skript unduhan anda yakni neighbors.py dan klik Open.

../_images/512.png
  1. Ketika skript sudah terbuka, mungkin anda ingin mengubah nilai _NAME_FIELD dan _SUM_FIELD untuk menyesuaikan attribut dari layer anda. Jika anda bekerja dengan layer ne_10m_admin_0_countries , anda dapat membiarkannya seperti sedia kala. Klik tombol Save pada panel Editor jika anda membuat perubahan. Sekarang klik tombol Run script untuk mengeksekusi skript.

../_images/611.png
  1. Ketika skript selesai , klik kanan layer ne_10m_admin_0_countries dan pilih Open Attribute Table.

../_images/711.png
  1. Anda akan melihat 2 attribut baru yang bernama NEIGHBORS dan SUM . Attribut ini ditambahkan oleh skript.

../_images/811.png

Di bawah adalah skript lengkap untuk referensi. Anda bisa memodifikasinya bergantung kebutuhan anda.

################################################################################
# Copyright 2014 Ujaval Gandhi
#
#This program is free software; you can redistribute it and/or
#modify it under the terms of the GNU General Public License
#as published by the Free Software Foundation; either version 2
#of the License, or (at your option) any later version.
#
#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with this program; if not, write to the Free Software
#Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#
################################################################################
from qgis.utils import iface
from PyQt4.QtCore import QVariant

# Replace the values below with values from your layer.
# For example, if your identifier field is called 'XYZ', then change the line
# below to _NAME_FIELD = 'XYZ'
_NAME_FIELD = 'NAME'
# Replace the value below with the field name that you want to sum up.
# For example, if the # field that you want to sum up is called 'VALUES', then
# change the line below to _SUM_FIELD = 'VALUES'
_SUM_FIELD = 'POP_EST'

# Names of the new fields to be added to the layer
_NEW_NEIGHBORS_FIELD = 'NEIGHBORS'
_NEW_SUM_FIELD = 'SUM'

layer = iface.activeLayer()

# Create 2 new fields in the layer that will hold the list of neighbors and sum
# of the chosen field.
layer.startEditing()
layer.dataProvider().addAttributes(
        [QgsField(_NEW_NEIGHBORS_FIELD, QVariant.String),
         QgsField(_NEW_SUM_FIELD, QVariant.Int)])
layer.updateFields()
# Create a dictionary of all features
feature_dict = {f.id(): f for f in layer.getFeatures()}

# Build a spatial index
index = QgsSpatialIndex()
for f in feature_dict.values():
    index.insertFeature(f)

# Loop through all features and find features that touch each feature
for f in feature_dict.values():
    print 'Working on %s' % f[_NAME_FIELD]
    geom = f.geometry()
    # Find all features that intersect the bounding box of the current feature.
    # We use spatial index to find the features intersecting the bounding box
    # of the current feature. This will narrow down the features that we need
    # to check neighboring features.
    intersecting_ids = index.intersects(geom.boundingBox())
    # Initalize neighbors list and sum
    neighbors = []
    neighbors_sum = 0
    for intersecting_id in intersecting_ids:
        # Look up the feature from the dictionary
        intersecting_f = feature_dict[intersecting_id]

        # For our purpose we consider a feature as 'neighbor' if it touches or
        # intersects a feature. We use the 'disjoint' predicate to satisfy
        # these conditions. So if a feature is not disjoint, it is a neighbor.
        if (f != intersecting_f and
            not intersecting_f.geometry().disjoint(geom)):
            neighbors.append(intersecting_f[_NAME_FIELD])
            neighbors_sum += intersecting_f[_SUM_FIELD]
    f[_NEW_NEIGHBORS_FIELD] = ','.join(neighbors)
    f[_NEW_SUM_FIELD] = neighbors_sum
    # Update the layer with new attribute values.
    layer.updateFeature(f)

layer.commitChanges()
print 'Processing complete.'
comments powered by Disqus

This work is licensed under a Creative Commons Attribution 4.0 International License