Как создать вывод с координатной сеткой в ​​​​R и исключить квадраты сетки, которые не находятся над сушей?

Я пытаюсь получить данные об осадках с координатной сеткой над Великобританией, используя алгоритм Thin Plate Spline, и исключить значения, которые не относятся к суше, в R - процесс, который я пока могу выполнить только вручную. Проблема сложная (для меня) и даже сложная для объяснения, поэтому я пройдусь по тому, что я сделал до сих пор. Любая помощь будет принята с благодарностью.

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

library(fields)
library(maptools)
library(gstat)

dat <- read.table("1961month1day1.csv", header=T, sep=",", quote = "")
names(dat) <- c("easting", "northing", "dailyrainfall","avaerageyearlyrainfall")

Вот пример данных:

dput(head(dat, 20))
structure(list(easting = c(130000L, 145000L, 155000L, 170000L, 
180000L, 180000L, 180000L, 180000L, 185000L, 200000L, 200000L, 
205000L, 210000L, 220000L, 225000L, 230000L, 230000L, 230000L, 
230000L, 235000L), northing = c(660000L, 30000L, 735000L, 40000L, 
30000L, 45000L, 60000L, 750000L, 725000L, 50000L, 845000L, 65000L, 
770000L, 105000L, 670000L, 100000L, 620000L, 680000L, 95000L, 
120000L), dailyrainfall = c(9.4, 4.1, 12.4, 2.8, 1.3, 3.6, 4.8, 26.7, 19.8, 
4.6, 1.7, 4.1, 12.7, 1.8, 3, 5.3, 1, 1.5, 1.5, 4.6), averageyearlyrainfall = c(1334.626923, 
1123.051923, 2072.030769, 1207.584615, 928, 1089.334615, 880.0884615, 
2810.323077, 1933.719231, 1215.642308, 2644.171154, 1235.913462, 
2140.111538, 1010.436538, 1778.432692, 1116.934615, 912.2807692, 
1579.386538, 1085.498077, 1250.601923)), .Names = c("easting", 
"northing", "dailyrainfall", "averageyearlyrainfall"), row.names = c(NA, 20L), class = "data.frame")

Затем я могу подогнать к данным сплайн тонкой пластины, чтобы получить поверхность с сеткой, и построить поверхность:

fit <- Tps(cbind(dat$easting,dat$northing),dat$dailyrainfall)
surface(fit)

Затем я могу создать сетку Великобритании с шагом в 1 км, используя:

xvals <- seq(0, 700000, by=1000)
yvals <- seq(0, 1250000, by=1000)

а затем нанесите поверхность на эту сетку и запишите данные в таблицу:

griddf <- expand.grid(xvals, yvals)
griddf$pred <- predict(fit, x=as.matrix(griddf))
write.table(griddf, file="1Jan1961grid.csv", sep=",", qmethod="double")

Отлично - пока все хорошо. Теперь я преобразовал свои точечные данные в данные с координатной сеткой 1 км по всей сетке от 0 до 700000 (E) и от 0 до 1250000 (N). Таблица письменных данных представляет собой список, содержащий индекс, восточное и северное положение и прогнозируемое количество осадков.

Теперь задача — я хочу исключить из этого списка все значения, которые не относятся к суше. Я могу добиться этого вручную, загрузив данные в Excel (или Access) и сравнив данные с другим файлом, который содержит ту же сетку и среднегодовое количество осадков (файл называется 1kmgridaveragerainfall.csv). Вот пример этого файла:

dput(head(dat1, 20))
structure(list(easting = c(-200000L, -200000L, -200000L, -200000L, 
-200000L, -200000L, -200000L, -200000L, -200000L, -200000L, -200000L, 
-200000L, -200000L, -200000L, -200000L, -200000L, -200000L, -200000L, 
-200000L, -200000L), northing = c(1245000L, 1240000L, 1235000L, 
 1230000L, 1225000L, 1220000L, 1215000L, 1210000L, 1205000L, 1200000L, 
 1195000L, 1190000L, 1185000L, 1180000L, 1175000L, 1170000L, 1165000L, 
 1160000L, 1155000L, 1150000L), averageyearlyrainfall = c(-9999, -9999, -9999, 
 -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, 
 -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999)), .Names = c("easting", 
 "northing", "averageyearlyrainfall"), row.names = c(NA, 20L), class = "data.frame")

В любом квадрате сетки, который не находится над сушей, среднегодовое количество осадков составляет -9999. Следовательно, после сопоставления (т. е. с использованием vlookup или запроса в Access) я могу отфильтровать значения, которые имеют это значение -9999, и это оставляет мне таблицу данных, в которой есть восточное и северное направление, дневные осадки и среднегодовые осадки только для значений земли. Затем я могу загрузить это обратно в R и построить это, используя:

quilt.plot(cbind(dat$easting,dat$northing),dat$mm, add.legend=TRUE, nx=654, ny=1209,xlim=c(0,700000),ylim=c(0,1200000))

и у меня остался участок осадков над территорией Великобритании (а не над морем).

Итак, может ли кто-нибудь предложить способ добиться того же, но без всей фильтрации и т. Д., Используя Excel или доступ, т. Е. Можно ли добиться того же, используя только R? Есть ли способ загрузить обе таблицы данных в R в начале и каким-то образом подогнать TPS точечных данных к средним данным, чтобы квадраты сетки, равные -9999, не отображались.

Я знаю, что TPS можно взвешивать с помощью ковариаты (Z) — это вообще помогает? то есть

fit <- Tps(cbind(dat$easting,dat$northing),dat$dailyrainfall, Z=dat$averageyearlyrainfall)

Кроме того, когда я выполняю поверхность (подгонку) исходного TPS, как мне расширить график до краев графика - я уверен, что читал это где-то, где вы вводите что-то вроде interp=TRUE, но это не так. работай.

Любая помощь будет принята с благодарностью

Спасибо, Тони