Расчет матрицы пересечений в ArcGIS или Exel

[Ostran2004]

Дано: шейп-файл и таблица полигонов 20 000 выделов ландшафтов, относящихся к 377 видам, отдельной колонкой прописано районирование - 16 районов (ряд видов повторяется в нескольких районах).
Можно ли подсчитать матрицу пересечений: сколько видов относящихся к тому или иному району представлены в других районах?

[trir]

https://postgis.net/docs/ST_Intersection.html

[gamm]

не, там чисто табличная операция, не пространственное пересечение. Для первой строки что-то типа
select reg,count(*) as cnt from tab where tab.type in (select type from tab where tab.reg = 1) group by reg
А дальше 16 раз union, ничего более умного в голову не приходит. Либо грузить в R/Python, и там спокойно посчитать.

[_taras_]

Ещё можно попробовать экселевский инструмент сводная таблица

[gamm]

В сводной таблице, насколько я помню, куб строится, это немного другое. Хотя могли и такое предусмотреть

[Ostran2004]

Ещё можно попробовать экселевский инструмент сводная таблица

К сожалению, в сводной таблице можно посчитать сумму площадей или количество контуров, а не уникальных значений.

Вопрос решен "ручным" методом , хотя хотелось програмно

[_taras_]

К сожалению, в сводной таблице можно посчитать сумму площадей или количество контуров, а не уникальных значений.

[/quote]

Сводная таблица - поля данных - параметр количество... Использую для подсчета уникальных значений (либра офис). Думаю, что и экселе тоже самое.

[gamm]

задача воспринималась буквально, считалось число уникальных видов, без учета сколько раз он встречался в регионе. Программу нужно рассматривать как прототип. На питоне, вроде питон в Арке есть.

Код: Выделить всё

import numpy as np

# generate random data 
n_reg    =    16
n_type   =   377
n_parcel = 20000
reg_column = 0
type_column = 1
dat_mtr = np.zeros((n_parcel,2),dtype=int)
dat_mtr[:,reg_column] = np.random.randint(1, high=n_reg+1, size=n_parcel, dtype=int)
dat_mtr[:,type_column] = np.random.randint(1, high=n_type+1, size=n_parcel, dtype=int)

# create table
reg_list = np.sort(np.unique(dat_mtr[:,reg_column]))
n_reg = len(reg_list)
result_mtr = np.zeros((n_reg,n_reg),dtype=int)
pair_mtr, pair_cnt = np.unique(dat_mtr,axis=0,return_counts=True)
for i_reg in range(n_reg):
  cur_reg = reg_list[i_reg]
  cur_reg_type_set = pair_mtr[pair_mtr[:,reg_column] == cur_reg,type_column]
  cur_ind = np.isin(pair_mtr[:,type_column],cur_reg_type_set)
  cur_reg_list, cur_cnt = np.unique(pair_mtr[cur_ind,reg_column],return_counts=True)
  cur_reg_pos = np.where(np.isin(reg_list,cur_reg_list))
  result_mtr[i_reg,cur_reg_pos] = cur_cnt
#endfor
print(result_mtr)

[Ostran2004]

К сожалению, в сводной таблице можно посчитать сумму площадей или количество контуров, а не уникальных значений.

Сводная таблица - поля данных - параметр количество... Использую для подсчета уникальных значений (либра офис). Думаю, что и экселе тоже самое.
[/quote]



Именно так и делаю, но тогда он выдает количество отдельных контуров с данным индексом...

[Ostran2004]

задача воспринималась буквально, считалось число уникальных видов, без учета сколько раз он встречался в регионе. Программу нужно рассматривать как прототип. На питоне, вроде питон в Арке есть.

Код: Выделить всё

import numpy as np

# generate random data 
n_reg    =    16
n_type   =   377
n_parcel = 20000
reg_column = 0
type_column = 1
dat_mtr = np.zeros((n_parcel,2),dtype=int)
dat_mtr[:,reg_column] = np.random.randint(1, high=n_reg+1, size=n_parcel, dtype=int)
dat_mtr[:,type_column] = np.random.randint(1, high=n_type+1, size=n_parcel, dtype=int)

# create table
reg_list = np.sort(np.unique(dat_mtr[:,reg_column]))
n_reg = len(reg_list)
result_mtr = np.zeros((n_reg,n_reg),dtype=int)
pair_mtr, pair_cnt = np.unique(dat_mtr,axis=0,return_counts=True)
for i_reg in range(n_reg):
  cur_reg = reg_list[i_reg]
  cur_reg_type_set = pair_mtr[pair_mtr[:,reg_column] == cur_reg,type_column]
  cur_ind = np.isin(pair_mtr[:,type_column],cur_reg_type_set)
  cur_reg_list, cur_cnt = np.unique(pair_mtr[cur_ind,reg_column],return_counts=True)
  cur_reg_pos = np.where(np.isin(reg_list,cur_reg_list))
  result_mtr[i_reg,cur_reg_pos] = cur_cnt
