计算R中融雪的开始日期[关闭]答案

【问题标题】：computing onset date of snowmelt in R [closed]计算R中融雪的开始日期[关闭]
【发布时间】：2018-06-14 19:27:37
【问题描述】：

从 1950 年到 2017 年，我有这种格式的每日温度 Data

我需要计算融雪开始日期定义为当日气温低于 0 摄氏度时，在 3 月至 5 月的最后 5 天之后，日气温高于 0 摄氏度的第一天. 到目前为止我的代码：

  df1<-read.csv("temp.csv")
  require(dplyr)
  # applying the condition to check each temperature value
  df1$boolean<- ifelse(df1$temp<0.0 , 1, 0)

  #computing the total sum < 0 and the start and end date
  snow<-df1 %>%
  mutate(boolean = ifelse(is.na(boolean), 0, boolean)) %>%
  group_by(group = cumsum(c(0, diff(boolean) != 0))) %>%
   filter(boolean == 1 & n() > 1) %>%
   summarize("Start Date"=min(as.character(date)),
        "End Date"=max(as.character(date)),
        "Length of Run"=n()) %>%
   ungroup() %>%
  select(-matches("group"))
colnames(snow)[3] <- 'length'

# subset length that greater >5
obs<-subset(snow,length >=5)

上面的代码为我提供了部分解决方案（如果进一步手动编辑，我将获得符合我定义的理想解决方案）我只对每年的一个发病日期感兴趣。我需要一些关于如何编辑此代码以根据上述定义计算开始日期的进一步指导。

我有多个位置，因此手动编辑这不是理想的解决方案。您的帮助将不胜感激

【问题讨论】：

更多实际包含您正在寻找的案例的数据会有所帮助
请阅读How to Create a Complete, Minimal, and Verifiable Example 并更新您的帖子。我们特别需要您发布足够的数据，以便我们弄清楚如何计算融雪开始日期。
我添加了一部分数据
您能否更具体地谈谈极端情况？比如三五月份没有低于 0 的五天？
由于位置的原因，我预计我的最后五天将在三月或四月或五月发生。因此不必太担心

标签： r dplyr threshold

【解决方案1】：

我们在 (1) 中假设融化日必须发生在 3 月、4 月或 5 月，并且在 (2) 中只有 5 个零下天发生在 3 月、4 月、5 月，但融化日可能发生在 6 月，例如.

1) 定义 df2，即 df1 加上附加列：月份、年份和代码，其中代码为 0，如果日期不在三月、四月、五月，否则为 1，如果 temp = 0。

如果最近 6 个日期的代码为 1、1、1、1、1、2，则现在使用 df2 对返回 TRUE 的代码运行 rollapplyr，否则为 FALSE。取 TRUE 行，每年只保留最后一行。将其加入所有年份的数据框，以便在输出中生成任何缺失年份的 NA。

library(zoo)

df2 <- df1 %>%
    mutate(Date = as.Date(Date), month = as.numeric(format(Date, "%m")), 
           year = as.numeric(format(Date, "%Y")),
           code = (month %in% 3:5) * ((temp < 0) + 2 * (temp >= 0)),
           OK = rollapplyr(code, 6, identical, c(1, 1, 1, 1, 1, 2), fill = FALSE))

df2 %>%
       filter(OK) %>%
       filter(!duplicated(year, fromLast = TRUE)) %>%
       right_join(unique(df2["year"]), by = "year") %>%
       select(year, Date)

给予：

   year       Date
1  1950 1950-05-24
2  1951 1951-05-21
3  1952 1952-05-28
4  1953 1953-05-15
5  1954 1954-05-28
6  1955 1955-05-14
7  1956 1956-05-27
8  1957 1957-05-17
9  1958 1958-05-21
10 1959       <NA>
11 1960 1960-05-26
12 1961 1961-05-16
13 1962 1962-05-19
14 1963 1963-05-13
15 1964 1964-05-27
16 1965 1965-05-20
17 1966 1966-05-26
18 1967 1967-05-26
19 1968 1968-05-27
20 1969 1969-05-30
21 1970 1970-05-21

2) 在 (1) 中，我们假设融化开始日必须在 3 月、4 月或 5 月，但这里我们假设只有零下以下的日子在该范围内，并且融化开始日可能进一步向外延伸。

计算与 (1) 中的相同，只是代码现在使得 1 表示 3 月、4 月或 5 月的零下温度，2 表示任何时间高于零的任何温度（不仅仅是在 3 月、4 月和 5 月） 0 是其他任何东西。我们将代码折叠成一个字符串（每个日期一个字符）并在其上使用正则表达式来查找由 5 个 1 组成的子字符串，后跟任何内容，直到我们到达下一个 2。我们按照 (1) 处理其余部分，除了现在我们不需要加入，因为总会有一个融化开始的日子。如果没有连接，我们现在可以将其表示为单个管道。

