У меня есть следующий фрейм данных с тремя столбцами, а именно «Широта», «Долгота» и «значение», а также растровый файл (прикрепленный с именем «ras.tif ') и файл формы области (как прикрепленный с именем "okl_shape.shp"):
Dataframe
df <- data.frame(Latitude=c(34.80833, 34.79851, 34.91418, 36.07204, 36.80253, 35.96305,
36.69256, 36.83129, 35.5915, 34.8497, 34.60896, 35.65282, 36.74813, 35.54615, 33.92075,
35.54848, 35.20494, 36.26353, 36.84053, 36.60183), Longitude=c(-98.02325, -96.66909, -98.29216,
-99.90308, -100.53012, -95.86621, -102.49713, -99.64101, -99.27059, -97.0033, -96.33309,
-96.80407, -98.36274, -99.7279, -96.32027, -98.03654, -99.80344, -98.49766, -96.42777, -101.6013
), value=c(0.84629845, NA, 0.916147287, 0.735707364, 0.432443798, 0.959728682, 0.123924419,
0.849589147, 0.307998062, 0.932215116, 0.939897287, 0.915825581, NA, 0.273945736, NA,
0.705023256, 0.447494186, 0.840686047, 0.901098837, 0.202523256))
Прикрепленные файлы
Вот файл формы моей области: файл Вот растровый файл: файл
Описание проблемы
Мне нужно создать график, сравнивающий точечные данные (т.е. фрейм данных, описанный выше) с прикрепленным растром (например, "ras.tif" ), и файл формы должен перекрывать эти графики. См. рисунок ниже для справки (на этом рисунке левая панель отображает значения точек, как определено во фрейме данных, правая панель отображает растр, а многоугольник представляет файл формы области): введите описание изображения здесь
Мой код
Я попытался решить это с помощью "spplot", но не смог получить желаемый график (что-то похожее на график выше). Кроме того, я не знаю, как мне вставить файл формы в свой график. Ниже приведен код, который я пробовал до сих пор:
library(raster)
library(sf)
library(maptools)
meso.shape <- st_read('okl_shape.shp')
ras1 <- raster('ras.tif')
df <- data.frame(Latitude=c(34.80833, 34.79851, 34.91418, 36.07204, 36.80253, 35.96305,
36.69256, 36.83129, 35.5915, 34.8497, 34.60896, 35.65282, 36.74813, 35.54615, 33.92075,
35.54848, 35.20494, 36.26353, 36.84053, 36.60183
), Longitude=c(-98.02325, -96.66909, -98.29216, -99.90308, -100.53012, -95.86621, -102.49713,
-99.64101, -99.27059, -97.0033, -96.33309, -96.80407, -98.36274, -99.7279, -96.32027, -98.03654,
-99.80344, -98.49766, -96.42777, -101.6013), value=c(0.84629845, NA, 0.916147287, 0.735707364,
0.432443798, 0.959728682, 0.123924419, 0.849589147, 0.307998062, 0.932215116, 0.939897287,
0.915825581, NA, 0.273945736, NA, 0.705023256, 0.447494186, 0.840686047, 0.901098837,
0.202523256))
coordinates(df) <- ~Longitude+Latitude
x <- vect2rast(df, fname = names(df)[1])
x <- raster(x)
extent(x) <- extent(ras1)
crs(x) <- crs(ras1)
x <- projectRaster(x,ras1)
s <- stack(ras1,x)
sp <- as(s, 'SpatialPolygons')
p <- spplot(s, names.attr = c("Simulated FWI","Observed FWI"))
show(p)