Obliczanie statystyk ogniskowych dla specjalnego sąsiedztwa?


18

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_cov

Niestety jest NIESAMOWICIE powolny.


1
to ciekawy problem. Możesz ustawić wszystkie komórki poza obszarem badawczym na NoData, ale nie wiem, jak uda ci się uzyskać sąsiedztwo w celu dostosowania i utrzymania tego samego rozmiaru 20 km kw. (Musiałoby to zmienić kształt).
jbchurchill

@ CSB jbchurchill ma rację, najlepszą rzeczą do zrobienia jest przypisanie wartości NoData poza obszarem studiów. Narzędzie Focal Stats może odpowiednio traktować te wartości węzłów. Zobacz „Przetwarzanie komórek NoData” tutaj resources.arcgis.com/en/help/main/10.1/index.html#//…
WhiteboxDev

@WhiteboxDev - Twoja sugestia nie rozwiąże mojego problemu. Przeredaguję powyższe i wyjaśnię, dlaczego to nie zadziała.
CSB

Czy widziałeś ten post, który omawia korzystanie z Focal Statistics ze zmiennym promieniem ( gis.stackexchange.com/questions/34306/… )? Wydaje się, że to twój problem - komórki na krawędzi powinny mieć duży promień i rozważać tylko półkoliste sąsiedztwo. Oczywiście, w zależności od wielkości komórki, może być konieczne utworzenie wielu, wielu rastrów do wyboru.
łyko

1
@CSB Zamierzasz uzyskać efekty krawędziowe bez względu na to, czy używasz NoData i zmniejszonej dzielnicy, czy zmienisz kształt / położenie swojej dzielnicy, aby zapewnić jej rozmiar. Przynajmniej w tym pierwszym przypadku nie będziesz nadmiernie próbkował / reprezentował danych o krawędzi w nieprzejrzysty sposób. Jest to część niesławnego problemu modyfikowalnej jednostki powierzchni.
WhiteboxDev

Odpowiedzi:


0

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ę ...


Czytanie kodu wygląda na to, że istnieje spora szansa, że ​​użycie indeksu przestrzennego może przyspieszyć kod, często dramatycznie. Robienie skrzyżowań z takim wielokątem jest bardzo wolne.
Snorfalorpagus,

@ Snorfalorpagus Tak, jeśli spojrzysz na odpowiedź, wskazuję dwa inne wątki związane z tym pytaniem. Oba omawiają za pomocą indeksu przestrzennego.
CSB,
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.