Submission to Applied Sciences

R/fig5.R

@@ -8,7 +8,6 @@
 
 rm(list = ls())
 
-
 simulo <- read_feather("data/clean/simulo/compute-canada/simulo_4_lambertian_sources.feather") %>% 
   filter(between(pixel_distance_to_center, 0, 50))
 
@@ -21,55 +20,81 @@ simulo <- simulo %>%
   summarise(value = mean(value)) %>%
   ungroup()
 
+# Fit a Guassian curve on raidance data to remove noise -------------------
+
+fit_gaussian <- function(df, depth) {
+
+  print(depth)
+
+  mod <- minpack.lm::nlsLM(
+    radiance ~ a1 * exp(-((mid_distance - b1) ^ 2 / (2 * c1 ^ 2))) + k,
+    data = df,
+    start = list(
+      a1 = mean(df$radiance),
+      b1 = 0,
+      c1 = 5,
+      k = mean(df$radiance)
+    ),
+    lower = c(
+      a1 = 0,
+      b1 = 0,
+      c1 = 0,
+      k = 0
+    ),
+    upper = c(
+      a1 = max(df$radiance),
+      b1 = 0,
+      c1 = 100,
+      k = max(df$radiance)
+    )
+  )
+
+  return(mod)
+}
+
+simulo <- simulo %>% 
+  bind_rows(mutate(simulo, mid_distance = -mid_distance)) %>% 
+  distinct()
+
 simulo <- simulo %>% 
   spread(source, value)
 
-res <- simulo %>%
+simulo <- simulo %>%
   group_by(depth) %>%
   nest() %>%
-  mutate(mod_gam = map(data,  ~gam(.$radiance ~ s(.$mid_distance, fx = FALSE, k=10, bs = "cr")) , data = .)) %>% 
-  mutate(pred_sf = map(data, function(x) {
-
-    sf <- smooth.spline(x$mid_distance, x$radiance, spar = 0.75)
-    tibble(pred_sf = sf$y)
-
-  })) %>% 
-  mutate(pred_gam = map(mod_gam, predict)) %>% 
-  unnest(data, pred_gam, pred_sf) 
-
-p1 <- res %>%
-  filter(depth == 5) %>% 
+  mutate(mod = map2(data, depth, fit_gaussian)) %>% 
+  mutate(pred = map2(data, mod, add_predictions)) %>% 
+  unnest(pred)
+
+p1 <- simulo %>% 
+  filter(depth %in% c(0.5, 5, 10, 15, 20, 25)) %>% 
   ggplot(aes(x = mid_distance, y = intensity)) +
-  geom_point() +
-  ylab("Number of irradiance photon") +
+  geom_line() +
+  facet_wrap(~depth, scales = "free") +
   xlab("Distance from the center of the melt pond (m)") +
-  theme(text = element_text(size = 10))
+  ylab("Number of irradiance photons") +
+  theme(legend.position = "none")
 
-p2 <- res %>%
-  filter(depth == 5) %>% 
+p2 <- simulo %>% 
+  filter(depth %in% c(0.5, 5, 10, 15, 20, 25)) %>% 
   ggplot(aes(x = mid_distance, y = radiance)) +
-  geom_point() +
-  geom_line(aes(y = pred_gam), color = "red") +
-  ylab("Number of radiance photon") +
+  geom_line(size = 0.1) +
+  geom_line(aes(y = pred, color = "red")) +
+  facet_wrap(~depth, scales = "free") +
   xlab("Distance from the center of the melt pond (m)") +
-  theme(text = element_text(size = 10))
+  ylab("Number of radiance photons") +
+  theme(legend.position = "none")
 
-p <- cowplot::plot_grid(p1, p2, labels = "AUTO", align = "hv")  
-ggsave("graphs/supp_fig_5.pdf", width = 7, height = 3, device = cairo_pdf)
+p <- cowplot::plot_grid(p1, p2, labels = "AUTO", align = "hv", ncol = 1)  
+ggsave("graphs/supp_fig_5.pdf", width = 7, height = 6, device = cairo_pdf)
 
