объект kriging из geoR

[Игорь Черниенко]

Рад приветствовать
Кто-нибудь может подсказать. можно ли вытащить грид из объекта kriging пакета geoR, с тем чтобы потом записать его на диск методом write.asciigrid?

[gamm]

Рад приветствовать
Кто-нибудь может подсказать. можно ли вытащить грид из объекта kriging пакета geoR, с тем чтобы потом записать его на диск методом write.asciigrid?

а в чем проблема? Координаты сетки у вас есть (они указывались в lacations), значения у вас лежат в результате (predict), сделайте объект SpatialGridDataFrame, и пишите. Ккак сделать - ниже пример из help (немного переделанный):
# toy example:
df = data.frame(z = rnorm(25), # значения (по строкам, в соответствии с моими координатами) = predict
xc = rep(1:5,5), # центры ячеек по x = locations[,1]
yc = rep(1:5,each=5)) # центры по y = locations[,2]
coordinates(df) = ~xc+yc
gridded(df) = TRUE
df = as(df, "SpatialGridDataFrame") # to full grid

write.asciigrid(df, "foo.asc")

[Игорь Черниенко]

В принципе я так и пытаюсь делать. Тонкость в том, что при расчете я использую опцию borders. Грид получается обрезанным. С координатами расчетные значения можно соединить, используя функцию locations.inside, но грид в этом случае получается не прямоугольным. Его можно объединить с исходным гридом функцией overlay, но опять-таки при попытке записи выдается сообщение Asciigrid does not support grids with non-square cells. Хотя ячейки в объединенном гриде самые что ни наесть квадратные.

[gamm]

В принципе я так и пытаюсь делать. Тонкость в том, что при расчете я использую опцию borders. Грид получается обрезанным. С координатами расчетные значения можно соединить, используя функцию locations.inside, но грид в этом случае получается не прямоугольным. Его можно объединить с исходным гридом функцией overlay, но опять-таки при попытке записи выдается сообщение Asciigrid does not support grids with non-square cells. Хотя ячейки в объединенном гриде самые что ни наесть квадратные.

"Не создавайте излишних сущностей" (с) Оккам

1) Если расчетные координаты не подмножество узлов прямоугольной сетки с квадратными ячеками, то в ArcGIS это передать в принципе нельзя (зачем тогда на них считать).

2) Если расчетные координаты есть подмножество узлов прямоугольной сетки с квадратными ячеками, то кто мешает создать эту сетку (см. предыдущее письмо), заполнить ее (-999999), а потом в нужные места поместить результаты кригинга? Позиция нужного места для (i,j) есть (i-1)*nx+j, если считать индексы с 1, и запись идет по строкам, слева направо, снизу вверх (см. предыдущее письмо).

[Игорь Черниенко]

2) Если расчетные координаты есть подмножество узлов прямоугольной сетки с квадратными ячеками, то кто мешает создать эту сетку (см. предыдущее письмо), заполнить ее (-999999), а потом в нужные места поместить результаты кригинга? Позиция нужного места для (i,j) есть (i-1)*nx+j, если считать индексы с 1, и запись идет по строкам, слева направо, снизу вверх (см. предыдущее письмо).

Но ведь именно так я и делаю! Задачу упростил, выполнил расчет на прямоугольный грид. Далее:
> dp.gr<-cbind(cmm.gr2,dp.kg$pred)#cmm.gr2 - прямоугольный грид,dp.kg-результат расчета
> names(dp.gr)<-list('lon','lat','dp')
> coordinates(dp.gr)=c('lon','lat')
> gridded(dp.gr)=T
> dp.gr<-as(dp.gr,'SpatialGridDataFrame')
> write.asciigrid(dp.gr,'dp_gr.asc')
Однако на диск не пишется даже результат расчета на прямоугольную сетку.

Я приложил результат расчета в виде csv (разделитель - ";") может посмотрите, что у меня не так?

[gamm]

Я приложил результат расчета в виде csv (разделитель - ";") может посмотрите, что у меня не так?

обычно "что не так" смотрят в программе

в данном случае - ошибка округления - можно посмотреть, где это выводится, набрав
write.asciigrid
или в отладчике, набрав debug(write.asciigrid) и запустив запись.
Заодно показал, как записывать произвольное подмножество значений сетки

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

p<-read.table("dp_gr.csv",header=TRUE,sep=";")
names(p)
# [1] "lon" "lat" "dp"

# небольшой контроль ...
p.x<-sort(unique(p$lon)); p.nx<-length(p.x)
p.y<-sort(unique(p$lat)); p.ny<-length(p.y)

