У меня есть база данных (netcdf) вихрей, движущихся через океан, и множество следов рыбы, делающих то же самое. Я пытаюсь добавить информацию о вихрях к следу рыбы, когда рыба находится в пределах досягаемости вихря (местоположения в пределах (изменяющегося) порогового расстояния и даты идентичны). Структура данных, сложности и воспроизводимый пример приведены ниже. Я подозреваю, что это не супер сложно, я просто завязываю себя в узлы, пытаясь найти правильное решение. Заранее спасибо за любые мысли!
Структура данных / воспроизводимый пример: 3 дня рыбы и вихревой дорожки
library(lubridate)
fish <- data.frame(lat = c(42.1, 42.6, 43.2), lon = c(-10, -10.1, -10.2), date = ymd(c("1990-01-01", "1990-01-02", "1990-01-02")))
eddy <- data.frame(lat = c(44, 42.3, 40), lon = c(-15, -10.1, -6), date = ymd(c("1990-01-01", "1990-01-02", "1990-01-02")), track = c(1,1,1), radius = c(81,82,83), vals = c(11,12,13))
Вот тут немного сложнее:
- круговой вихревой буфер с переменным радиусом
Размер вихря (из которого рассчитывается перекрытие или расстояние) задается в столбце радиуса, поэтому:
sf_eddy <- st_as_sf(eddy)
st_buffer(sf_eddy, dist = radius)
Для st_buffer «Все операции работают для каждого объекта», поэтому я не уверен, будет ли это рассматривать sf_eddy как один объект (возможно, (много) линейная строка) или ряд отдельных точек. Я полагаю, что могу поиграть с этим, пока не получу его на работу, так что каждый вихревой день (n = 151937) представляет собой буферный круг sf, к которому я могу предположительно добавить его vals, если они еще не включены.
- Заданный на расстоянии латунный прямоугольный буфер рыбы
Фиксированное расстояние ошибки для положений рыбы составляет + -1,9 ° широты, + -0,77 ° долготы. Так что концептуально я хочу создать высокий тонкий многоугольник вокруг каждого дня рыбы. st_buffer не позволяет этого, но могут быть и другие варианты , возможно, путем создания многоугольника вокруг каждого центроида согласно square примеру здесь .
- Если рыба и любой вихрь перекрываются в пространстве и времени, добавьте вихревые значения к рыбе
Я мог бы потенциально сократить рабочую нагрузку, выполнив только вышеуказанные шаги, когда eddy $ date == fish $ date? Хотя расчет вихрей будет разовым, и я могу сохранить результаты. Я подозреваю, что расчет рыбных ящиков не будет особенно обременительным. В любом случае:
fish$vals <- rep(NA, nrow(fish))
sf_fish <- sf_as_sf(fish)
for (i in seq(along = sf_fish)){
if(st_intersects(sf_fish[i,],sf_eddy)) sf_fish[i,"vals"] <- sf_eddy$vals}
Это, вероятно, ужасный способ сделать это, очевидно.
Но является ли какой-либо из этих способов хоть немного разумным способом решения этой проблемы, или я упускаю гораздо более изящные решения из-за своего любительского знания sf, sp, rgeos и т. Д.? У меня 166 рыбных файлов со средним числом мест 286 дней = всего 47476 поисков. Было бы хорошо / быстро создать индекс пар совпадений дат рыба / вихрь (может быть несколько вихрей в один и тот же день), а затем проверять только пересечения на них? Поскольку мне нужно добавить eddy$vals (на самом деле 2 или 3 столбца) к fish, было бы лучше использовать пространственное соединение, возможно, st_join?
Редактировать: комментарии и выходные странности карты:
600 000 выводов ggplot рыбы / вихревой строки:
Похоже, мне нужно исправить проекцию, так как вихри не печатаются к западу от 0. Теперь они точно в том же CRS; возможно, оригиналы 0-360, а не -180: 180. Кроме того, в то время как коробки с рыбой печатаются хорошо, вихри слишком малы: имеют значение r = 80 км = ~ 1 градус @ 45N, и они определенно печатают намного меньше, чем это. Я тоже подумаю об этом. Также: возникнет ли проблема для for (), если в один день существуют 2 вихря?
Edit2: оригиналы были 0: 360, которые я исправил. Радиусы вихрей все еще слишком малы:
Согласно обсуждению кода, оригинальный вихревой файл с латлонами был преобразован с st_as_sf в crs4326 (wgs84), а затем st_transform(sf_eddy,6931), который является Северной Атлантикой, согласно изображению выше. Столбец геометрии в результирующем файле имеет размер POLYGON, размер XY, значения в проецируемых координатах, а не в латлонах, например, sf_eddy_buffered$geometry[1] bbox: xmin -4894790 ymin -2418555 xmax -4894648 ymax -2418413. speed_radius равно 71, поэтому разница между xmin и xmax составляет 142, что составляет 71 * 2, что является правильным. Разве единица радиуса не равна единице проекционного расстояния? Примерами Стюарта были 81000. Единица радиуса вихря: км. Единица проекции / геометрия: https://epsg.io/6931 говорит, что единицы измерения действительно являются метрами. Умножьте радиус на 1000 и попробуйте снова:
Такэто все хорошо. Только что попробовал обновленный код st_join, думаю, что он может использовать некоторую работу - возможно, это будет длиться вечно, вероятно, из-за 3639209 рядов вихрей против 461 ряда рыб. Также он создал несколько строк для нескольких вихрей для одной и той же рыбы в один и тот же день. Потенциально функция относительно большой возможности перекрытия, учитывая буфер ошибок вокруг рыбы.
Я подозреваю, что смогу исправить это, отметив дубликаты и удалив дубликаты, которые находятся дальше всего. Текущий способ (через v и rbind) лучше, чем предыдущий (прямое присвоение ячейки значения) IMO, так как в противном случае несколько значений будут отправлены в одну ячейку, либо сбой, либо перезапись без вывода сообщений.
Для ускорения кода я подумал, что я мог бы установить для файла вихря только совпадение в диапазоне дат текущей рыбы и использовать это подмножество в цикле соединения. Цикл соединения уже подразделяет файлы fish и eddy по дате в любом случае, так что, может быть, это избыточно, или, возможно, уменьшение размера вихревых данных, над которыми нужно работать, будет каждый раз уменьшать нагрузку на CPU / RAM?
Edit3: это УДОВОЛЬСТВЕННО ускоряет вещи. Добавление простого счетчика циклов для thisDate показало, что оно завершено на 275 из 461 строки. И создал файл длиной 481. Возможно, есть много дней, когда рыба не была вихрем (цикл thisDate приводит к отсутствию v и rbind), и много дней, когда рыба находилась в 2+ вихрях (* 1077) * цикл приводит к многорядным v и rbind). По-прежнему ощущается, что он должен обработать 461 цикл, однако ...
Edit4: небольшая настройка:
v <- st_join(thisFishRow, thisEddyRow, left = F, largest = TRUE)
largest означает, что он присоединится к вихрю только с наибольшей площадью перекрытия, в результате чего максимальная длина строки v будет равна 1, и, следовательно, 1 значение в день будет максимально. Тем не менее, нюансы путей означают, что рыба может двигаться вперед и назад между различными продолжающимися вихрями. Мое подозрение (основанное на аналогичном исследовании) состоит в том, что они, скорее всего, останутся в сильной антициклонической вихревой системе. Радиус (и, следовательно, площадь) вихревых элементов управления перекрывают, возможно, больше, чем положение. Для полуторамесячного раздела я нашел наибольшее число перекрывающихся треков вихрей, отразившихся между 3 вихрями. Вихревой трек № 217008 является наиболее мощным антициклоническим. Вот дерьмовый сюжет:
Как вы можете видеть, 217008 находится прямо в центре окна ошибки рыбы, тогда как более крупные, более рассеянные, более слабые вихри 217343 и 39625 смещены к краю. Однако их больший размер часто приводит к тому, что они сталкиваются с вершиной, потому что у них больше места, и близость центроида не учитывается.
Итак: если fishbox перекрывается с eddy в тот же день, то включите eddy в короткий список
(thisFishRow & thisEddyRow остаются прежними). Затем: выберите вихрь из списка, основываясь на минимальном расстоянии от центроидов (что-то вместо st_join). Продолжение следует!
Edit5:
fishNearEddies <- NULL
for (thisDate in sf_df_nona$date) {
thisFishRow <- sf_df_nona[sf_df_nona$date == thisDate, ]
thisEddyRow <- sebDateSub[sebDateSub$date == thisDate, ]
overlapToday <- st_join(thisFishRow, thisEddyRow, left = F)
if (nrow(overlapToday) > 0) {
overlapEddies <- thisEddyRow[thisEddyRow$track %in% overlapToday$track,]
st_distance(x = thisFishRow, y = overlapEddies, by_element = TRUE) #0 0 due to overlap min dist = 0
} #close if
} #close for
st_distance не будет работать, поскольку пространственные объекты (буферная ячейка и буферная окружность) перекрываются, поэтому минимальное расстояние = 0. Мне также нужны центроиды для проверки, я полагаю?
Edit6: Final Edit с рабочим кодом и ответ наградить Стюарта за всю его помощь. Еще раз спасибо, сэр.
# overlap join loop####
fishNearEddies <- NULL
EdCentroidDist <- NULL
counter <- 1
ofhowmany <- length(sf_df_nona$date)
for (thisDate in sf_df_nona$date) {
thisFishRow <- sf_df_nona[sf_df_nona$date == thisDate, ] #will be 1 row per day
thisEddyRow <- sebDateSub[sebDateSub$date == thisDate, ] #all eddies that day, could be multi rows
overlapToday <- st_join(thisFishRow, thisEddyRow, left = F) #will join only if they spatial intersect and (already) time match
if (nrow(overlapToday) > 0) {
# now need to join based on closest centroid and NOT on highest overlap
overlapEddies <- thisEddyRow[thisEddyRow$track %in% overlapToday$track,] #subset TER by overlap tracks
fishcentroid <- st_centroid(thisFishRow)
eddycentroid <- st_centroid(overlapEddies)
fishEdDists <- st_distance(x = fishcentroid, y = eddycentroid, by_element = TRUE) #result vector corresponding to overlapEddies rownumbers
fishNearEddies <- rbind(fishNearEddies, overlapEddies[which.min(fishEdDists),]) #row index for overlapEddies, no index but has date
EdCentroidDist <- rbind(EdCentroidDist, as.numeric(min(fishEdDists)))
} #close if
print(paste0(counter, " of ", ofhowmany, " fish days"))
counter <- counter + 1
} #close for
if (!is.null(fishNearEddies)) { #if there are overlaps, do processing. If not it'll fail
fishNearEddies %<>% as.data.frame() %>% # convert to nonspatial so I can remove buffer column
select(-geometry) %>% # remove geometry column. Attributes still remain. Whatever?
cbind(EdCentroidDist) #add to FNE as column