Snow indicators (NDSI & SWE) Vs ASI

Afghanistan

WIP

Intro

In code below we try to assess if monthly remotely sensed snow indicators have predictive power on Agriculture Stress Index (ASI) over the March-April-May (MAM) planting season in Afghanistan.

Data Sources

• Normalized Difference Snow Index (NDSI): MODIS Terra Satellite
• Snow Water Equivalent (SWE): Famine Early Warning Systems Network (FEWS NET) Land Data Assimilation System (FLDAS)
• Agricultural Stress Index (ASI): FAO

We run the analysis for each province in Afghanistan.

Results

• These results concur with the previous analysis of Faryab province that indicated a low correlation/predictive between driest SWE years and least healthy ASI years within
• Feb-Mar snow indicators (SWE & NDSI) show a more intuitive relationship with ASI than Dec-Jan snow indicators. This could indicate that later season hydro-dynamics play a more important role in agricultural MAM drought than the actual the total quantity of snow detected in the heart of the winter months (D-J).
• SWE & NDSI are moderately correlated with each other as expected

Discussion & Next Steps

• Simple N-D-J aggregated snow indicators (NDSI & SWE) does not appear to be useful in trigger design. However,
• To confirm or support the importance of late-season hydro-dynamics we may consider further integrating snow-melt dynamics into the analysis as started output here

strip_text_size <- 8
pt_size <- 0.5
box::use(
  ../R/blob_connect
)

box::use(
  dplyr[...],
  forcats[...],
  ggplot2[...],
  gghdx[...],
  janitor[...],
  lubridate[...],
  purrr[...],
  readr[...],
  scales[...],
  stringr[...],
  glue[...],
  tidyr[...],
  scales
)

gghdx()

df_ndsi <- blob_connect$read_blob_file("DF_ADM1_MODIS_SNOW") |> 
  filter(
    # ADM1_NAME == "Faryab",
    parameter == "NDSI_Snow_Cover_mean"
  ) |> 
  clean_names()

df_swe <- blob_connect$read_blob_file("DF_ADM1_FLDAS_SWE") 



df_ndsi_proc <- df_ndsi |> 
  group_by(
    adm1_name, adm1_code,
    mo = month(date)
  ) |> 
  mutate(
    no_snow_frac = 100 - value,
    ndsi_anom = (no_snow_frac - mean(no_snow_frac)) ,
    ndsi_z = ndsi_anom /sd(no_snow_frac)
  ) |> 
  group_by(adm1_name,adm1_code) |> 
  arrange(adm1_name,adm1_code,date) |> 
  mutate(
    snow_frac_change= value -lag(value)
  ) |> 
  ungroup()

df_swe_proc <- df_swe |> 
  group_by(
    adm1_name, adm1_code,
    mo = month(date)
  ) |> 
  mutate(
    swe_anom = value - mean(value) ,
    swe_z = swe_anom /sd(value)
  ) |> 
  group_by(adm1_name,adm1_code) |> 
  arrange(adm1_name,adm1_code,date) |> 
  mutate(
    swe_change= value -lag(value)
  ) |> 
  ungroup()


df_asi <- read_csv("https://www.fao.org/giews/earthobservation/asis/data/country/AFG/MAP_ASI/DATA/ASI_Dekad_Season1_data.csv") |> 
  clean_names()

# Filter out year 2000 because it doesn't include winter season from 1999
df_ndsi_long <- df_ndsi_proc |> 
  select(
    adm1_name,
    adm1_code,
    date, 
    value, 
    no_snow_frac,
    ndsi_anom,
    ndsi_z,
    snow_frac_change
  ) |> 
  mutate(
    yr = year(date),
    mo = month(date),
    yr = ifelse(mo> 5, yr + 1, yr),
    mo = month(mo, abbr = T,label =T)
  ) #|> 
    # filter(yr>=2001)


df_asi_yr <- df_asi |> 
  mutate(
    mo = month(date)
  ) |> 
  filter(
    mo %in% c(3,4,5)
  ) |>
  group_by(
    adm1_code, mo
  ) |> 
  mutate(
    asi_anom = data-mean(data),
    asi_z = asi_anom/sd(data)
  ) |> 
  group_by(
    adm1_code,
    yr_date = floor_date(date,"year")
  ) |> 
  arrange(date) |> 
  summarise(
    asi_end = last(data), # this is 3rd dekad in MAy
    asi_max = max(data),
    asi_mean = mean(data),
    asi_max_anom = max(asi_anom),
    asi_max_z = max(asi_z),.groups = "drop"
  ) |> 
  mutate(
    yr = year(yr_date)
  )





