Подложите векторное изображение к сетке для кригинга в R

Question

После долгих поисков, расспросов и написания некоторого кода я как бы получил минимум для выполнения кригинга в gstat R.

Используя 4 очка (я знаю, совершенно плохо), я кригал несэмплированные точки, расположенные между ними. Но на самом деле мне не нужны все эти моменты. Внутри этой области есть меньший подрайон... эта область - то, что мне действительно нужно.

Короче говоря... У меня есть измерения, взятые с 4 метеостанций, которые сообщают данные об осадках. Координаты широты и долготы для этих точек:

lat    long
7.16   124.21
8.6    123.35
8.43   124.28
8.15   125.08

Мой путь к кригингу можно увидеть через мои предыдущие вопросы о Stackru.

Это: Создать вариограмму в пакете R Gstat

И это: создать Grid в 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")

РЕДАКТИРОВАТЬ: код изменился с прошлого раза. Но это не имеет значения, я думаю, я просто хочу знать, возможно ли это, и кто-нибудь может привести самый простой пример того, как это возможно. Я могу найти решение для примера моей собственной проблемы оттуда.

Answer 1

Дайте мне знать, если вы найдете это полезным:

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

Для этого вы загружаете свои векторные данные:

donut <- rgdal::readOGR('/variogram', 'donut')
day_1@proj4string@projargs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" # Becouse donut shape have this CRS

plot(donut, axes = TRUE, col = 3)
plot(day_1, col = 2, pch = 20, add = TRUE)

Затем вы удаляете "внешнее кольцо" и сохраняете инсайдера. Также указывает, что второе больше не дыра:

hole <- donut # for keep original shape
hole@polygons[1][[1]]@Polygons[1] <- NULL
hole@polygons[1][[1]]@Polygons[1][[1]]@hole <- FALSE

plot(hole, axes = TRUE, col = 4, add = TRUE)

После этого вы проверяете, какие точки находятся внутри "дырочного" нового синего векторного слоя:

over.pts <- over(day_1, hole)
day_1_subset <- day_1[!is.na(over.pts$Id), ]

plot(donut, axes = TRUE, col = 3)
plot(hole, col = 4, add = TRUE)
plot(day_1, col = 2, pch = 20, add = TRUE)
plot(day_1_subset, col = 'white', pch = 1, cex = 2, add = TRUE)

write.csv(day_1_subset@data, 'myfile.csv') # write intersected points table
write.csv(as.data.frame(coordinates(day_1_subset)), 'myfile.csv') # write intersected points coords
writeOGR(day_1_subset, 'path', 'mysubsetlayer', driver = 'ESRI Shapefile') # write intersected points shape

С помощью этого кода вы можете заполнить кольцевую или кольцевую дыру, если у вас уже есть шейп-файл. Если у вас есть изображение и вы хотите его обрезать, попробуйте следующее:

В случае, если вы загружаете растр (получаете изображение базовой карты из Интернета):

coordDf <- as.data.frame(coordinates(day_1)) # get basemap from points
# coordDf <- data.frame(hole@polygons[1][[1]]@Polygons[1][[1]]@coords) # get basemap from hole
colnames(coordDf) <- c('x', 'y') 
imag <- dismo::gmap(coordDf, lonlat = TRUE)
myimag <- raster::crop(day_1.kriged, hole)
plot(myimag)
plot(day_1, add = TRUE, col = 2)

Если вы используете day_1.kriged:

myCropKrig<- raster::crop(day_1.kriged, hole)

  myCropKrig %>% 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) +
  geom_point(data=coordDf[!is.na(over.pts$Id), ], aes(x=x, y=y), color="blue", size=3, shape=20) +
  theme_bw()

И "Наконец, я беру значения области, покрытой изображением, и сохраняю их в формате CSV или в базе данных".

write.csv(as.data.frame(myCropKrig), 'myCropKrig.csv')

Надеюсь, вы найдете это полезным, и я отвечу на ваш смысл

Answer 2

Чтобы упростить ваш вопрос:

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

Требуется несколько шагов

1. Вы должны использовать QGIS для привязки вашего изображения (Raster > Georeferencer). Вам нужна карта с географической привязкой в ​​фоновом режиме, чтобы помочь. Это создает растровый объект с пространственной информацией.
2. Две возможности
   • 2.а Центральная часть вашего изображения имеет цвет, который можно напрямую использовать в качестве маски в R (например, все зеленые пиксели в середине красных пикселей).
   • 2.b. Если нет, вам нужно использовать QGIS, чтобы вручную очертить полигон области (Layer > Create Layer > New Shapefile > Polygon)
3. Импортируйте свой растровый или полигональный шейп-файл в R
4. Использовать функцию raster::mask чтобы извлечь значения вашей интерполяции, используя растровое изображение или SpatialPolygon.