Я моделирую распространение афалины в лагуне острова южной части Тихого океана. Я хотел бы использовать более гладкую мыльную пленку для моделирования вероятности присутствия дельфина на двумерной поверхности (долгота х широта) с учетом границ суши (поскольку, очевидно, дельфины не могут попасть на сушу).
Мне бы хотелосьзнаю, как зафиксировать границу моего учебного района (суши и прибрежных вод) до состояния, равного нулю, так как я ожидаю нулевую вероятность нахождения дельфина на суше, а также в большинстве оффшорных вод моего исследованияобласть (этот вид дельфина встречается только на мелководье лагуны). До сих пор я протестировал несколько подходов, которые я опишу ниже, но карты, предсказанные моими моделями, не соответствуют моим ожиданиям относительно того, как следует обрабатывать граничное условие.
Вот как выглядит набор данных, когдананесены на карту. Местоположения дельфинов показаны синим цветом, а места отсутствия - розовым. Земля показана черным цветом, а рифы - серой лагуной.
Шаг 1: создайте границу и завяжите объекты
Я создал полигон границы с именем soap_bnd , преобразовал его в список списков с именем bound и заполнил его узлами soap_knots . Граница включает в себя две внутренние "дыры", которые образованы двумя островками, найденными в районе исследования. Следующий код был вдохновлен: отличный пост Гэвина Симпсона https://www.fromthebottomoftheheap.net/2016/03/27/soap-film-smoothers/ и https://journals.plos.org/plosone/article/file?type=supplementary&id=info:doi/10.1371/journal.pone.0205921.s001 Я использую функцию autocrunch , созданную Дэвидом Л. Миллером (см. https://github.com/dill/soap_checker).
Позиции находятся в спроектированной системе координат UTM (следовательно, столбцы координат называются utmx и utmy)
## boundary of the soap film (one polygon with two inner loops)
# soap_bnd is initially a SpatialPolygonDataFrame
crds <- tidy(soap_bnd)
crds <- crds %>% dplyr::select(long, lat, piece)
names(crds) <- c("x", "y", "piece")
bound <- split(crds, crds$piece)
bound <-lapply(bound,`[`, c(1, 2))
nr <- seq(1, 3)
bound <-lapply(nr, function(n) as.list.data.frame(bound[[n]]))
## bound is a list of 3 lists (because 1 large polygon and 2 loops), each including 3 elements (utmx, utmy and f)
## an element f is set to a vector of zeros to set the boundary condition
bound <- lapply(nr, function(n) bound[[n]] <- c(bound[[n]], list(f = rep(0, length(bound[[n]]$x)))))
## created a grid of knots within the boundary
soap_knots <- make.soapgrid(bound[[1]][c("x", "y")], 20)
# removing the knots that overlap with the inner loops (the islands)
x <- soap_knots$x
y <- soap_knots$y
soap_knots <- soap_knots[!inSide(x = x, y = y, bnd = bound[[2]]) &
!inSide(x = x, y = y, bnd = bound[[3]]),]
# just changing names so everything matches
bound <- llply(bound, function(b){names(b)[1:2] <- c("utmx", "utmy"); return(b)})
names(soap_knots) <- c("utmx", "utmy")
# remove points too close to the boundary
crunch_ind <- autocruncher(bound, soap_knots, k = 20, xname = "utmx", yname = "utmy")
soap_knots <- soap_knots[-crunch_ind, ]
## plotting the final product
plot(soap_bnd)
points(soap_knots)
Шаг 2: запустить GAM
Модель GAM является биномиальной, поскольку я моделирую распределение местоположений дельфинов (1) с другими местоположениями, в которых дельфины отсутствовали (0).
m1.1 <- mgcv::gam(response ~
s(utmx, utmy, bs = "so", xt = list(bnd = bound)),
knots = soap_knots,
family = binomial, method = "REML",
data = target_df)
# predict the probability of presence using a raster stack called envLS_stack (which includes other environmental variables which are not used here)
pred_ras <- raster::predict(envLS_stack[[c("utmx", "utmy")]], m1.1, type = "response")
# plot the predictions - included between zero (low probability of presence) and one (high probability of presence)
plot(pred_ras, col = rev(brewer.pal(11, "Spectral")), axes = F)
Испытания и проблемы
Вот вероятность присутствия дельфина (суши и прибрежные воды белым цветом).
Я не понимаю, почему вероятность на границах оказывается самой высокой, хотя я установил границу f равной нулю.
Я попробовал другой подход, в котором я удалил f элементов в моих списках bound . Я запускаю GAM сСледующий код:
bound_nof <- llply(bound, function(d){list(utmx = d$utmx, utmy = d$utmy)})
# I tried this model with k = 20 or k = 60
m1.2 <- mgcv::gam(response ~
s(utmx, utmy, k = 20, bs = "so", xt = list(bnd = bound_nof)),
knots = soap_knots,
family = binomial, method = "REML",
data = target_df)
И вот что я получаю, соответственно, при использовании k = 20 и k = 60:
Я не могу понять, что здесь происходит? На первой карте (k = 20) это выглядит так, как будто модель сильно подобранна на одном наблюдении за дельфинами, проведенном в крайней западной части исследуемой области и дающем большой участок с высокой вероятностью присутствия там ... говоря статистически, я делаюНе думаю, что эта модель актуальна. Вторая карта (k = 60) еще более озадачивает ...
Я думаю, что эта схема может быть связана с отсутствием данных к западу и югу от моего учебного района. Необходимо ли уменьшить размер конечной области, где производится мыльная пленка? Недостатком является то, что это предотвратит любую дальнейшую экстраполяцию модели в большую область, в которой я заинтересован.