После долгих поисков, расспросов и написания кода я как бы получил необходимый минимум для кригинга в gstat R.
Используя 4 точки (знаю, совсем плохо), я провел кригинг несэмплированных точек, расположенных между ними. Но на самом деле мне не нужны все эти точки. Внутри этой области есть меньшая подобласть... это та область, которая мне действительно нужна.
Короче говоря. У меня есть измерения, сделанные с 4 метеостанций, которые сообщают данные об осадках. Широта и долгота для этих точек:
lat long
7.16 124.21
8.6 123.35
8.43 124.28
8.15 125.08
Мой путь к кригингу можно увидеть из моих предыдущих вопросов на StackOverflow.
Это: Создайте вариограмму в пакете gstat R
А это: Создать сетку в R для кригинга в gstat
Я знаю, что изображение имеет координаты (по крайней мере, по моим оценкам):
Leftmost: 124 13ish 0 E(DMS)
Rightmost : 124 20ish 0 E
Topmost corrdinates: 124 17ish 0 E
Bottommost coordinates: 124 16ish 0 E
Для этого будет иметь место преобразование, но это не имеет значения, я думаю, или с этим будет проще разобраться позже.
Изображение также нерегулярное (но не все).
Думайте об этом как о пончике: вы криируете всю круглую форму пончика, но вам нужна только область, покрытая отверстием, поэтому вы удаляете или, по крайней мере, игнорируете значения, полученные от самого пончика.
У меня есть изображение (.jpg) рассматриваемой области, мне нужно будет преобразовать изображение в шейп-файл или какой-либо другой векторный формат с помощью QGIS или аналогичного программного обеспечения. После этого мне нужно будет вставить это векторное изображение внутрь области 4-точечного кригинга, чтобы я знал, какие координаты учитывать, а какие удалить.
Наконец, я беру значения области, покрытой изображением, и сохраняю их в CSV или базе данных.
Кто-нибудь знает, как я могу начать с этого? Полный нуб в R и статистике. Спасибо всем, кто ответит.
Я просто хочу знать, возможно ли это, и если это дать несколько советов. Спасибо еще раз.
Можно также опубликовать мой сценарий:
suppressPackageStartupMessages({
library(sp)
library(gstat)
library(RPostgreSQL)
library(dplyr) # for "glimpse"
library(ggplot2)
library(scales) # for "comma"
library(magrittr)
library(gridExtra)
library(rgdal)
library(raster)
library(leaflet)
library(mapview)
})
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, dbname="Rainfall Data", host="localhost", port=5432,
user="postgres", password="postgres")
day_1 <- dbGetQuery(con, "SELECT lat, long, rainfall FROM cotabato.sample")
coordinates(day_1) <- ~ lat + long
plot(day_1)
x.range <- as.integer(c(7.0,9.0))
y.range <- as.integer(c(123.0,126.0))
grid <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=0.05),
y=seq(from=y.range[1], to=y.range[2], by=0.05))
coordinates(grid) <- ~x+y
plot(grid, cex=1.5)
points(day_1, col='red')
title("Interpolation Grid and Sample Points")
day_1.vgm <- variogram(rainfall~1, day_1, width = 0.02, cutoff = 1.8)
day_1.fit <- fit.variogram(day_1.vgm, model=vgm("Sph", psill = 8000, range = 1))
plot(day_1.vgm, day_1.fit)
plot1 <- day_1 %>% as.data.frame %>%
ggplot(aes(lat, long)) + geom_point(size=1) + coord_equal() +
ggtitle("Points with measurements")
plot(plot1)
############################
plot2 <- grid %>% as.data.frame %>%
ggplot(aes(x, y)) + geom_point(size=1) + coord_equal() +
ggtitle("Points at which to estimate")
plot(plot2)
grid.arrange(plot1, plot2, ncol = 2)
coordinates(grid) <- ~ x + y
############################
day_1.kriged <- krige(rainfall~1, day_1, grid, model=day_1.fit)
day_1.kriged %>% as.data.frame %>%
ggplot(aes(x=x, y=y)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
scale_fill_gradient(low = "yellow", high="red") +
scale_x_continuous(labels=comma) + scale_y_continuous(labels=comma) +
theme_bw()
write.csv(day_1.kriged, file = "Day_1.csv")
РЕДАКТИРОВАТЬ: код изменился с прошлого раза. Но это не имеет значения, я думаю, я просто хочу знать, возможно ли это, и может ли кто-нибудь привести простейший пример того, как это возможно. Отсюда я могу вывести решение примера для своей проблемы.
