Попарное сравнение растра в R: альтернатива для -l oop? - PullRequest
Новая
       87

Попарное сравнение растра в R: альтернатива для -l oop?

2 голосов
/ 28 января 2020

Как эффективно сравнить пары растров распределения (raster слои, содержащие только 0 и 1)? Мне нужно измерить сходство между ~ 6500 отдельных глобальных растров. Istat из SDMTools должен сделать эту работу.

Вот мой код: 

library(raster)
library(SDMTools)

Создание воспроизводимых примеров данных: растры со значениями 0 и 1 

# first raster
r1 <- raster(nrow=1800, ncol=3600, xmn=-18000000, xmx=18000000, ymn=-9000000, ymx=9000000, 
             crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs", 
             resolution=10000, vals=0)
r2 <- raster(nrow=1800, ncol=3600, xmn=-18000000, xmx=0, ymn=0, ymx=9000000, 
             crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs", 
             resolution=10000, vals=2)
r12 <- mosaic(r1, r2, fun=mean)

# second raster
r3 <- raster(nrow=1800, ncol=3600, xmn=-18000000, xmx=18000000, ymn=-9000000, ymx=9000000, 
             crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs", 
             resolution=10000, vals=0)
r4 <- raster(nrow=1800, ncol=3600, xmn=-12000000, xmx=15000000, ymn=2000000, ymx=3000000, 
             crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs", 
             resolution=10000, vals=2)
r34 <- mosaic(r3, r4, fun=mean)

Список растров 

files_list <- list(r12, r34)

Создать пустую матрицу для заполнения данными из l oop

ras_comp <- matrix(NA, nrow=length(files_list), ncol=length(files_list))
ras_comp
# label rows and columns of matrix
rownames(ras_comp) <- c("r12", "r34")
colnames(ras_comp) <- c("r12", "r34")
ras_comp

L oop, чтобы сравнить все возможные пары матриц / растров 

for (i in 1:length(files_list)) {
  # load raster i
  ras_i <- as.matrix(files_list[[i]])

  for (j in 1:length(files_list)) {
    # load raster j
    ras_j <- as.matrix(files_list[[j]])

    # compare both rasters
    ras_Istat <- Istat(ras_i, ras_j, old=F)

    # write value into matrix
    ras_comp[i,j] <- ras_Istat
  }
}

Проверить итоговую матрицу 

ras_comp
> ras_comp
          r12       r34
r12 1.0000000 0.1814437
r34 0.1814437 1.0000000

Преобразование растров в матрицы с помощью as.matrix значительно сокращает время расчета, и мне нужна итоговая финальная таблица, но выполнение этой операции для тысяч растров занимает целую вечность. Как оптимизировать код для более эффективного сравнения растров?

1 Ответ

1 голос
/ 28 января 2020

Istat выполняет кучу тестов и масштабирования перед простым вычислением. Если вы знаете, что эти тесты пройдены, вы можете выполнять масштабирование как единое целое и работать с масштабированными значениями. Он делает: 

if (length(which(dim(x) == dim(y))) != 2) 
    stop("matrix / raster objects must be of the same extent")
if (min(c(x, y), na.rm = T) < 0) 
    stop("all values must be positive")

Затем он проверяет, где оба растра являются "конечными", что включает NA значения: 

pos = which(is.finite(x) & is.finite(y))

Затем он вычисляет масштабированные значения растров: 

px = x[pos]/sum(x[pos])
py = y[pos]/sum(y[pos])
H = sqrt(sum((sqrt(px) - sqrt(py))^2))

Если old=FALSE как у вас, то он возвращает: 

    return(1 - (H^2)/2)

> Istat(r12,r34)
[1] 0.1814437

Если я вырежу тесты и напишу функцию, которая работает с масштабированными значениями, я могу получить ее до : 

fIstat = function(px,py){
    1 - (sum((sqrt(px) - sqrt(py))^2))/2
}

Проверка путем масштабирования растров и запуска: 

r12px = r12[]/sum(r12[])
r34px = r34[]/sum(r34[])
fIstat(r12px, r34px)
# [1] 0.1814437

То же значение. Отлично, но так ли это быстрее? 

> microbenchmark(fIstat(r12px, r34px), Istat(r12,r34))
Unit: milliseconds
                 expr        min         lq       mean     median         uq
 fIstat(r12px, r34px)   49.95867   78.28649   78.10863   79.45235   80.85234
      Istat(r12, r34) 1084.84825 1181.31116 1217.64122 1212.93180 1263.50811
       max neval
  106.6803   100
 1349.0239   100

Да, по большому счету.

Итак ... если в ваших данных отсутствуют пропущенные значения или бесконечности , создайте files_list из этих масштабированных растровых значений, назовите my fIstat, только l oop над верхний треугольник, и вы должны ускорить это в 10 раз.

...