df_swe_long <- df_swe_proc |> 
  select(
    adm1_name,
    adm1_code,
    date, 
    value,
    swe_z,
    swe_anom,
    swe_change
  ) |> 
  mutate(
    yr = year(date),
    mo = month(date), 
    yr = ifelse(mo> 5, yr + 1, yr),
    mo = month(mo, abbr = T,label =T)
  ) 
  # filter to same time range as asi full
    # filter(yr>=1984)
  
# filter as to same time range as modis
# df_asi$date |> range()
# df_asi_yr_filt_modis <- df_asi_yr |> 
#   filter(yr >=2001)
#   

df_ndsi_long_simp <- df_ndsi_long |> 
  
  select(
    adm1_name,
    adm1_code, 
    date , 
    yr,
    mo,
    raw_value = value,
    anom = ndsi_anom,
    z= ndsi_z 
  ) |> 
  mutate(
    parameter = "ndsi"
  )

df_swe_long_simp <-  df_swe_long |> 
  select(
    adm1_name,
    adm1_code, 
    date , 
    yr,
    mo,
    raw_value = value,
    anom = swe_anom,
    z= swe_z 
  ) |> 
  mutate(
    parameter = "swe"
  )

df_snow_long <- bind_rows(
  df_ndsi_long_simp,
  df_swe_long_simp
) 

df_snow_long |> 
    pivot_wider(
    id_cols = c("date","adm1_name","adm1_code"),
    values_from = raw_value,
    names_from = parameter
  ) |> 
  ggplot(
    aes(
      x = ndsi, 
      y= swe,
    )
  )+
  geom_point(
    size = pt_size
  )+
  facet_wrap(~adm1_name, scales = "free")+
  labs(
    title = "NDSI vs SWE"
  )+
  theme(
    strip.text = element_text(size = 10)
  )

df_ndsi_long |> 
  filter(yr > 2000) |> 
  left_join(
    df_asi_yr
  ) |> 
  filter(
    mo %in% c("Nov","Dec","Jan", "Feb","Mar")
    # adm1_name == "Kabul"
    ) |> 
  ggplot(
    aes(x= no_snow_frac/100, y= asi_end/100)
  )+
  geom_point(size = pt_size)+
  scale_x_continuous(
    labels = scales$label_percent()
  )+
  scale_y_continuous(
    labels = scales$label_percent()
  )+
  facet_grid(
    rows = vars(adm1_name),
    cols = vars(mo),scale = "free"
  )+
  theme(
    panel.border = element_rect(colour = "black", fill=NA),
    strip.text = element_text(size= strip_text_size),
    axis.text.x = element_text(angle = 90),
    axis.text.y = element_text(size =7)
  )+
  labs( 
    x=  "Bare - no snow (%)",
    y=  "Cumulative ASI (% stressed)",
    title = "Monthly Snow Fraction (NDSI) vs Cumulative MAM ASI"
    )

df_swe_long |> 
  left_join(
    df_asi_yr
  ) |> 
  filter(
    mo %in% c("Nov","Dec","Jan", "Feb","Mar")
    # adm1_name == "Kabul"
    ) |> 
  ggplot(
    aes(x= value, y= asi_end/100)
  )+
  geom_point(size = pt_size)+
  scale_y_continuous(
    labels = scales$label_percent()
  )+
  facet_grid(
    rows = vars(adm1_name),
    cols = vars(mo),scale = "free"
  )+
  theme(
    panel.border = element_rect(colour = "black", fill=NA),
    strip.text = element_text(size= strip_text_size),
    axis.text.x = element_text(angle = 90),
    axis.text.y = element_text(size =7)
  )+
  labs( 
    x=  "SWE Value - FLDAS",
    y=  "Cumulative ASI (% stressed)",
    title = "Monthly Snow Water Equivalent (SWE) vs Cumulative MAM ASI"
    
    )

# df_swe_monthly_wide <- df_swe_long |>
#   pivot_wider(id_cols = c("adm1_name","adm1_code","yr"),
#     names_from = mo, 
#     values_from = value # same here ... using raw value or anomaly 
#   )
# 
# df_swe_asi_wide <- left_join(df_asi_yr,df_swe_monthly_wide, by = c("adm1_code","yr")  )

ndsi_params <- 
  list(
    value = "value",
    no_snow_frac ="no_snow_frac",
    ndsi_anom = "ndsi_anom" , 
    ndsi_z="ndsi_z", 
    snow_frac_change=  "snow_frac_change"
  )

