Liczba obiektów punktowych wewnątrz oczek siatki

OPROGRAMOWANIE:

Geomedia Professional: 06.01.11.13

PROBLEM:

Siatka kwadratowa o boku 500 m.

Siatka punktów o różnych wartościach atrybutu ranga — Fig. 1. Siatka punktów o różnych wartościach atrybutu `ranga`

Wewnątrz oczek siatki znajdują się obiekty punktowe o różnej wartości atrybutu ranga. Dla każdego oczka siatki mamy wykonać proste obliczenie arytmetyczne dane wzorem:

$$ \begin{equation} H_{jed} =-\sum _{i=1}^{m} q_{i} \cdotp \ln q_{i} \end{equation} $$ (1)

$$ \begin{equation} q_{i} =\frac{liczba\ jednostek\ i-tego\ typu}{m} \end{equation} $$ (2)

gdzie:
m - liczba wszystkich jednostek w danym oczku siatki,

np: w przypadku przedstawionym na Fig. 1:

Pole		Liczba elementów
Nr	lokalizacja	czerwone	zielone	niebieskie	wszystkie
I	lewe górne	1	2	3	6
II	prawe górne	0	0	0	0
III	lewe dolne	0	1	0	1
IV	prawe dolne	0	0	0	0

Oczko lewe górne (I):

$$ \begin{equation} H_{I} =-\sum _{i=1}^{m} q_{I} \cdotp \ln q_{I} =-\left[\left(\frac{2}{6} \cdotp \ln\frac{2}{6}\right) +\left(\frac{1}{6} \cdotp \ln\frac{1}{6}\right) +\left(\frac{3}{6} \cdotp \ln\frac{3}{6}\right)\right] \end{equation} $$

Oczko prawe górne (II):

$$ \begin{equation} H_{II} =-\sum _{i=1}^{m} q_{II} \cdotp \ln q_{II} =0 \end{equation} $$

Oczko lewe dolne (III):

$$ \begin{equation} H_{III} =-\sum _{i=1}^{m} q_{III} \cdotp \ln q_{III} =-\left(\frac{1}{1} \cdotp \ln\frac{1}{1}\right) =0 \end{equation} $$

Oczko prawe dolne (IV):

$$ \begin{equation} H_{IV} =-\sum _{i=1}^{m} q_{IV} \cdotp \ln q_{IV} =0 \end{equation} $$

ROZWIĄZANIE:

[dziękuję za olbrzymią pomoc P. Leszkowi Nowickiemu]:

W oknie dialogowym Feature Class, otwieramy kolejno klasę punktów (np: geologia_punktowe) i klasę siatki (np: siatka) i zmieniamy mazwy ich kluczy podstawowych na intuicyjne, kolejno: ID_geologia_punkty oraz ID_siatka. Krok ten nie jest konieczny ale ułatwi nam późniejsze prace.
Obliczamy iloczyn przestrzenny Spartial intersection klasy punktów (geologia_punktowe) i klasy siatki (siatka) (zapytanie wynikowe: geol_pkt_siatka). Każdy z punktów otrzymuje w ten sposób dodatkowy atrybut opisujący jego położenie względem oczek siatki (ID_siatka).
Za pomoca polecenia Aggregation, dla każdego okna siatki, dokonujemy obliczeń liczby wszystkich elementów. Załóżmy, że mamy 8 typów obiektów punktowych (rang) (np: kamieniołomy, piaskownie, jaskinie, ostańce denudacyjne itd.). Sumujemy najpierw obiekty danej rangi, a później sumy wszystkich rang. Należy zauważyć, że w podanym poniżej wzorze wyeliminowano obiekty punktowe o randze "4" i nie biorą one udziału w obliczeniach. Aggregate to summary features in: siatka; From detail features in: geol_pkt_siatka; Spartial Aggregation: touch; Attribute Aggregation: Id_siatka/Id_siatka; Output: NEW. Tworzymy atrybut funkcyjny dla każdego okna siatki obliczający liczbę wszystich punktów, tzn. m (L_j):
SUM(SUM(IF(Detail.ranga=1; 1; 0)); SUM(IF(Detail.ranga=2; 1; 0)); SUM(IF(Detail.ranga=3; 1; 0)); SUM(IF(Detail.ranga=5; 1; 0)); SUM(IF(Detail.ranga=6; 1; 0)); SUM(IF(Detail.ranga=7; 1; 0)); SUM(IF(Detail.ranga=8; 1; 0)))
Utworzony atrybut funkcyjny zapisujemy pod nazwą: SUMAelementów_Lj.
Przechodzimy do kolejnego kroku, t.j. obliczenia dla każdego oczka siatki, liczby elementów danego typu (rangi). Jeżeli nie wychodziliśmy z okna Aggregation, wystarczy w zakładce: Output; obok okna Output functional attributes, kliknąć new i wprowadzić formułę:
SUM(IF(Detail.ranga=1; 1; 0))
zapisując ją pod nową nazwą atrybutu funkcyjnego np: SUMA1. To samo robimy dla kolejnych rang punktów:
SUM(IF(Detail.ranga=2; 1; 0)) -> nazwa: SUMA2
SUM(IF(Detail.ranga=3; 1; 0)) -> nazwa: SUMA3
SUM(IF(Detail.ranga=5; 1; 0)) -> nazwa: SUMA5
SUM(IF(Detail.ranga=6; 1; 0)) -> nazwa: SUMA6
SUM(IF(Detail.ranga=7; 1; 0)) -> nazwa: SUMA7
SUM(IF(Detail.ranga=8; 1; 0)) -> nazwa: SUMA8

Utworzone atrybuty funkcyjne, co prawda nie rozwiązują zasadniczego problemu ale staną się pomocne przy testowaniu poprawności głównej formuły.
Przechodzimy do najważniejszej cześci, t.j. opracowania atrybutu funkcyjnego dla każdego oczka siatki obliczającego wartość parametru H_jed. Możemy przyjąć, że formuła będzie miała następującą postać ogólną:
-SUM( obliczenia dla rangi 1; obliczenia dla rangi 2; obliczenia dla rangi 3; obliczenia dla rangi 5; obliczenia dla rangi 6; obliczenia dla rangi 7; obliczenia dla rangi 8; )
Tworzymy obliczenia dla rangi 1. W obliczeniach wykorzystujemy działanie utworzonego wcześniej atrybutu: SUMAelementów_Lj

(CAST(SUM(IF(Detail.ranga=1; 1; 0))/Output.SUMAelementów_Lj; double))*(IF(SUM(IF(Detail.ranga=1; 1; 0))>0; (LN(CAST(SUM(IF(Detail.ranga=1; 1; 0))/Output.SUMAelementów_Lj; Double))); 0))

