Какой ВИД вариограммы я должен построить перед выполнением кригинга в R?

У меня есть CSV-файл с именем seoul3112, содержащий концентрации PM10. пожалуйста, загрузите.. Я попытался построить образец вариограммы и подобрать модель в теме.

library(sp)
    library(gstat)
    library(rgdal)
    seoul3112<-read.csv("seoul3112.csv",row.name=1)
    seoul3112<-na.omit(seoul3112)

#назначить CRS и перепроецировать

coordinates(seoul3112)=~LON+LAT
proj4string(seoul3112) =  "+proj=longlat +datum=WGS84" 
seoul3112<-spTransform(seoul3112, CRS("+proj=utm +north +zone=52 +datum=WGS84"))

#построить полувариограмму

g<-gstat(id="PM10",formula=PM10~LON+LAT, data=seoul3112)
seoul3112.var<-variogram(g, cutoff=70000, width=6000)
seoul3112.var
plot(seoul3112.var, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 3112")

#Подходит для модели

model.3112<- fit.variogram(seoul3112.var,vgm(700,"Gau",40000,400),fit.method = 2)
                       
plot(seoul3112.var,model=model.3112, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 3112")

После написания этого кода я получил вот такую ​​полувариограмму.

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

Поскольку я очень новичок в геостатистике, я не понимаю, подходит ли моя вышеприведенная вариограмма для моего набора данных. Потому что в типичной вариограмме значение семидисперсии становится горизонтальным на пороге. Но эта вариограмма идет вверх! Должен ли я внести некоторые исправления в свой код?

Другое дело, что на самом деле моей конечной целью является проведение интерполяции методом кригинга на моем наборе данных (seoul3112). Я не понимаю, достаточно ли для проведения кригинга этой пробной вариограммы, или мне следует построить дирекционную вариограмму или что-то еще? Может ли кто-нибудь объяснить это подробно?


person Orpheus    schedule 28.07.2015    source источник
comment
Поскольку ваш вопрос касается скорее выбора и интерпретации статистических методов, чем программирования, я считаю, что сайт stats.stackexchange.com является более подходящим. Я проголосовал за перенос.   -  person Henrik    schedule 28.07.2015


Ответы (1)


Если вы посмотрите на свою подогнанную модель, она будет иметь параметр sill (т. е. nugget + psill), но он выходит за рамки ваших образцов.

 sill = sum(model.3112$psill)

Возможно, у вас не было достаточного расстояния между парой точек, чтобы добраться от range до sill. Я не думаю, что это не проблема, пока вы используете эту вариограмму для прогнозирования в пределах пространственной области ваших данных + 60000 м (расстояние, для которого у вас есть поддержка данных). Я придерживаюсь этой позиции, потому что наиболее важной частью вариограммы является ее начало, а по объему данных соответствие хорошее.

Что-то, что нужно проверить, это то, сколько пар точек (np) у вас есть, поддерживающих bins, которые вы начертили (чем больше пар точек, тем лучше).

Другое предложение заключается в поиске анизотропии с использованием графика уровней.

seoul3112.var_map<-variogram(g, cutoff=70000, width=6000, map=TRUE)
plot(seoul3112.var_map, col.regions=terrain.colors(20))

График уровней (вариограммы направлений)

или просмотрев подгонку во многих направлениях, установив alpha и нарисовав подгонку модели. Вам нужно исследовать только 180 градусов, потому что он симметричен (вы можете видеть это на графике ниже, где 0 и 180 совпадают).

seoul3112.var_angles<-variogram(g, cutoff=70000, width=6000, alpha=seq(0,180,30))
plot(seoul3112.var_angles, model=model.3112, pch=16, ylim=c(0,3000))

Направленные вариограммы

Если есть различия по направлениям, вы можете смоделировать большую и малую оси анизотропии, используя fit.variogram с anis, определенными для модели vgm(). Пример:

vgm(700,"Gau",40000,400, anis=c(someangle, someratio))
person Jared Smith    schedule 28.07.2015
comment
Как я могу понять по 7 направленной вариограмме выше есть анизотропия или нет? Я знаю о геометрической и зональной анизотропии. Но могу ли я сказать из вышеприведенных цифр, что подоконник почти одинаков во всех направлениях? а так же не могу понять что в этих цифрах разный диапазон или нет! Не могли бы вы прокомментировать анизотропию на основе этого рисунка? - person Orpheus; 29.07.2015
comment
Что вы действительно ищете, так это то, в каком направлении range ближе к sill, и в каком направлении длиннее всего. Из графика уровня и дирекционных вариограмм кажется, что 30 градусов могут быть малой осью (имеет некоторое увеличение полудисперсии, потенциально до порога), что делает 120 градусов главной осью. Я рекомендую построить эти углы отдельно и посмотреть, как выглядят эти вариограммы. - person Jared Smith; 29.07.2015
comment
Спасибо. Но я видел пример, где они показывают числовое значение порога и диапазона в таблице для каждого направления. Благодаря чему я могу четко определить большую и малую оси. Ссылка: techmat.vgtu.lt/~art/proc/file/ BudrLi.pdf (стр. 365). Есть ли встроенная функция в R, с помощью которой я могу видеть числовое значение диапазона и порога для каждого направления? - person Orpheus; 29.07.2015
comment
Если бы sills был достигнут во всех направлениях, то да, вы могли бы создать таблицу, подобную той, на которую вы ссылаетесь, используя fit.variogram на vgm(), определенном для variogram() с интересующим alpha (с разумным tol.hor на alpha). Затем используйте что-то вроде sill = sum(model.3112$psill) и range = model.3112$range для их извлечения. - person Jared Smith; 29.07.2015
comment
Вы имели в виду это? seoul3112.var1<-variogram(g,width=8000,cutoff=80000, alpha=seq(0,180,45),tol.hor=15) // model<- fit.variogram(seoul3112.var1,vgm(800,"Sph",40000,200),fit.method = 2) // sill = sum(model$psill) // range = model$range // sill // range Пробовал. но не получил желаемого результата. (здесь объект g, который я создал в своем вопросе) - person Orpheus; 29.07.2015
comment
Я имел в виду сначала использовать что-то вроде seoul3112.var1<-variogram(g,width=8000,cutoff=80000, alpha=MajorAngle,tol.hor=15), чтобы создать вариограмму только для основного направления анизотропии. Затем повторите для меньшего угла (или любого другого угла) и сравните. - person Jared Smith; 29.07.2015