ldf_ndsi_asi_wide <- map(
  ndsi_params,\(tp){
    df_ndsi_wide <- df_ndsi_long |>
      pivot_wider(id_cols = c("adm1_name","adm1_code","yr"),
                  names_from = mo, 
                  values_from = all_of(tp) # this is the critical value -- could use z-score, anomaly, raw -value etc or bareground
      )
    
    left_join(df_asi_yr,df_ndsi_wide, by = c("adm1_code","yr")  ) |> 
    filter(
      if_all(.cols = any_of(month(c(1:12), label = T, abbr=T)),~!is.na(.x))
    )
  }
    )

swe_params <- 
  list(
    value = "value",
    anom = "swe_anom" , 
    z="swe_z", 
    change = "swe_change"
  )

ldf_swe_asi_wide <- map(
  swe_params,\(tp){
    df_swe_wide <- df_swe_long |>
      pivot_wider(id_cols = c("adm1_name","adm1_code","yr"),
                  names_from = mo, 
                  values_from = all_of(tp) # this is the critical value -- could use z-score, anomaly, raw -value etc or bareground
      )
    
    left_join(df_asi_yr,df_swe_wide, by = c("adm1_code","yr")  ) |> 
    filter(
      if_all(.cols = any_of(month(c(1:12), label = T, abbr=T)),~!is.na(.x))
    )
  }
    )

Snow vs ASI

asi_stats = list(
  `EOS` = "asi_end",
  max = "asi_max",
  `anomaly` = "asi_max_anom",
  `z-score`= "asi_max_z",
  mean = "asi_mean"
  )


lp_ndsi_heat <-  
  imap(ldf_ndsi_asi_wide,\(dft,ndsi_stat){
        map(
          asi_stats,
          \(stat_temp){
            dft_corrs_wide <- dft |> 
              group_by(
                adm1_code, adm1_name
              ) |> 
              summarise(
                across(
                  .cols = c("Jun","Jul","Aug","Sep","Oct","Nov","Dec","Jan","Feb","Mar","Apr","May"),
                  .fns = ~cor(.,!!sym(stat_temp))
                ),.groups="drop"
              )
          
            dft_long <- dft_corrs_wide |> 
              pivot_longer(cols = Jun: May) |> 
              mutate(
                name = fct_relevel(name, "Oct","Nov","Dec","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep")
              )
  
            p_title <- glue("Correlation between monthly NDSI ({ndsi_stat}) and ASI ({stat_temp}) over M-A-M")
            
            dft_long |> 
              filter(!name %in% month(6:9, abbr= T, label=T)) |> 
              ggplot(
                aes(x= name,
                    y= adm1_name,
                    fill = value)
              )+
              geom_tile()+
              scale_fill_gradient2(
                low = hdx_hex("tomato-hdx"),      # Color for negative values
                mid = "white",    # Color for zero
                high = hdx_hex("mint-hdx"),   # Color for positive values
                midpoint = 0      # Set midpoint at zero
                
              )+
              geom_tile(
                data= dft_long |> 
                  filter(adm1_name == "Faryab", 
                         name %in% c("Oct","Nov","Dec","Jan","Feb","Mar","Apr","May")),
                fill = NA, color ="black", lwd = 1.5
              )+
              
              geom_text(
                aes(label = round(value,2))
              )+
              labs(
                title = p_title,
                subtitle = "Afghanistan by Province",
                y= "Province"
              )+
              theme(
                axis.title.x = element_blank(),
                legend.title= element_blank(),
                plot.title = element_text(size = 12),
                plot.subtitle = element_text(size = 12),
                legend.text = element_text(angle=90)
              )
          }
        )
    }
  )
      
      

# lp_ndsi_heat$ndsi_anom$EOS
# lp_ndsi_heat$ndsi_z$EOS
# lp_ndsi_heat$ndsi_z$max

Below we have correlated monthly NDSI values against End of Season (EOS). As NDSI is a measure of % snow we’d expect an negatibve relationship in key winter months where high % snow fraction should be associated with low EOS ASI and visa versa.

However the heat map below shows generally a positive relationship between % snow in Dec & Jan & EOS ASI. This is the opposite of what we’d intuitively expect, but is also confirmed by SW.

lp_ndsi_heat$value$EOS

