Regression on when there will be no sea ice

SeaIceRegression.r

#get the the most recent day
#1. Get Data: I downloaded the data from these ftps
#ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/north/daily/data/NH_seaice_extent_final_v2.csv
#ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/north/daily/data/NH_seaice_extent_nrt_v2.csv
nisdc1 = read.csv("NH_seaice_extent_final_v2.csv",skip=2,col.names=col.names)
nisdc2 = read.csv("NH_seaice_extent_nrt_v2.csv",skip=2,col.names=col.names)

#2. load libraries
library(dplyr)
library(tidyr)
library(ggplot2)
library(animation)
library(ggthemes)
#3. get data in right format
col.names <- c("Year", "Month", "Day", "Extent", "Missing", "Source_Data")

nisdc <- data.frame(Year=c(nisdc1$Year,nisdc2$Year),
                    Month=c(nisdc1$Month,nisdc2$Month),
                    Day=c(nisdc1$Day,nisdc2$Day),
                    Extent=c(nisdc1$Extent,nisdc2$Extent),
                    Missing=c(nisdc1$Missing,nisdc2$Missing))

#get today in different years                    
today <- filter(nisdc, Day == 17 & Month == 9)
head(today)
#fill in the blanks from early on when they only got every second day
today2 <- filter(nisdc, Day == 18 & Month == 9 & Year==1981)
today3 <- rbind(today, today2)
today2 <- filter(nisdc, Day == 18 & Month == 9 & Year==1983)
today3 <- rbind(today3, today2)
today2 <- filter(nisdc, Day == 18 & Month == 9 & Year==1984)
today3 <- rbind(today3, today2)
today2 <- filter(nisdc, Day == 18 & Month == 9 & Year==1986)
today3 <- rbind(today3, today2)

# Basic scatterplot
p1 <- ggplot(today3, aes(x = Year, y = Extent))
p1 + geom_point()
# Linear regression to predict when there will be 0 sea ice on this date (in 2068). Actual predictions are much
# more clever then this. So don't believe this date. This is just showing what will happen if things keep going the same way
# but all sorts of non linear feedback systems come into play.
# this article https://judithcurry.com/2016/09/18/is-the-arctic-sea-ice-spiral-of-death-dead/ claims there is no positive feedback loop in melting which wold be good news
p2 <- p1 + geom_point(color="black") + xlim(1979,2068) +stat_smooth(method="lm",fullrange=TRUE)

p3 <- p2 + labs(x="Year", 
                y = "Ice Extent",parse = TRUE) 
#p4 <- p3 + theme_bw() 

 p5 <- p3 + annotate(x=1985, y=8.2, geom="text", label = "10^6*km^2",parse = TRUE, size = 4) +
                ggtitle(label = "Arctic Sea Ice Size\n on Sept 17th each year") 
p6 <- p5 + annotate(x=2048, y=-.5, geom="text", label = "Chart by @iamreddave | Data source: nsidc.org", size = 3)
p7<-p6 + theme_economist() + scale_colour_economist() 
ggsave("Arctic.png", width=6, height=6, dpi=100)