объект kriging из geoR
- Игорь Черниенко
- Активный участник
- Сообщения: 137
- Зарегистрирован: 28 мар 2009, 01:05
- Репутация: 11
- Откуда: Хабаровск, Южно-Сахалинск
объект kriging из geoR
Рад приветствовать
Кто-нибудь может подсказать. можно ли вытащить грид из объекта kriging пакета geoR, с тем чтобы потом записать его на диск методом write.asciigrid?
Кто-нибудь может подсказать. можно ли вытащить грид из объекта kriging пакета geoR, с тем чтобы потом записать его на диск методом write.asciigrid?
-
- Гуру
- Сообщения: 4056
- Зарегистрирован: 15 окт 2010, 08:33
- Репутация: 1054
- Ваше звание: программист
- Откуда: Казань
Re: объект kriging из geoR
а в чем проблема? Координаты сетки у вас есть (они указывались в lacations), значения у вас лежат в результате (predict), сделайте объект SpatialGridDataFrame, и пишите. Ккак сделать - ниже пример из help (немного переделанный):Игорь Черниенко писал(а):Рад приветствовать
Кто-нибудь может подсказать. можно ли вытащить грид из объекта kriging пакета geoR, с тем чтобы потом записать его на диск методом write.asciigrid?
# 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")
- Игорь Черниенко
- Активный участник
- Сообщения: 137
- Зарегистрирован: 28 мар 2009, 01:05
- Репутация: 11
- Откуда: Хабаровск, Южно-Сахалинск
Re: объект kriging из geoR
В принципе я так и пытаюсь делать. Тонкость в том, что при расчете я использую опцию borders. Грид получается обрезанным. С координатами расчетные значения можно соединить, используя функцию locations.inside, но грид в этом случае получается не прямоугольным. Его можно объединить с исходным гридом функцией overlay, но опять-таки при попытке записи выдается сообщение Asciigrid does not support grids with non-square cells. Хотя ячейки в объединенном гриде самые что ни наесть квадратные.
-
- Гуру
- Сообщения: 4056
- Зарегистрирован: 15 окт 2010, 08:33
- Репутация: 1054
- Ваше звание: программист
- Откуда: Казань
Re: объект kriging из geoR
"Не создавайте излишних сущностей" (с) ОккамИгорь Черниенко писал(а):В принципе я так и пытаюсь делать. Тонкость в том, что при расчете я использую опцию borders. Грид получается обрезанным. С координатами расчетные значения можно соединить, используя функцию locations.inside, но грид в этом случае получается не прямоугольным. Его можно объединить с исходным гридом функцией overlay, но опять-таки при попытке записи выдается сообщение Asciigrid does not support grids with non-square cells. Хотя ячейки в объединенном гриде самые что ни наесть квадратные.
1) Если расчетные координаты не подмножество узлов прямоугольной сетки с квадратными ячеками, то в ArcGIS это передать в принципе нельзя (зачем тогда на них считать).
2) Если расчетные координаты есть подмножество узлов прямоугольной сетки с квадратными ячеками, то кто мешает создать эту сетку (см. предыдущее письмо), заполнить ее (-999999), а потом в нужные места поместить результаты кригинга? Позиция нужного места для (i,j) есть (i-1)*nx+j, если считать индексы с 1, и запись идет по строкам, слева направо, снизу вверх (см. предыдущее письмо).
- Игорь Черниенко
- Активный участник
- Сообщения: 137
- Зарегистрирован: 28 мар 2009, 01:05
- Репутация: 11
- Откуда: Хабаровск, Южно-Сахалинск
Re: объект kriging из geoR
Но ведь именно так я и делаю! Задачу упростил, выполнил расчет на прямоугольный грид. Далее:gamm писал(а): 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 (разделитель - ";") может посмотрите, что у меня не так?
- Вложения
-
- dp_gr.csv
- (525.2 КБ) 568 скачиваний
-
- Гуру
- Сообщения: 4056
- Зарегистрирован: 15 окт 2010, 08:33
- Репутация: 1054
- Ваше звание: программист
- Откуда: Казань
Re: объект kriging из geoR
обычно "что не так" смотрят в программеИгорь Черниенко писал(а): Я приложил результат расчета в виде 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")
- Игорь Черниенко
- Активный участник
- Сообщения: 137
- Зарегистрирован: 28 мар 2009, 01:05
- Репутация: 11
- Откуда: Хабаровск, Южно-Сахалинск
Re: объект kriging из geoR
Спасибо!
Что характерно, ни за что бы не подумал на ошибку округления, хотя именно она очевидна - создавал грид методом pred_grid(dph.geo$borders,by=1000). Буду иметь ввиду.
Отдельно рукопожимаю за "вспомогательные примеры :о)|||
Что характерно, ни за что бы не подумал на ошибку округления, хотя именно она очевидна - создавал грид методом pred_grid(dph.geo$borders,by=1000). Буду иметь ввиду.
Отдельно рукопожимаю за "вспомогательные примеры :о)|||
-
- Гуру
- Сообщения: 4056
- Зарегистрирован: 15 окт 2010, 08:33
- Репутация: 1054
- Ваше звание: программист
- Откуда: Казань
Re: объект kriging из geoR
поэтому все лучше делать руками, самому; особенно мне не нравится стиль "black box" gstat. В геостатистике главное - моделирование вариограммы, для чего нужно убрать тренд и использовать робастные естиматоры, ни того ни другого в автомате нет. Если мне не изменяет мой склероз, там даже геометрической анизотропии нет, и дисперсия всегда константа.Игорь Черниенко писал(а):Спасибо!
Что характерно, ни за что бы не подумал на ошибку округления, хотя именно она очевидна - создавал грид методом pred_grid(dph.geo$borders,by=1000). Буду иметь ввиду.
Отдельно рукопожимаю за "вспомогательные примеры :о)|||
Вообще, гипотеза стационарности слишком сильная штука, да и "тренд и остатки" - классическая "курица и яйцо", не зря придумывали всякие робастные median polish, чтобы тренд убрать. Если данных хватает, можно попробовать вычислять параметры вариограммы как функцию от координат, но в R такого нет, разве в geoR может добавят когда.
Так что лучше оценить вид вариограммы, и использовать gamm(), в который можно добавить и вариограмму, и модель дисперсии. Хотя возиться приходится достаточно долго, особенно если одновременно строится модель и оцениваются параметры вариограммы, особенно на больших данных. Зато есть оценки значимости, возможность сравнивать модели, и пр. "вкусности".
- Игорь Черниенко
- Активный участник
- Сообщения: 137
- Зарегистрирован: 28 мар 2009, 01:05
- Репутация: 11
- Откуда: Хабаровск, Южно-Сахалинск
Re: объект kriging из geoR
О, оказывается, "не прямоугольный" грид тоже пишется на диск :о)
- Игорь Черниенко
- Активный участник
- Сообщения: 137
- Зарегистрирован: 28 мар 2009, 01:05
- Репутация: 11
- Откуда: Хабаровск, Южно-Сахалинск
Re: объект kriging из geoR
Кстати, что такое gamm(), который Вы советуете использовать? В сети находятся только сссылки на Ваши посты на ГисЛабе :о)|||
-
- Гуру
- Сообщения: 4056
- Зарегистрирован: 15 окт 2010, 08:33
- Репутация: 1054
- Ваше звание: программист
- Откуда: Казань
Re: объект kriging из geoR
не в той сети ищите, я русские сайты вообще отключил в поискеИгорь Черниенко писал(а):Кстати, что такое 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).
Если перечисленные выше слова не знакомы, ее лучше почитать. И вообще, книжка - одна из лучших, что я видел.
Кто сейчас на конференции
Сейчас этот форум просматривают: нет зарегистрированных пользователей и 25 гостей