Chcę obliczyć statystyki ogniskowe dla każdej komórki rastra, w granicach określonych kryteriów.
Tło - Mam trzy binarne rastry, z których każdy reprezentuje jeden interesujący rodzaj roślinności. Chciałbym obliczyć procent pokrycia każdego rodzaju roślinności w obrębie (np.) 20 km ^ 2 dowolnej komórki w moim obszarze badań (suma / całkowita liczba komórek w sąsiedztwie). Problem polega na tym, że nie mogę użyć prostego okręgu lub kwadratu wokół każdej komórki, ponieważ gdybym to zrobił, obszar wyszukiwania użyty do obliczenia sumy obejmowałby obszary poza moim obszarem badań. Ten wyjątek jest ważny, ponieważ statystyki zostaną wykorzystane jako dane wejściowe dla modelu siedliska, a obszarów poza moim obszarem badań nie można uznać za możliwe siedliska - są zurbanizowane. Włączenie ich dałoby mi błędne statystyki. Więc co jan zależy od liczby komórek wymaganych do pokrycia obszaru równego mojej pożądanej wielkości sąsiedztwa), które spełniają moje kryteria. Kryteria są takie, że nie mieszczą się w obszarze zurbanizowanym. Myślę, że należy zastosować jakąś formę automatów komórkowych. Jednak nigdy nie pracowałem z CA.
Myślę, że chciałbym coś w rodzaju kodu startowego lub punktu we właściwym kierunku.
ODPOWIEDŹ NA KOMENTARZ PONIŻEJ:
Powiedzmy, że obliczam tę statystykę dla komórki na granicy mojej witryny badawczej. Jeśli przypiszę wszystkie obszary poza obszarem badań do zera (lub zignoruję NoData), uzyskam statystykę reprezentującą mniej więcej połowę pokrycia obszarowego, który mnie interesuje. Zatem procent pokrycia w obszarze ~ 10 km ^ 2 , zamiast obszaru 20 km ^ 2. Ponieważ badam rozmiary asortymentu domowego, jest to ważne. Okolica musi zmienić kształt, ponieważ w ten sposób zwierzę widzi / wykorzystuje krajobraz. Jeśli potrzebują 20 km ^ 2, zmienią kształt lub terytorium rodzinne. Jeśli nie zaznaczę zignoruj NoData, komórką będzie NoData - a NoData nie pomoże.
„PROGRESS” Z DNIA 24.10.2014 r
Oto kod, który do tej pory wymyśliłem, używając Shapely i Fiony:
import numpy as np
import pprint
import shapely
from shapely.geometry import*
import fiona
from fiona import collection
import math
traps = fiona.open('C:/Users/Curtis/Documents/ArcGIS/GIS_Data/occurrence/ss_occ.shp', 'r')
study_area = fiona.open('C:/Users/Curtis/Documents/ArcGIS/GIS_Data/Study_Area.shp', 'r')
for i in study_area: #for every record in 'study_area'
sa = shape(i['geometry']) #make a variable called 'sa' that is a polygon
grassland = fiona.open('C:/Users/Curtis/Documents/ArcGIS/GIS_Data/land_cover/polys_for_aa/class3_aa.shp', 'r')
pol = grassland.next()
gl = MultiPolygon([shape(pol['geometry']) for pol in grassland])
areaKM2 = 20
with traps as input:
r = (math.sqrt(areaKM2/math.pi))*1000
for point in input:
pt = shape(point['geometry'])
pt_buff = pt.buffer(r)
avail_area = pt_buff.intersection(sa).area
# works to here
while avail_area < areaKM2:
r += 10
pt_buff = pt.buffer(r)
avail_area = pt_buff.intersection(sa).area
perc_cov = pt_buff.intersection(gl).area//areaKM2
print perc_covNiestety jest NIESAMOWICIE powolny.
Odpowiedzi:
Powyższy kod jest ostateczną i niedoskonałą odpowiedzią na ten problem. W końcu pomyślałem, że najlepszym rozwiązaniem byłoby użycie okrągłego sąsiedztwa i obliczenie obszaru przecinającego mój „dostępny” obszar. (Okrągłe sąsiedztwo i tak dałoby n-najbliższym celom - więc nie trzeba się zbytnio zachowywać przy użyciu Automatów Komórkowych.) Jeśli obszar był zbyt mały, po prostu powiększałem krąg, dopóki nie był.
Działało dobrze, ale, jak zauważyłem, było bardzo wolne. Zobacz ten wątek, aby uzyskać sugestie dotyczące przyspieszenia. Maksymalizacja wydajności kodu dla Shapely . Postępowałem zgodnie z sugestiami, które doprowadziły do tego wątku Zrozumienie stosowania indeksów przestrzennych . Nie skończyło się to na zastosowaniu r-drzewa, ponieważ tak naprawdę nigdy nie skończyło się na użyciu kodu. Odkryłem, że mogę podejść do problemu z zupełnie innej perspektywy i zaoszczędzić sobie dużo czasu / energii, więc zrobiłem to. Może kiedyś to skończę ...
źródło