-res <- res %>% 
-  select(-pred_sf, -radiance) %>% 
-  rename(radiance = pred_gam) %>% 
-  gather(source, value, intensity, radiance)
 
-# res %>%
-#   ggplot(aes(x = mid_distance, y = value)) +
-#   geom_line() +
-#   # geom_line(data = res, aes(x = mid_distance, y = pred_gam, color = "pred_gam")) +
-#   # geom_line(data = res, aes(x = mid_distance, y = pred_sf, color = "pred_sf")) +
-#   facet_wrap(~depth + source, scales = "free_y")
-
-simulo <- res %>% 
+## Replace radiance data with smoothed value (Gaussian fits)
+simulo <- simulo %>% 
+  select(-radiance) %>% 
+  rename(radiance = pred) %>% 
+  gather(source, value, intensity, radiance) %>% 
+  filter(mid_distance >= 0) %>% 
   mutate(source = ifelse(source == "radiance", "Radiance (Lu)", "Irradiance (Ed)"))
 
 ## So we can reuse it later
@@ -84,9 +109,18 @@ df <- simulo %>%
   mutate(interpolated = map(interpolated, ~akima::interp2xyz(., data.frame = TRUE))) %>%
   unnest(interpolated)
 
+# df <- df %>%
+#   group_by(source) %>%
+#   mutate(z = z / max(z))
+
+## Idea of Marcel
+df <- df %>%
+  group_by(source) %>%
+  mutate(z = z / z[x == 50 & y == 0.5])
+
 df <- df %>%
   group_by(source) %>%
-  mutate(z = z / max(z))
+  mutate(z = log(z))
 
 ## Just "mirror" the x-ditances, this makes a better looking plot
 df <- df %>%

R/fig6-9.R

@@ -2,29 +2,7 @@ rm(list = ls())
 
 source("R/simulo/simulo_utils.R")
 
-# simulo <-
-#   read_feather("data/clean/simulo/compute-canada/simulo_45degrees.feather") %>%
-#   mutate(source = ifelse(source == "radiance", "Radiance (Lu)", "Irradiance (Ed)")) %>%
-#   filter(pixel_distance_to_center <= 50)
-
-# simulo <- read_feather("data/clean/simulo/compute-canada/simulations-with-15m-melt-pond-4-lambertian-sources.feather")
 simulo <- read_feather("data/clean/simulo/compute-canada/simulo_4_lambertian_sources_smoothed_radiance.feather")
-# simulo <- simulo %>%
-#   simulo_xy_to_polar(min_distance = 10, max_distance = 60) %>% 
-#   smooth_pixels()
-
-# simulo <- simulo %>%
-#   mutate(iso_distance = cut_width(pixel_distance_to_center, width = 0.1)) %>%
-#   separate(iso_distance, into = c("start_distance", "end_distance"), sep = ",") %>%
-#   mutate_at(vars(start_distance, end_distance), parse_number) %>%
-#   mutate(mid_distance = start_distance + (end_distance - start_distance) / 2) %>%
-#   group_by(depth, source, mid_distance) %>%
-#   summarise(value = mean(value)) %>%
-#   ungroup()
-
-# simulo <- simulo %>% 
-#   select(-value) %>% 
-#   rename(value = smoothed_value)
 
 reference_profile <- simulo %>%
   mutate(class_25_percent = as.character(cut(mid_distance, breaks = c(0, sqrt(25 / 0.25)), include.lowest = TRUE, right = TRUE, dig.lab = 5))) %>%