print(c(p.nx*p.ny,nrow(p)))

print(p.x[-1] - p.x[-p.nx])
print(p.y[-1] - p.y[-p.ny])

# укладываем значения в нужные места, независимо от того, вся сетка заполнена, или нет 
pos.x <-match(p$lon,p.x)
pos.y <-match(p$lat,p.y)
pos.xy<-(pos.y-1)*p.nx + pos.x
pos.z <-1:(p.nx*p.ny)
pos.xyz<-match(pos.xy,pos.z)

p.dp = data.frame(z = rep(NA,p.nx*p.ny),
                  xc = rep(p.x,p.ny),       
                  yc = rep(p.y,each=p.nx))  
p.dp$z[pos.xyz]<-p$dp

# превращаем в сетку
coordinates(p.dp) = ~xc+yc
gridded(p.dp) = TRUE
p.dp = as(p.dp, "SpatialGridDataFrame") # to full grid
print(p.dp@grid@cellsize[1]-p.dp@grid@cellsize[2])

# ошибка округления - в программе проверка на равно, поэтому делаем сетку строго квадратной
print(p.cellsize<-0.5*(p.dp@grid@cellsize[1]+p.dp@grid@cellsize[2]))
p.dp@grid@cellsize[1:2]<-p.cellsize

write.asciigrid(p.dp, "foo.asc")

[Игорь Черниенко]

Спасибо!
Что характерно, ни за что бы не подумал на ошибку округления, хотя именно она очевидна - создавал грид методом pred_grid(dph.geo$borders,by=1000). Буду иметь ввиду.
Отдельно рукопожимаю за "вспомогательные примеры :о)|||

[gamm]

Спасибо!
Что характерно, ни за что бы не подумал на ошибку округления, хотя именно она очевидна - создавал грид методом pred_grid(dph.geo$borders,by=1000). Буду иметь ввиду.
Отдельно рукопожимаю за "вспомогательные примеры :о)|||

поэтому все лучше делать руками, самому; особенно мне не нравится стиль "black box" gstat. В геостатистике главное - моделирование вариограммы, для чего нужно убрать тренд и использовать робастные естиматоры, ни того ни другого в автомате нет. Если мне не изменяет мой склероз, там даже геометрической анизотропии нет, и дисперсия всегда константа.

Вообще, гипотеза стационарности слишком сильная штука, да и "тренд и остатки" - классическая "курица и яйцо", не зря придумывали всякие робастные median polish, чтобы тренд убрать. Если данных хватает, можно попробовать вычислять параметры вариограммы как функцию от координат, но в R такого нет, разве в geoR может добавят когда.

Так что лучше оценить вид вариограммы, и использовать gamm(), в который можно добавить и вариограмму, и модель дисперсии. Хотя возиться приходится достаточно долго, особенно если одновременно строится модель и оцениваются параметры вариограммы, особенно на больших данных. Зато есть оценки значимости, возможность сравнивать модели, и пр. "вкусности".

[Игорь Черниенко]

О, оказывается, "не прямоугольный" грид тоже пишется на диск :о)

[Игорь Черниенко]

Кстати, что такое gamm(), который Вы советуете использовать? В сети находятся только сссылки на Ваши посты на ГисЛабе :о)|||

[gamm]

Кстати, что такое gamm(), который Вы советуете использовать? В сети находятся только сссылки на Ваши посты на ГисЛабе :о)|||

не в той сети ищите, я русские сайты вообще отключил в поиске

GAM/GAMM - Genaralized Additive (Mixed) Model, обобщенная аддитивная модель (со случайными эффектами)

Обобщенная означает, что распределение может отличаться от Гауссовского (логистическая модель, Пуассоновское, Гамма, и т.д. распределения), аддитивная означает, что линейный предиктор обобщенной модели есть сумма гладких функций. Со случайными эффектами означает, что дисперсия моделируется не только в ошибке линейного предиктора, но и в других местах, через случайные эффекты моделируется сглаживание. В российской транскрипции это модели со смешанными эффектами, но их нигде, кроме Statistica, я их не встречал, (но более безобразной документации, чем у Statistica, я еще не видел, даже всякие реферат.ру - кладезь разума на ее фоне).

Есть очень хорошая книжка, Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press. Я ее купил. Wood - автор пакета mgcv, который есть в R (даже в стандартной поставке, вместе с ядром, т.е. recommended).

Если перечисленные выше слова не знакомы, ее лучше почитать. И вообще, книжка - одна из лучших, что я видел.