Converting Time to Date to create the Year / Month variable Each Month Corresponds to a season and it’s numeric Month Value. ex. January = 1
df$Time <- as.Date(df$Time)
df$Year <- year(df$Time)
df$Month <- month(df$Time)
df$YearMonth <- floor_date(df$Time, "month")
df$Season <- case_when(
df$Month %in% c(12) ~ "Winter1",
df$Month %in% c(1) ~ "Winter2",
df$Month %in% c(2) ~ "Winter3",
df$Month %in% c(3) ~ "Spring1",
df$Month %in% c(4) ~ "Spring2",
df$Month %in% c(5) ~ "Spring3",
df$Month %in% c(6) ~ "Summer1",
df$Month %in% c(7) ~ "Summer2",
df$Month %in% c(8) ~ "Summer3",
df$Month %in% c(9) ~ "Fall1",
df$Month %in% c(10) ~ "Fall2",
df$Month %in% c(11) ~ "Fall3",
TRUE ~ NA_character_
)Breaking down monthly averages for each station, and for numeric columns.
# Calculating the monthly average per station
df_monthly <- df %>%
group_by(IDStations, YearMonth) %>%
summarise(across(where(is.numeric), mean, na.rm = TRUE), .groups = "drop")
# The monthly average for numeric columns
df_monthly <- df %>%
group_by(YearMonth) %>%
summarise(across(where(is.numeric), mean, na.rm = TRUE), .groups = "drop")WE_temp_2m and WE_precipitation_t to time series and seasonal adjustments frequency = 12 ; means monthly
temp_ts <- ts(df_monthly$WE_temp_2m, start = c(min(df_monthly$Year), min(df_monthly$Month)), frequency = 12)
precip_ts <- ts(df_monthly$WE_precipitation_t, start = c(min(df_monthly$Year), min(df_monthly$Month)), frequency = 12)
#Seasonally Adjust
adj_temp <- seas(temp_ts)
adj_precip <- seas(precip_ts)
#Extracting adjusted data
adj_temp_data <- final(adj_temp)
adj_precip_data <- final(adj_precip)
# moving the adjusted data back to monthly df
df_monthly$WE_temp_2m_adj <- adj_temp_data
df_monthly$WE_precipitation_t_adj <- adj_precip_data~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~ Importing Agricultural Output Data From Lombardia Regions. Yearly_Production_Quintals is the agricultural output value for crops collected by ISTAT.it in Italy for Lombardia. Crops for this analysis were carefully chosen to match specific seasons. Namely, Barley, Common Wheat, Common Ryegrass, Grain Maize and Rice.
Yearly_Production_Quintals <- read_excel("/Users/redaabouzaid/Desktop/All Files/Projects/Agriculture_Polution/Yearly_Production_Quintals.xlsx")
Yearly_Forage_Quintals <- read_excel("/Users/redaabouzaid/Desktop/All Files/Projects/Agriculture_Polution/Seasonality/Forage_Quintals.xlsx", sheet = 2)## New names:
## • `forage barley` -> `forage barley...4`
## • `forage barley` -> `forage barley...6`Yearly_Quintals <- merge(Yearly_Forage_Quintals, Yearly_Production_Quintals)
# Convert Barley into time series, yearly frequency
low_freq_Barley <- ts(Yearly_Quintals$`Barley`, start = min(df$Year), frequency = 1)
low_freq_Wheat <- ts(Yearly_Quintals$`Common wheat`, start = min(df$Year), frequency = 1)
low_freq_Forage <- ts(Yearly_Quintals$`common ryegrass`, start = min(df$Year), frequency = 1)
low_freq_Maize <- ts(Yearly_Quintals$`Grain maize`, start = min(df$Year), frequency = 1)
low_freq_Rice <- ts(Yearly_Quintals$`Rice`, start = min(df$Year), frequency = 1) ~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~ Apply Chow Lin Method to disaggregate time series (using td() temporal disaggregation function) The Chow-lin or Chow Lin Maxlog method is applied to the Yearly outputs for each crop in order to disaggregate the Crop Data. The Predictor used in order to disaggregate is proxy for crop growth and is correlated with crop growth. In this case, we utilize temperature.
# Use Precipitation and Temperature as predictors
high_freq_series_combined <- cbind(df_monthly$WE_temp_2m_adj, df_monthly$WE_precipitation_t_adj)
result_Barley <- td(low_freq_Barley ~ high_freq_series_combined, to = "month", method = "chow-lin-maxlog")
result_Wheat <- td(low_freq_Wheat ~ high_freq_series_combined, to = "month", method = "chow-lin-maxlog")
result_Forage <- td(low_freq_Forage ~ high_freq_series_combined, to = "month", method = "chow-lin-maxlog")
result_Maize <- td(low_freq_Maize ~ high_freq_series_combined, to = "month", method = "chow-lin-maxlog")
result_Rice <- td(low_freq_Rice ~ high_freq_series_combined, to = "month", method = "chow-lin-maxlog")
# Extract from time series then convert to data frame + filter
disaggregated_Barley <- result_Barley$values
disaggregated_Wheat <- result_Wheat$values
disaggregated_Forage <- result_Forage$values
disaggregated_Maize <- result_Maize$values
disaggregated_Rice <- result_Rice$values
disaggregated_agrioutputdf <- data.frame(Date = as.Date(as.yearmon(time(disaggregated_Barley), "%Y-%m")), Series = coredata(disaggregated_Barley))
disaggregated_agrioutputdf <- cbind(disaggregated_agrioutputdf, disaggregated_Wheat, disaggregated_Forage,
disaggregated_Maize, disaggregated_Rice)
disaggregated_agrioutputdf <- rename(disaggregated_agrioutputdf, Time = Date, Barley = Series, Wheat = disaggregated_Wheat,
Forage = disaggregated_Forage, Maize = disaggregated_Maize, Rice = disaggregated_Rice)
# Adding YearMonth variable to agricultural output
disaggregated_agrioutputdf$YearMonth <- floor_date(disaggregated_agrioutputdf$Time, "month")
# Adding a Season column to the df_monthly then joining
df_monthly$Season <- df$Season[match(df_monthly$YearMonth, df$YearMonth)]
disaggregated_agrioutputdf <- left_join(disaggregated_agrioutputdf, df_monthly %>% select(YearMonth, Season), by = "YearMonth")Weighting By Season ~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~ Finally, weights for each season are applied to each crop. Not all crops grow or are harvested all throughout the year. As such, for each crop, there must be weights applied per month in order to achieve realistic monthly outputs. The weights are selected from
Assessing in-season crop classification performance using satellite data: A test case in Northern Italy by: Azar, Villa Stroppiana and Crema. Linked Below: https://www.researchgate.net/figure/Typical-crop-calendar-for-major-cultivation-in-Lombardy-Northern-Italy_fig2_305784624
This paper provides a crop calendar for the Lombardy - Northern Italy region. Weights are assigned only for harvesting seasons for each crop. All months have weekly blocks. The total number of harvesting weeks for each crops becomes it’s weight. Every month without a weekly block has a value of 0.
Further analysis can be done by assigning wants to developing weeks and finding relationships between relationship and development weeks. Or # of development weeks. Given the Crop calendar is updated yearly.
WintCrop_adj_factor <- c("Winter1" = 0, "Winter2" = 0, "Winter3" = 0, "Spring1" = 0, "Spring2" = 0, "Spring3" = 4/6, "Summer1" = 2/6, "Summer2" = 0, "Summer3" = 0, "Fall1" = 0, "Fall2" = 0, "Fall3" = 0)
ForCrop_adj_factor <- c("Winter1"=0,"Winter2"=0,"Winter3"=0,"Spring1"=0,"Spring2"=0,"Spring3"=1/4,"Summer1"=1/4,"Summer2"=1/4,"Summer3"=1/4,"Fall1"=0,"Fall2"=0,"Fall3"=0)
Maize_adj_factor <- c("Winter1"=0,"Winter2"=0,"Winter3"=0,"Spring1"=0,"Spring2"=0,"Spring3"=0,"Summer1"=0,"Summer2"=0,"Summer3"=1/5,"Fall1"=4/5,"Fall2"=0,"Fall3"=0)
Rice_adj_factor <- c("Winter1"=0,"Winter2"=0,"Winter3"=0,"Spring1"=0,"Spring2"=0,"Spring3"=0,"Summer1"=0,"Summer2"=0,"Summer3"=0,"Fall1"=1/5,"Fall2"=4/5,"Fall3"=0)
# Adjust the seasonal_adjustment_factor to sum = 1
WintCrop_adj_factor <- WintCrop_adj_factor / sum(WintCrop_adj_factor)
ForCrop_adj_factor <- ForCrop_adj_factor / sum(ForCrop_adj_factor)
Maize_adj_factor <- Maize_adj_factor / sum(Maize_adj_factor)
Rice_adj_factor <- Rice_adj_factor / sum(Rice_adj_factor)
# Create a Year column in agrioutput after disaggreagation
disaggregated_agrioutputdf$Year <- year(disaggregated_agrioutputdf$Time)
disaggregated_agrioutputdf$Month <- month(disaggregated_agrioutputdf$Time)
# Calculate total annual output for each year
annual_output <- disaggregated_agrioutputdf %>%
group_by(Year) %>%
summarise(Total_Annual_Output_Barley = sum(Barley), Total_Annual_Output_Wheat = sum(Wheat),
Total_Annual_Output_Rice = sum(Rice), Total_Annual_Output_Maize = sum(Maize),
Total_Annual_Output_Forage = sum(Forage), .groups = "drop")
# Merge the annual output to the original dataframe
disaggregated_agrioutputdf <- left_join(disaggregated_agrioutputdf, annual_output, by = "Year")
# Calculate the seasonally adjusted SumAll
disaggregated_agrioutputdf <- disaggregated_agrioutputdf %>%
mutate(Barley_adj_season = Total_Annual_Output_Barley * WintCrop_adj_factor[Season],
Rice_adj_season = Total_Annual_Output_Rice * Rice_adj_factor[Season],
Maize_adj_season = Total_Annual_Output_Maize * Maize_adj_factor[Season],
Forage_adj_season = Total_Annual_Output_Forage * ForCrop_adj_factor[Season],
Wheat_adj_season = Total_Annual_Output_Wheat * WintCrop_adj_factor[Season])
# Putting Seasonally Adjusted Data into one dataset
dis_seas_data <- data.frame(disaggregated_agrioutputdf$Time, disaggregated_agrioutputdf$Month, disaggregated_agrioutputdf$Season, disaggregated_agrioutputdf$Barley_adj_season,
disaggregated_agrioutputdf$Rice_adj_season, disaggregated_agrioutputdf$Maize_adj_season,
disaggregated_agrioutputdf$Forage_adj_season, disaggregated_agrioutputdf$Wheat_adj_season)
dis_seas_data <- rename(dis_seas_data, Time = disaggregated_agrioutputdf.Time, Season = disaggregated_agrioutputdf.Season,
Barley = disaggregated_agrioutputdf.Barley_adj_season, Rice = disaggregated_agrioutputdf.Rice_adj_season,
Maize = disaggregated_agrioutputdf.Maize_adj_season, Forage = disaggregated_agrioutputdf.Forage_adj_season,
Wheat = disaggregated_agrioutputdf.Wheat_adj_season, Month = disaggregated_agrioutputdf.Month)
dis_seas_data <- dis_seas_data %>%
mutate(sumAll = Maize + Wheat + Forage + Rice + Barley)
dis_seas_data <- dis_seas_data %>%
filter(Month >= 5 & Month <= 10)
# Add a YearMonth column to dis_seas_data for merging
dis_seas_data$YearMonth <- floor_date(dis_seas_data$Time, "month")~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~ Merging All Data Frames to create data file necessary for analysis.
df2 <- subset(df_monthly, subset = Month %in% c(1,2,3,4,5,6,7,8,9,10,11,12))
# Merging dis_seas_data and df2 based on YearMonth
Amerged_df <- merge(df2, dis_seas_data, by = "YearMonth", all.x = TRUE)
write.csv(df2, file = "/Users/redaabouzaid/Desktop/All Files/1.csv")
write.csv(dis_seas_data, file = "/Users/redaabouzaid/Desktop/All Files/2.csv")
write.csv(Amerged_df, file = "/Users/redaabouzaid/Desktop/All Files/Merged_Agrimonia.csv")
Amerged_df2 <- Amerged_df %>%
arrange(YearMonth) %>%
mutate(
AQ_pm10 = lag(AQ_pm10, 5),
AQ_pm25 = lag(AQ_pm25, 5),
AQ_co = lag(AQ_co, 5),
AQ_nh3 = lag(AQ_nh3, 5),
)
view(Amerged_df2)