lp_swe_heat <-  
  imap(ldf_swe_asi_wide,\(dft,swe_stat){
        map(
          asi_stats,
          \(stat_temp){
            dft_corrs_wide <- dft |> 
              group_by(
                adm1_code, adm1_name
              ) |> 
              summarise(
                across(
                  .cols = c("Nov","Dec","Jan","Feb","Mar"),
                  .fns = ~cor(.,!!sym(stat_temp))
                ),.groups="drop"
              )
            
            dft_long <- dft_corrs_wide |> 
              pivot_longer(cols = Nov:Mar) |> 
              mutate(
                name = fct_relevel(name, "Nov","Dec","Jan","Feb","Mar")
              )
  
            p_title <- glue("Correlation between monthly SWE ({swe_stat}) and ASI ({stat_temp}) over M-A-M")
            
            dft_long |> 
              ggplot(
                aes(x= name,
                    y= adm1_name,
                    fill = value)
              )+
              geom_tile()+
              scale_fill_gradient2(
                low = hdx_hex("tomato-hdx"),      # Color for negative values
                mid = "white",    # Color for zero
                high = hdx_hex("mint-hdx"),   # Color for positive values
                midpoint = 0      # Set midpoint at zero
                
              )+
              geom_tile(
                data= dft_long |> 
                  filter(adm1_name == "Faryab", 
                         name %in% c("Nov","Dec","Jan","Feb","Mar")),
                fill = NA, color ="black", lwd = 1.5
              )+
              
              geom_text(
                aes(label = round(value,2))
              )+
              labs(
                title = p_title,
                subtitle = "Afghanistan by Province",
                y= "Province"
              )+
              theme(
                axis.title.x = element_blank(),
                legend.title= element_blank(),
                plot.title = element_text(size = 12),
                plot.subtitle = element_text(size = 12),
                legend.text = element_text(angle=90)
              )
          }
        )
    }
  )

where we see the same relationship with high water/snow values in Dec & Jan correlated with higher Ag stress (EOS).

Potential explanation

Of note is the fact that by Feb-Mar the relationship shifts back to the direction we would expect: More snow-less stress. This could indicate that it’s not simply the total snowfall impacting water availability and crop health, but perhaps there is a later season dynamic (snow melt?) that impacts water availability. Perhaps if there is some (x?) amount of snow available, it keeps a consistent supply of water to crops, but if it disappears too early, than the agriculture suffers

lp_swe_heat$value$EOS      

ASI explore

df_asi |> 
  filter(
    !year(date) %in% 2024
  ) |> 
  # filter(province == "Faryab") |> 
  group_by(
    adm1_code, province, mo = month(date, abbr = T, label = T)
  ) |> 
  # summarise(
  #   asi = median(data)
  # ) |> 
  ggplot(
    aes(x= mo, y= data)
  )+
  geom_boxplot()+ 
  facet_wrap(~province)

df_snow |> 
  filter(
    !year(date) %in% 2024
  ) |> 
  group_by(
    adm1_name, adm1_code, mo = month(date, abbr = T, label = T)
  ) |> 

  ggplot(
    aes(x= mo, y= value/100)
  )+
  geom_boxplot()+ 
  labs(
    title = "SWE by Province & Month"
  )+
  facet_wrap(~adm1_name,scale = "free")+
  scale_y_continuous(label = scales::label_percent())+
  theme(
    axis.text.x=element_text(angle=90)
  )

df_snow |> 
  filter(
    !year(date) %in% 2024
  ) |> 
  mutate(
    mo = month(date, abbr = T, label = T)
  ) |> 

  ggplot(
    aes(x= mo, y= value/100)
  )+
  geom_boxplot()+ 
  labs(
    title = "SWE by Province & Month"
  )+
  facet_wrap(~adm1_name,scale = "free")+
  scale_y_continuous(label = scales::label_percent())+
  theme(
    axis.text.x=element_text(angle=90)
  )

df_snow |> 
  filter(
    !year(date) %in% 2024
  ) |> 
  mutate(
    mo = month(date, abbr = T, label = T)
  ) |> 

  ggplot(
    aes(x= mo, y= no_snow_frac/100)
  )+
  geom_boxplot()+ 
  labs(
    title = "SWE by Province & Month"
  )+
  facet_wrap(~adm1_name,scale = "free")+
  scale_y_continuous(label = scales::label_percent())+
  theme(
    axis.text.x=element_text(angle=90)
  )

box::use(terra[...])
box::use(sf[...])

bbox <- st_bbox(gdf)

r <- rast(
  x = "/vsicurl/https://agricultural-production-hotspots.ec.europa.eu/data/indicators_sm/sm_combined/sm_combined_20231221.tif",
  win = bbox
  )

plot(r)
plot(gdf,add= T, col = NA)