R: Рассчитать порог, диапазон и самородок из растрового объекта

Мне нужно рассчитать порог, диапазон и самородок из растрового слоя. Я изучил пакеты gstat, usdm, в которых можно создать вариограмму, однако я не смог найти функцию, которая с учетом растрового слоя будет оценивать эти параметры. В большинстве функций эти параметры должны быть определены, например. кригинг.

У меня есть слои растровых данных для разной высоты, похожие на введите здесь описание изображения

Я хотел бы получить порог, самородок и диапазон из параметров вариограммы, подогнанных к этим слоям данных, чтобы создать график, подобный этому: введите здесь описание изображения

Исходные слои данных доступны здесь в виде многоканального TIFF. Вот рисунок из этот документ, который дополнительно иллюстрирует концепцию.

введите описание изображения здесь


person Arihant    schedule 16.03.2016    source источник
comment
Приведите пример. Gstat работает с пространственными точками, если вы хотите работать с сетками, вы можете прочитать их как фрейм данных пространственной сетки или вы также можете использовать функцию интерполяции в растровом пакете. растровый пакет предоставляет несколько примеров.   -  person Geo-sp    schedule 17.03.2016
comment
Я отредактировал свой вопрос, чтобы добавить больше деталей.   -  person Arihant    schedule 17.03.2016
comment
Я не думаю, что это проблема пространственной статистики. Я думаю, вы должны оценить полудисперсию по слоям и построить график относительно средних значений каждого слоя. Хотя не уверен, что это может помочь.   -  person Geo-sp    schedule 17.03.2016
comment
когда я рисую это, высота колеблется от 0 до 1, а форма кривой противоположна тому, что вы показываете здесь; показывая самую низкую полудисперсию на высоте 0,5.   -  person Geo-sp    schedule 17.03.2016
comment
Эти изображения предназначены для иллюстративных целей. Данные имеют 99 слоев, каждый слой расположен на высоте 0,5 метра. График зависимости фактических данных от высоты будет показывать противоположную тенденцию по сравнению с графиком зависимости высоты от полудисперсии. Я не знаю, как будет выглядеть соотношение высоты и полудисперсии для данных, которые я разместил здесь, но я предполагаю, что это будет похоже на изображения выше. можно ли будет увидеть ваш сюжет?   -  person Arihant    schedule 18.03.2016
comment
Я также добавил еще одну диаграмму, чтобы проиллюстрировать, что я ищу и как должен выглядеть общий шаблон.   -  person Arihant    schedule 18.03.2016


Ответы (2)


Используя gstat, вот пример:

library(raster)
library(gstat)
demo(meuse, ask = FALSE, echo = FALSE)
set.seed(131) # make random numbers reproducible
# add some noise with .1 variance
meuse.grid$dist = meuse.grid$dist + rnorm(nrow(meuse.grid), sd=sqrt(.1))
r = raster(meuse.grid["dist"])
v = variogram(dist~1, as(r, "SpatialPixelsDataFrame"))

(f = fit.variogram(v, vgm("Sph")))
#   model      psill    range
# 1   Nug 0.09035948    0.000
# 2   Sph 0.06709838 1216.737

f$psill[2] # sill
# [1] 0.06709838

f$range[2] # range
# [1] 1216.737

f$psill[1] # nugget
# [1] 0.09035948

Подключите свой собственный растр для r, и он должен работать. Измените Sph, чтобы он соответствовал другой модели вариограммы, попробуйте plot(v,f), чтобы проверить график.

person Edzer Pebesma    schedule 25.03.2016
comment
Спасибо, Эдзер. По какой-то причине я получаю эту ошибку: Ошибка в fit.variogram(v, vgm(Sph)) : модель должна относиться к классу variogramModel (использовать vgm) в строке (f = fit.variogram(v, vgm(Sph ))) - person Arihant; 25.03.2016
comment
Может быть, обновить свой пакет gstat? Какую версию вы используете? - person Edzer Pebesma; 25.03.2016
comment
Спасибо, это помогло, но при использовании моего собственного растрового слоя из примера многоканального TIFF, прикрепленного к вопросу, я получаю еще одну ошибку. › r= s[[81]] › Класс r: канал растра: 81 (из 99 каналов) размеры: 49, 50, 2450 (nrow, ncol, ncell) разрешение: 5, 5 (x, y) экстент: 364284, 364534, 4305537, 4305782 (xmin, xmax, ymin, ymax) координат. исх. : +proj=tmerc +lat_0=0 +lon_0=-75 +k=0.9995999932289124 +x_0=500000 +y_0=0 +ellps=WGS84 +units=m +no_defs мин Макс) - person Arihant; 26.03.2016
comment
› r= as(r, SpatialPixelsDataFrame) › v = variogram(dist~1, r) Ошибка в model.frame.default(terms(formula), as(data, data.frame), na.action = na.fail) : объект не является матрицей - person Arihant; 26.03.2016
comment
в ваших демонстрационных данных, используемых в вашем примере, замените dist на old_gap_. names(r) подскажет, какое имя использовать. - person Edzer Pebesma; 27.03.2016

Это всего лишь предположение. Вот как я оцениваю полудисперсию

где n — количество слоев, среднее значение которых меньше общего среднего. m — это общее среднее значение по всем слоям. r — это среднее значение каждого слоя, значение которого ниже общего среднего значения.

s <- stack("old_gap_.tif")
m <- cellStats(mean(s), stat="mean", na.rm=T) # 0.5620522
r <- m[m < 0.5620522]
sem <- 1/53 * (0.5620522 - r)^2
plot(sem, r)
person Geo-sp    schedule 18.03.2016