df1 %>%
    mutate(Date = as.Date(Date), month = as.numeric(format(Date, "%m")), 
           year = as.numeric(format(Date, "%Y")),
           code = (month %in% 3:5) * (temp < 0) + 2 * (temp >= 0),
           OK = { g <- gregexpr("1{5}.*?2", paste(code, collapse = ""))[[1]]
                  seq_along(code) %in% (g + attr(g, "match.length") - 1) }) %>%
    filter(OK) %>%
    filter(!duplicated(year, fromLast = TRUE)) %>%
    select(year, Date)

给予：

   year       Date
1  1950 1950-05-24
2  1951 1951-06-01
3  1952 1952-05-28
4  1953 1953-05-15
5  1954 1954-05-28
6  1955 1955-05-14
7  1956 1956-05-27
8  1957 1957-05-17
9  1958 1958-05-21
10 1959 1959-06-02
11 1960 1960-05-26
12 1961 1961-05-16
13 1962 1962-05-19
14 1963 1963-06-01
15 1964 1964-05-27
16 1965 1965-05-20
17 1966 1966-05-26
18 1967 1967-05-26
19 1968 1968-05-27
20 1969 1969-05-30
21 1970 1970-05-21

【讨论】：

我今天从你的回答中学到了很多。我只是想确认 c(1, 1, 1, 1, 1, 2) 是相同的可选参数，并且代码中的每 6 个元素都使用相同的可选参数进行检查。对吗？
同rollapplyr(codes, 6, function(x) identical(x, c(1, 1, 1, 1, 1, 2), fill = FALSE)。这是identical 的第一个参数是codes 上的窗口，第二个是 c(1, 1, 1, 1, 1, 2)`。
非常感谢您的解释。代码让我确认了我的疑惑。
添加了 (2)，它改变了 (1) 中的假设，因此融化开始日可能会发生得更远。
谢谢！ 2) 可以很好地满足我的需求！

【解决方案2】：

tidyverse 中的简单解决方案。

library(tidyverse)
library(lubridate)


readxl::read_excel("temp.xlsx") -> df1

df1 %>%
  mutate(year = year(Date),
         month = month(Date)) %>%
  group_by(year) %>%
  mutate(
    below_0 = as.numeric(temp < 0),
    streak5 = cumsum(below_0) - cumsum(lag(below_0, 5, 0)),
    onset = month %in% c(3, 4, 5) & lag(streak5) == 5 & below_0 == 0) %>% 
  filter(onset) %>%
  summarise(Date = last(Date))

给予

# A tibble: 20 x 2
    year       Date
   <dbl>     <dttm>
 1  1950 1950-05-24
 2  1951 1951-05-21
 3  1952 1952-05-28
 4  1953 1953-05-15
 5  1954 1954-05-28
 6  1955 1955-05-14
 7  1956 1956-05-27
 8  1957 1957-05-17
 9  1958 1958-05-21
10  1960 1960-05-26
11  1961 1961-05-16
12  1962 1962-05-19
13  1963 1963-05-13
14  1964 1964-05-27
15  1965 1965-05-20
16  1966 1966-05-26
17  1967 1967-05-26
18  1968 1968-05-27
19  1969 1969-05-30
20  1970 1970-05-21

我希望代码或多或少能解释自己，streak5 是之前温度低于 0 的天数，onset 实现了问题中给出的标准，summarise 选择了给定年份的最后一个日期。

【讨论】：

【解决方案3】：

rle() 来救援！

library(broom)
library(tidyverse)

temp <- read_csv("temp.csv")

在阅读这个辅助函数之前最好先阅读下面的管道。

我们每年：

采用高于/低于 0 的行程编码
第一个为 TRUE (
获取下一个索引
如果过多（没有符合条件的日期）返回 NA
否则返回那个日期

因此：

mk_runs <- function(xdf) {

  r <- rle(xdf$below_0) take the T/F RLE
  pos <- which(r$values & r$length>=5)[1] # find the first one meeting criteria
  idx <- (sum(r$lengths[1:pos]))+1 # sum the lengths up until this point and add 1 to get to the first > 0 day

  if (idx > nrow(xdf)) { # if past our date range return NA
    data_frame(year=xdf$year[1], date=NA)
  } else {
    xdf[idx, c("year", "date")]
  }

}

我们需要整理数据：

separate(temp, Date, c("month", "day", "year")) %>%
  mutate_all(as.numeric) %>% 
  mutate(year = ifelse(year >=50, 1900+year, 2000+year)) %>% 
  mutate(date = as.Date(sprintf("%04d-%02d-%02d", year, month, day))) %>% 
  mutate(month = lubridate::month(date)) %>% 
  mutate(below_0 = temp < 0) %>% 
  filter(month >= 3 & month <=5) %>% 
  group_by(year) %>% # year groups
  arrange(date) %>%  # in order
  do(mk_runs(.)) %>% # see above function
  print(n=21)
## # A tibble: 21 x 2
## # Groups:   year [21]
##     year       date
##    <dbl>     <date>
##  1  1950 1950-04-30
##  2  1951 1951-05-21
##  3  1952 1952-05-28
##  4  1953 1953-05-15
##  5  1954 1954-05-28
##  6  1955 1955-05-14
##  7  1956 1956-05-02
##  8  1957 1957-05-07
##  9  1958 1958-04-27
## 10  1959         NA
## 11  1960 1960-04-24
## 12  1961 1961-05-16
## 13  1962 1962-05-19
## 14  1963 1963-05-13
## 15  1964 1964-05-20
## 16  1965 1965-05-20
## 17  1966 1966-05-07
## 18  1967 1967-04-27
## 19  1968 1968-05-10
## 20  1969 1969-05-22
## 21  1970 1970-05-21

【讨论】：

我无法使用 tidyverse 方式处理所有问题。为您提供整洁的解决方案。
恐怕结果不正确（例如1950）。

【解决方案4】：

这是另一个尝试。在我的第一步中，我首先创建了两个新列（即年和月）。然后，我过滤了 3 月到 5 月之间的数据。然后，我为温度高于 0 摄氏度的行创建了索引号。这个过程每年进行一次。由于您需要在温度高于零的那些日子之前连续五天，因此需要忽略等于/小于 5 的索引号。这是在外部if_else() 的真实条件下完成的if_else()。

在我的第二步中，我选择使用SOfun 的包，它是由splitstackshape 的作者开发的。你可以从 github 下载这个包。 getMyRows() 正在做什么； 1) 它通过指定pattern 来确定应考虑哪些行，2) 从 1) 中标记的行中获取一定范围的行，以及 3) 创建一个列表。这里range = -5:0 表示我选择了目标行的前五行，以及目标行本身。

在第三步中，我使用两个逻辑条件对mylist 进行了子集化。 !is.na(x$ind[6]) 检查 ind 的第 6 个元素是否不是 NA，all(x$temp[1:5] < 0) 检查 temp（温度）的第 1-5 个元素是否都小于零。 Filter() 选择满足两个逻辑条件的列表元素。然后，我从每个数据帧中提取了第 6 行，因为那是目标行。我绑定了列表，按年份对数据进行了分组，并使用slice() 选择了每年的第一个观察值。

library(devtools)

install_github("mrdwab/overflow-mrdwab")
install_github("mrdwab/SOfun")

library(overflow)
library(SOfun)
library(readxl)
library(dplyr)

# Part 1
mydf <- read_excel("temp.xlsx") %>%
        mutate(year = as.numeric(format(Date, "%Y")),
               month = as.numeric(format(Date, "%m"))) %>%
        filter(between(month, 3, 5)) %>%
        group_by(year) %>%
        mutate(ind = if_else(temp > 0, 
                     {ind <- row_number()
                      if_else(ind <= 5, NA_integer_, ind)},
                      NA_integer_)) %>%
        ungroup

# Part 2
mylist <- getMyRows(mydf,
                    pattern = which(complete.cases(mydf$ind)),
                    range = -5:0, isNumeric = TRUE)

# Part 3
Filter(function(x) !is.na(x$ind[6]) & all(x$temp[1:5] < 0), mylist) %>%
lapply(function(x) x[6, ]) %>%
bind_rows %>%
group_by(year) %>%
slice(1) %>%
select(Date)

    year Date               
   <dbl> <dttm>             
 1  1950 1950-04-30 00:00:00
 2  1951 1951-05-21 00:00:00
 3  1952 1952-05-28 00:00:00
 4  1953 1953-05-15 00:00:00
 5  1954 1954-05-28 00:00:00
 6  1955 1955-05-14 00:00:00
 7  1956 1956-05-02 00:00:00
 8  1957 1957-05-07 00:00:00
 9  1958 1958-04-27 00:00:00
10  1960 1960-04-24 00:00:00
11  1961 1961-05-16 00:00:00
12  1962 1962-05-19 00:00:00
13  1963 1963-05-13 00:00:00
14  1964 1964-05-20 00:00:00
15  1965 1965-05-20 00:00:00
16  1966 1966-05-07 00:00:00
17  1967 1967-04-27 00:00:00
18  1968 1968-05-10 00:00:00
19  1969 1969-05-22 00:00:00
20  1970 1970-05-21 00:00:00

【讨论】：

恐怕结果不正确（例如1950）。
@liborm 感谢您的评论。我查了资料。 1950 年 4 月 30 日，气温为 1.50。前连续五天温度低于零。该日期是 1950 年 3 月至 5 月期间温度首次保持在零以上的日期。您能用您的数据确认吗？
@liborm the first day when daily temperature is above 0 C, following the last five-day period between March and May, when the daily temperature is below 0 C 似乎表示第一天温度高于零，而前连续五天的温度低于零。如果这个解释不正确，我认为 OP 需要澄清这个定义。
@liborm 我明白了。我认为我们直观的想法可能不会反映在数据集中。如果是这样，我们要质疑的是数据集的可靠性，而不是判断什么结果是正确的或不正确的。
heh .. 对不起，我反应匆忙，以为我在问题下方的线程中。您的解决方案的问题很简单，与 1959 无关 - OP 要求 following the last five-day period。您未能确定 1950 年的最后一个时期，您只找到了一个 5 天的时期。