Penggabungan Spasial yang Lebih Efisien dengan Python tanpa QGIS, ArcGIS, PostGIS, dll

Saya mencoba melakukan join spasial seperti contoh di sini: Apakah ada opsi python untuk "menggabungkan atribut berdasarkan lokasi"? . Namun, pendekatan itu tampaknya sangat tidak efisien / lambat. Bahkan menjalankan ini dengan 250 poin sederhana membutuhkan waktu hampir 2 menit dan gagal sepenuhnya pada shapefile dengan> 1.000 poin. Apakah ada pendekatan yang lebih baik? Saya ingin melakukan ini sepenuhnya dalam Python tanpa menggunakan ArcGIS, QGIS, dll.

Saya juga tertarik untuk mengetahui apakah itu mungkin untuk atribut SUM (yaitu populasi) dari semua titik yang termasuk dalam poligon dan menggabungkan jumlah itu ke poligon.

Ini kode yang saya coba konversi. Saya mendapatkan kesalahan di saluran 9:

poly['properties']['score'] += point['properties']['score']

yang mengatakan:

TypeError: jenis operan yang tidak didukung untuk + =: 'NoneType' dan 'float'.

Jika saya mengganti "+ =" dengan "=" itu berfungsi dengan baik tapi itu tidak menjumlahkan bidang. Saya juga mencoba menjadikan ini sebagai bilangan bulat tetapi gagal juga.

with fiona.open(poly_shp, 'r') as n: 
  with fiona.open(point_shp,'r') as s:
    outSchema = {'geometry': 'Polygon','properties':{'region':'str','score':'float'}}
    with fiona.open (out_shp, 'w', 'ESRI Shapefile', outSchema, crs) as output:
        for point in s:
            for poly in n:
                if shape(point['geometry']).within(shape(poly['geometry'])):  
                    poly['properties']['score']) += point['properties']['score'])
                    output.write({
                        'properties':{
                            'region':poly['properties']['NAME'],
                            'score':poly['properties']['score']},
                        'geometry':poly['geometry']})

Fiona mengembalikan kamus Python dan Anda tidak bisa menggunakan poly['properties']['score']) += point['properties']['score'])kamus.

Contoh menjumlahkan atribut menggunakan referensi yang diberikan oleh Mike T:

# read the shapefiles 
import fiona
from shapely.geometry import shape
polygons = [pol for pol in fiona.open('poly.shp')]
points = [pt for pt in fiona.open('point.shp')]
# attributes of the polygons
for poly in polygons:
   print poly['properties'] 
OrderedDict([(u'score', 0)])
OrderedDict([(u'score', 0)])
OrderedDict([(u'score', 0)])

# attributes of the points
for pt in points:
    print i['properties']
 OrderedDict([(u'score', 1)]) 
 .... # (same for the 8 points)

Sekarang, kita dapat menggunakan dua metode, dengan atau tanpa indeks spasial:

1) tanpa

# iterate through points 
for i, pt in enumerate(points):
     point = shape(pt['geometry'])
     #iterate through polygons
     for j, poly in enumerate(polygons):
        if point.within(shape(poly['geometry'])):
             # sum of attributes values
             polygons[j]['properties']['score'] = polygons[j]['properties']['score'] + points[i]['properties']['score']

2) dengan indeks R-tree (Anda dapat menggunakan pyrtree atau rtree )

# Create the R-tree index and store the features in it (bounding box)
 from rtree import index
 idx = index.Index()
 for pos, poly in enumerate(polygons):
       idx.insert(pos, shape(poly['geometry']).bounds)

#iterate through points
for i,pt in enumerate(points):
  point = shape(pt['geometry'])
  # iterate through spatial index
  for j in idx.intersection(point.coords[0]):
      if point.within(shape(multi[j]['geometry'])):
            polygons[j]['properties']['score'] = polygons[j]['properties']['score'] + points[i]['properties']['score']

Hasil dengan dua solusi:

for poly in polygons:
   print poly['properties']    
 OrderedDict([(u'score', 2)]) # 2 points in the polygon
 OrderedDict([(u'score', 1)]) # 1 point in the polygon
 OrderedDict([(u'score', 1)]) # 1 point in the polygon

Apa bedanya ?

• Tanpa indeks, Anda harus mengulangi semua geometri (poligon dan titik).
• Dengan indeks spasial pembatas ( Indeks Tata Ruang RTree ), Anda beralih hanya melalui geometri yang memiliki peluang untuk berpotongan dengan geometri Anda saat ini ('filter' yang dapat menghemat banyak perhitungan dan waktu ...)
• tetapi Indeks Spasial bukan tongkat ajaib. Ketika sebagian besar dataset harus diambil, Indeks Spasial tidak dapat memberikan manfaat kecepatan apa pun.

Setelah:

schema = fiona.open('poly.shp').schema
with fiona.open ('output.shp', 'w', 'ESRI Shapefile', schema) as output:
    for poly in polygons:
        output.write(poly)

Untuk melangkah lebih jauh, lihat Menggunakan Rtree Spatial Indexing Dengan OGR, Shapely, Fiona

Selain itu - geopanda sekarang secara opsional termasuk rtreesebagai ketergantungan, lihat repo github

Jadi, alih-alih mengikuti semua kode (sangat bagus) di atas, Anda bisa melakukan sesuatu seperti:

import geopandas
from geopandas.tools import sjoin
point = geopandas.GeoDataFrame.from_file('point.shp') # or geojson etc
poly = geopandas.GeoDataFrame.from_file('poly.shp')
pointInPolys = sjoin(point, poly, how='left')
pointSumByPoly = pointInPolys.groupby('PolyGroupByField')['fields', 'in', 'grouped', 'output'].agg(['sum'])

Untuk mendapatkan fungsionalitas yang bagus ini, pastikan untuk menginstal libspatialindex C-library terlebih dahulu

EDIT: impor paket yang diperbaiki

Gunakan Rtree sebagai indeks untuk melakukan penggabungan yang jauh lebih cepat, lalu Shapely melakukan predikat spasial untuk menentukan apakah suatu titik benar-benar berada dalam poligon. Jika dilakukan dengan benar, ini bisa lebih cepat daripada kebanyakan GIS lainnya.

Lihat contoh di sini atau di sini .

Bagian kedua dari pertanyaan Anda tentang 'SUM', gunakan dictobjek untuk mengumpulkan populasi menggunakan id poligon sebagai kuncinya. Meskipun, hal semacam ini dilakukan jauh lebih baik dengan PostGIS.

Halaman web ini menunjukkan bagaimana menggunakan pencarian point-in-polygon Bounding Box sebelum yang lebih mahal dalam kueri spasial dari Shapely.

http://rexdouglass.com/fast-spatial-joins-in-python-with-a-spatial-index/