@@ -34,8 +12,6 @@ reference_profile <- simulo %>%
   mutate(class_05_percent = as.character(cut(mid_distance, breaks = c(0, sqrt(25 / 0.05)), include.lowest = TRUE, right = TRUE, dig.lab = 5))) %>%
   mutate(class_01_percent = as.character(cut(mid_distance, breaks = c(0, sqrt(25 / 0.01)), include.lowest = TRUE, right = TRUE, dig.lab = 5)))
 
-# unique(reference_profile$class_0_40)
-
 reference_profile <- reference_profile %>%
   gather(class_distance, range, starts_with("class")) %>%
   drop_na(range)
@@ -109,20 +85,70 @@ averaged_simulo <- simulo %>%
   summarise(value = mean(value)) %>%
   ungroup()
 
+# Fit a Guassian curve on raidance data to remove noise -------------------
+
+fit_gaussian <- function(df, depth) {
+
+  print(depth)
+
+  mod <- minpack.lm::nlsLM(
+    radiance ~ a1 * exp(-((mid_distance - b1) ^ 2 / (2 * c1 ^ 2))) + k,
+    data = df,
+    start = list(
+      a1 = mean(df$radiance),
+      b1 = 0,
+      c1 = 5,
+      k = mean(df$radiance)
+    ),
+    lower = c(
+      a1 = 0,
+      b1 = 0,
+      c1 = 5, # Limit to 5 because some have very high value outliers at 0 m (center of the melt pond)
+      k = 0
+    ),
+    upper = c(
+      a1 = max(df$radiance),
+      b1 = 0,
+      c1 = Inf,
+      k = max(df$radiance)
+    ),
+    control = nls.lm.control(maxiter = 1024)
+  )
+
+  return(mod)
+}
+
+
+averaged_simulo <- averaged_simulo %>% 
+  bind_rows(mutate(averaged_simulo, mid_distance = -mid_distance)) %>% 
+  distinct()
+
 averaged_simulo <- averaged_simulo %>% 
   spread(source, value)
 
 averaged_simulo <- averaged_simulo %>%
   group_by(depth) %>%
   nest() %>%
-  mutate(mod_gam = map(data,  ~gam(.$radiance ~ s(.$mid_distance, fx = FALSE, k=10, bs = "cr")) , data = .)) %>% 
-  mutate(pred_gam = map(mod_gam, predict)) %>% 
-  unnest(data, pred_gam) 
-
+  mutate(mod = map2(data, depth, fit_gaussian)) %>% 
+  mutate(pred = map2(data, mod, add_predictions)) %>% 
+  unnest(pred)
+
+averaged_simulo %>% 
+  filter(depth %in% c(0.5, 5, 10, 15, seq(20, 25, by = 0.5))) %>% 
+  ggplot(aes(x = mid_distance, y = radiance)) +
+  geom_line(size = 0.1) +
+  geom_line(aes(y = pred, color = "red")) +
+  facet_wrap(~depth, scales = "free") +
+  xlab("Distance from the center of the melt pond (m)") +
+  ylab("Number of radiance photons") +
+  theme(legend.position = "none")
+
+## Replace radiance data with smoothed value (Gaussian fits)
 averaged_simulo <- averaged_simulo %>% 
   select(-radiance) %>% 
-  rename(radiance = pred_gam) %>% 
+  rename(radiance = pred) %>% 
   gather(source, value, intensity, radiance) %>% 
+  filter(mid_distance >= 0) %>% 
   mutate(source = ifelse(source == "radiance", "Radiance (Lu)", "Irradiance (Ed)"))
 
 ## Calculate K only starting at 5 meters (ice ridge)
@@ -258,14 +284,12 @@ ggsave("graphs/fig9.pdf", plot = p, device = cairo_pdf, height = 3, width = 7)
 
 ## Stats for the paper
 