i dalej kolejno:
- (CAST(SUM(IF(Detail.ranga=2; 1; 0))/Output.SUMAelementów_Lj; double))*(IF(SUM(IF(Detail.ranga=2; 1; 0))>0; (LN(CAST(SUM(IF(Detail.ranga=2; 1; 0))/Output.SUMAelementów_Lj; Double))); 0))
- (CAST(SUM(IF(Detail.ranga=3; 1; 0))/Output.SUMAelementów_Lj; double))*(IF(SUM(IF(Detail.ranga=3; 1; 0))>0; (LN(CAST(SUM(IF(Detail.ranga=3; 1; 0))/Output.SUMAelementów_Lj; Double))); 0))
- (CAST(SUM(IF(Detail.ranga=5; 1; 0))/Output.SUMAelementów_Lj; double))*(IF(SUM(IF(Detail.ranga=5; 1; 0))>0; (LN(CAST(SUM(IF(Detail.ranga=5; 1; 0))/Output.SUMAelementów_Lj; Double))); 0))
- (CAST(SUM(IF(Detail.ranga=6; 1; 0))/Output.SUMAelementów_Lj; double))*(IF(SUM(IF(Detail.ranga=6; 1; 0))>0; (LN(CAST(SUM(IF(Detail.ranga=6; 1; 0))/Output.SUMAelementów_Lj; Double))); 0))
- (CAST(SUM(IF(Detail.ranga=7; 1; 0))/Output.SUMAelementów_Lj; double))*(IF(SUM(IF(Detail.ranga=7; 1; 0))>0; (LN(CAST(SUM(IF(Detail.ranga=7; 1; 0))/Output.SUMAelementów_Lj; Double))); 0))
- (CAST(SUM(IF(Detail.ranga=8; 1; 0))/Output.SUMAelementów_Lj; double))*(IF(SUM(IF(Detail.ranga=8; 1; 0))>0; (LN(CAST(SUM(IF(Detail.ranga=8; 1; 0))/Output.SUMAelementów_Lj; Double))); 0))
Łącząc utworzone elementy w sumę, tworzymy atrybut funkcyjny obliczający H_jed:
-SUM(CAST(SUM(IF(Detail.ranga=1; 1; 0))/Output.SUMAelementów_Lj; Double)*IF(SUM(IF(Detail.ranga=1; 1; 0))>0; LN(CAST(SUM(IF(Detail.ranga=1; 1; 0))/Output.SUMAelementów_Lj; Double)); 0); CAST(SUM(IF(Detail.ranga=2; 1; 0))/Output.SUMAelementów_Lj; Double)*IF(SUM(IF(Detail.ranga=2; 1; 0))>0; LN(CAST(SUM(IF(Detail.ranga=2; 1; 0))/Output.SUMAelementów_Lj; Double)); 0); CAST(SUM(IF(Detail.ranga=3; 1; 0))/Output.SUMAelementów_Lj; Double)*IF(SUM(IF(Detail.ranga=3; 1; 0))>0; LN(CAST(SUM(IF(Detail.ranga=3; 1; 0))/Output.SUMAelementów_Lj; Double)); 0); CAST(SUM(IF(Detail.ranga=5; 1; 0))/Output.SUMAelementów_Lj; Double)*IF(SUM(IF(Detail.ranga=5; 1; 0))>0; LN(CAST(SUM(IF(Detail.ranga=5; 1; 0))/Output.SUMAelementów_Lj; Double)); 0); CAST(SUM(IF(Detail.ranga=6; 1; 0))/Output.SUMAelementów_Lj; Double)*IF(SUM(IF(Detail.ranga=6; 1; 0))>0; LN(CAST(SUM(IF(Detail.ranga=6; 1; 0))/Output.SUMAelementów_Lj; Double)); 0); CAST(SUM(IF(Detail.ranga=7; 1; 0))/Output.SUMAelementów_Lj; Double)*IF(SUM(IF(Detail.ranga=7; 1; 0))>0; LN(CAST(SUM(IF(Detail.ranga=7; 1; 0))/Output.SUMAelementów_Lj; Double)); 0); CAST(SUM(IF(Detail.ranga=8; 1; 0))/Output.SUMAelementów_Lj; Double)*IF(SUM(IF(Detail.ranga=8; 1; 0))>0; LN(CAST(SUM(IF(Detail.ranga=8; 1; 0))/Output.SUMAelementów_Lj; Double)); 0))
Atrybut nazywamy: Hjed.
Na koniec, choć to oczywiście nie jest niekonieczne, dla odróżnienia sytuacji, w której w oczku siatki brak jakichkolwiek elementów punktowych (H_jed=0), od sytuacji, w której mamy dowolną ilość punktów tej samej rangi (H_jed=0), możemy dołożyć warunek:
IF(Output.SUMAelementów_Lj=0; 99; formuła atrybutu funkcyjnego Hjed)
Teraz w oczkach siatki, w których brak jakichkolwiek elementów, w wyniku działania atrybutu, pojawi się wartość "99".

Powiązane tematy:

GIS Geomedia Professional georóżnorodność geodiversity Lj liczba punktów