#endfor
print(result_mtr)

Спасибо, но, к сожалению для меня это сложнее китайского , так как я не инфотматик работающий в географии, а географ, использующий ГИС.

[gamm]

выложите сюда ваш файл, в текстовом виде, два столбца - район и вид, разделитель ; в архиве, чтобы был короче. Прогоним, поглядим, что получилось.

Программу легко переделать, чтобы она читала файл и писала текстовую таблицу. Если есть Питон (а при наличии Арки он точно есть), программу можно будет запустить снаружи.

[Ostran2004]

выложите сюда ваш файл, в текстовом виде, два столбца - район и вид, разделитель ; в архиве, чтобы был короче. Прогоним, поглядим, что получилось.

Программу легко переделать, чтобы она читала файл и писала текстовую таблицу. Если есть Питон (а при наличии Арки он точно есть), программу можно будет запустить снаружи.

Файл выложил, попробуем сравнить результаты...

[gamm]

Запуск (годится любой 3-й питон), вывод результата в консоль и файл

Код: Выделить всё

python test_cross_matrix.py

текст программы ниже, имена файлов в тексте. В первую строку и первый столбец добавлены номера районов (они могут быть не последовательными числами).

Код: Выделить всё

import numpy as np

# load data
in_fn = "d:/temp/test/Landscape_for_analis.txt"
out_fn = "d:/temp/test/Landscape_for_analis_mtr.txt"
reg_column = 0
type_column = 1
dat_mtr = np.loadtxt(in_fn,dtype=int,skiprows=1,delimiter=";")

# create table
reg_list = np.sort(np.unique(dat_mtr[:,reg_column]))
n_reg = len(reg_list)
result_mtr = np.zeros((n_reg,n_reg),dtype=int)
pair_mtr, pair_cnt = np.unique(dat_mtr,axis=0,return_counts=True)
for i_reg in range(n_reg):
  cur_reg = reg_list[i_reg]
  cur_reg_type_set = pair_mtr[pair_mtr[:,reg_column] == cur_reg,type_column]
  cur_ind = np.isin(pair_mtr[:,type_column],cur_reg_type_set)
  cur_reg_list, cur_cnt = np.unique(pair_mtr[cur_ind,reg_column],return_counts=True)
  cur_reg_pos = np.where(np.isin(reg_list,cur_reg_list))
  result_mtr[i_reg,cur_reg_pos] = cur_cnt
#endfor
final_mtr = np.zeros((n_reg+1,n_reg+1),dtype=int)
final_mtr[1:,1:] = result_mtr
final_mtr[0,1:] = reg_list
final_mtr[1:,0] = reg_list
np.set_printoptions(linewidth=200)
print(final_mtr)
np.savetxt(out_fn,final_mtr,fmt="%3d")

Результат

Код: Выделить всё

  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16
  1  48  38  36  17  26  36  16   2   2   0   2   2  13   4  11   2
  2  38  69  41  36  44  53  24   4  11   0   3   6  26  16  22   3
  3  36  41 129  32  39  74  38  22  38   8  17  23  50  32  50  19
  4  17  36  32  73  43  46  24   3  12   0   2   5  24  31  26   2
  5  26  44  39  43  80  55  38   3  19   0   2   5  35  30  34   4
  6  36  53  74  46  55 131  59  10  42   4   6  11  58  45  61   6
  7  16  24  38  24  38  59  91   9  40   2   8   8  64  53  62   9
  8   2   4  22   3   3  10   9  47  36  18  35  26  38  11  42  32
  9   2  11  38  12  19  42  40  36 103  18  30  30  84  52  91  38
 10   0   0   8   0   0   4   2  18  18  24  18  11  18   5  19  17
 11   2   3  17   2   2   6   8  35  30  18  41  22  31   8  36  29
 12   2   6  23   5   5  11   8  26  30  11  22  35  31  14  31  22
 13  13  26  50  24  35  58  64  38  84  18  31  31 162  76 134  46
 14   4  16  32  31  30  45  53  11  52   5   8  14  76 109  81  16
 15  11  22  50  26  34  61  62  42  91  19  36  31 134  81 222  67
 16   2   3  19   2   4   6   9  32  38  17  29  22  46  16  67  76

[Ostran2004]

Отлично, 100% совпадение результатов!

[gamm]

Отлично

"отлично" находится на кнопке с пальцев вверх