-# res %>%
-#   spread(type, integral) %>% 
-#   mutate(relative_error = (reference - predicted) / reference) %>% 
-#   select(-predicted, -reference) %>% 
-#   spread(source, relative_error) %>% 
-#   janitor::clean_names() %>% 
-  # filter(str_detect(range, "25%|1%")) %>%
-#   mutate()
+res %>%
+  spread(type, integral) %>% 
+  mutate(relative_error = (reference - predicted) / reference) %>% 
+  group_by(source) %>% 
+  summarise(mean(relative_error) * 100)
+
 
 res %>% 
   spread(type, integral) %>%

R/main.R

@@ -21,6 +21,7 @@ library(sf)
 library(ggsn)
 library(modelr)
 library(glue)
+library(scales)
 
 rm(list = ls())
 graphics.off()

R/supp_fig_1.R

@@ -32,13 +32,13 @@ p <- ggplot(baffin) +
   ) +
   scale_x_continuous(breaks = seq(-90, -30, by = 10)) +
   scale_y_continuous(breaks = seq(50, 80, by = 3)) +
-  annotate("text", x = -795630.7, y = -2338551, label = "Qikiqtarjuaq ice camp", size = 2.5, vjust = -1, hjust = -0.1, family = "IBM Plex Sans") +
+  annotate("text", x = -795630.7, y = -2438551, label = "Qikiqtarjuaq\n ice camp", size = 3, vjust = -1, hjust = -0.1, family = "IBM Plex Sans") +
   annotate("text", x = -2174308, y = -2533284, label = "Hudson Bay", family = "IBM Plex Sans") +
   annotate("text", x = -76357.57, y = -2184929, label = "Greenland", angle = 90, family = "IBM Plex Sans") +
   annotate("text", x = -1720983, y = -3435623, label = "Québec", family = "IBM Plex Sans") +
   annotate("text", x = -607354.6, y = -1915067, label = "Baffin Bay", family = "IBM Plex Sans") +
-  theme(axis.title = element_blank()) +
-  scalebar(x.min = -1720983, x.max = -70357.57, y.min = -3535623, y.max = -1584929, dd2km = FALSE, dist = 100, st.size = 2)
+  theme(axis.title = element_blank())
+# +
+#   scalebar(x.min = -1720983, x.max = -70357.57, y.min = -3535623, y.max = -1584929, dd2km = FALSE, dist = 100, st.size = 2)
 
-ggsave("graphs/supp_fig_1.pdf")
-embed_fonts("graphs/supp_fig_1.pdf")
+ggsave("graphs/supp_fig_1.pdf", device = cairo_pdf)

R/supp_fig_4.R

@@ -8,4 +8,31 @@
 
 rm(list = ls())
 
+df1 <- readxl::read_excel("data/raw/angular_distribution.xlsx", n_max = 2, na = "NaN") %>% 
+  janitor::remove_empty() %>% 
+  gather(angle, value) %>% 
+  set_names(c("angle", "Measured under ice\ndownwelling radiance distribution")) %>% 
+  mutate_all(parse_number)
 
+df2 <- readxl::read_excel("data/raw/angular_distribution.xlsx", skip = 4, n_max = 2) %>% 
+  janitor::remove_empty() %>% 
+  gather(angle, value) %>% 
+  set_names(c("angle", "Emission source chosen\nfor the simulation")) %>% 
+  mutate_all(parse_number)
+
+df <- bind_rows(df1, df2) %>% 
+  gather(source, value, -angle) %>% 
+  drop_na()
+
+df %>% 
+  ggplot(aes(x = angle, y = value, color = source)) +
+  geom_path() +
+  xlab(expression("Azimuthal angle"~(degree))) +
+  ylab("Normalized raidance") +
+  theme(legend.title = element_blank()) +
+  theme(legend.justification = c(0, 0)) +
+  theme(legend.position = c(0.01, 0.01)) +
+  scale_x_continuous(limits = c(0, 90), breaks = seq(0, 90, by = 10)) +
+  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.1))
+
+ggsave("graphs/supp_fig_4.pdf", device = cairo_pdf, width = 3 * 1.61803398875, height = 3)