Update with GeoPressure/GeoPressureR#52

Rafnuss · Rafnuss · commit 4d95950d5f3b · 2022-10-02T16:36:19.000-04:00
diff --git a/analysis/1-pressure.R b/analysis/1-pressure.R
@@ -134,13 +134,17 @@ if (debug) {
 
   # Query timeserie of pressure based on these path
   pressure_timeserie <- geopressure_ts_path(path, pam_short$pressure, include_flight = c(0, 1))
+  pressure_ts_bind <- do.call("rbind", pressure_timeserie) %>%
+    filter(!is.na(sta_id))
 
   # Test 3 ----
   p <- ggplot() +
     geom_line(data = pam$pressure, aes(x = date, y = obs), colour = "grey") +
     geom_point(data = subset(pam$pressure, isoutlier), aes(x = date, y = obs), colour = "black") +
-    geom_line(data = pressure_na, aes(x = date, y = obs, color = factor(sta_id)), size = 0.5) +
-    geom_line(data = do.call("rbind", pressure_timeserie) %>% filter(!is.na(sta_id)), aes(x = date, y = pressure0, col = factor(sta_id)), linetype = 2) +
+    geom_line(data = pam$pressure %>%
+                mutate(obs = ifelse(isoutlier | sta_id == 0, NA, obs)),
+              aes(x = date, y = obs, color = factor(sta_id)), size = 0.5) +
+    geom_line(data = pressure_ts_bind %>% filter(sta_id > 0), aes(x = date, y = pressure0, col = factor(sta_id)), linetype = 2) +
     theme_bw() +
     scale_colour_manual(values = col) +
     scale_y_continuous(name = "Pressure(hPa)")
@@ -149,16 +153,19 @@ if (debug) {
 
 
   # Test 4 ----
-  par(mfrow = c(5, 6), mar = c(1, 1, 3, 1))
-  for (i_r in seq_len(length(pressure_timeserie))) {
-    if (!is.null(pressure_timeserie[[i_r]])) {
-      i_s <- unique(pressure_timeserie[[i_r]]$sta_id)
-      df3 <- merge(pressure_timeserie[[i_r]], subset(pam$pressure, !isoutlier & sta_id == i_s), by = "date")
-      df3$error <- df3$pressure0 - df3$obs.x
-      hist(df3$error, main = i_s, xlab = "", ylab = "")
-      abline(v = 0, col = "red")
-    }
-  }
+  # You might also want to adjust the value of s in `geopressure_prob_map()` based on the SD value
+  # of these histogram gpr$prob_map_s
+  pam$pressure %>%
+    left_join(pressure_ts_bind %>% dplyr::select(c("date","pressure0")), by="date") %>%
+    mutate(diff=ifelse(is.na(pressure0), 0, obs-pressure0)) %>%
+    filter(sta_id > 0 & !isoutlier) %>%
+    group_by(sta_id) %>%
+    mutate(sta_id = paste0(sta_id, " (SD=",round(sd(diff),2)," ; N=",n(),")")) %>%
+    ggplot( aes(x=diff)) +
+    geom_histogram(aes(y=(..count..)/tapply(..count..,..PANEL..,sum)[..PANEL..]), binwidth=.4) +
+    facet_wrap(~sta_id) +
+    scale_x_continuous(name = "Pressure Geolocator - best match ERA5 (hPa)") +
+    scale_y_continuous(name = "Normalized histogram")
 
   # Map the most likely position
   sta_duration <- unlist(lapply(pressure_prob, function(x) {
diff --git a/analysis/2-light.R b/analysis/2-light.R
@@ -130,38 +130,43 @@ if (debug) {
   hist(z, freq = F)
   lines(fit_z, col = "red")
 
-  # Compare with pressure
-  load(paste0("data/1_pressure/", gdl, "_pressure_prob.Rdata"))
-  dur <- unlist(lapply(pressure_prob, function(x) difftime(metadata(x)$temporal_extent[2], metadata(x)$temporal_extent[1], units = "days")))
-  long_id <- which(dur > 5)
-
-  par(mfrow = c(2, 3))
-  for (i_s in long_id) {
-    twl_fl <- twl %>%
-      filter(!deleted) %>%
-      filter(twilight > pressure_timeserie[[i_s]]$date[1] & twilight < tail(pressure_timeserie[[i_s]]$date, 1))
-    sun <- solar(twl_fl$twilight)
-    z_i <- refracted(zenith(sun, pressure_timeserie[[i_s]]$lon[1], pressure_timeserie[[i_s]]$lat[1]))
-    hist(z_i, freq = F, main = paste0("sta_id=", i_s, " n=", nrow(twl_fl)))
-    lines(fit_z, col = "red")
-  }
+  # Compute twilight at the best location of pressure
+  path <- geopressure_map2path(pressure_prob)
+  twl_path <- left_join(twl, path) %>%
+    mutate(
+      twilight = twilight(twilight,
+                          lon = lon, lat = lat, rise = rise, zenith = 96)
+    ) %>%
+    filter(!is.na(twilight))
 
-  # Light comparison
+  # Twilight comparison with pressure based position (probably not correct for short stationay period)
   lightImage(
     tagdata = raw_geolight,
     offset = gpr$shift_k / 60 / 60
   )
   tsimagePoints(twl$twilight,
-    offset = gpr$shift_k / 60 / 60, pch = 16, cex = 1.2,
-    col = ifelse(twl$deleted, "grey20", ifelse(twl$rise, "firebrick", "cornflowerblue"))
+                offset = -gpr$shift_k / 60 / 60, pch = 16, cex = 1.2,
+                col = ifelse(twl$deleted, "grey20", ifelse(twl$rise, "firebrick", "cornflowerblue"))
+  )
+  tsimageDeploymentLines(twl_path$twilight,
+                         lon = twl_path$lon, twl_path$lat,
+                         offset = -gpr$shift_k / 60 / 60, lwd = 3, col = adjustcolor("orange", alpha.f = 0.5)
   )
-  for (i_s in long_id) {
+
+  # Check histogram of zenith angle for each stationary period based on pressure.
+  # Filter long stationary period
+  sta_id_long <- pam$sta %>%
+    filter(difftime(end,start,units = "days")>5) %>%
+    .$sta_id
+
+  par(mfrow = c(2, 3))
+  for (i_s in sta_id_long) {
     twl_fl <- twl %>%
-      filter(twilight > pressure_timeserie[[i_s]]$date[1] & twilight < tail(pressure_timeserie[[i_s]]$date, 1))
-    tsimageDeploymentLines(twl_fl$twilight,
-      lon = pressure_timeserie[[i_s]]$lon[1], pressure_timeserie[[i_s]]$lat[1],
-      offset = gpr$shift_k / 60 / 60, lwd = 3, col = adjustcolor("orange", alpha.f = 0.5)
-    )
+      filter(!deleted) %>%
+      filter(sta_id==i_s)
+    z_i <- refracted(zenith(solar(twl_fl$twilight), path$lon[path$sta_id==i_s], path$lat[path$sta_id==i_s]))
+    hist(z_i, freq = F, main = paste0("sta_id=", i_s, " n=", nrow(twl_fl)))
+    lines(fit_z, col = "red")
   }
 }
 
@@ -236,7 +241,8 @@ if (debug) {
 }
 
 # Save ----
-save(twl,
+save(
+  twl,
   light_prob,
   z,
   fit_z,
diff --git a/analysis/3-static.R b/analysis/3-static.R
@@ -84,7 +84,7 @@ tmp[!is.na(tmp)] <- 0
 tmp[lat_calib_id, lon_calib_id] <- 1
 values(static_prob[[1]]) <- tmp
 
-if (!is.na(gpr$calib_2_start)) {
+if (!is.na(gpr$calib_2_start) & abs(difftime(gpr$calib_2_end, gpr$crop_end, units = "days"))<3) {
   if (!is.na(gpr$calib_2_lat)) {
     lon_calib_id <- which.min(abs(gpr$calib_2_lon - lon))
     lat_calib_id <- which.min(abs(gpr$calib_2_lat - lat))
@@ -101,22 +101,67 @@ if (!is.na(gpr$calib_2_start)) {
 path <- geopressure_map2path(static_prob)
 static_timeserie <- geopressure_ts_path(path, pam$pressure)
 
-
 if (debug) {
   # GeopressureViz
-  geopressureviz <- list(
+  path_modified <- geopressureviz(
     pam = pam,
     static_prob = static_prob,
     pressure_prob = pressure_prob,
     light_prob = light_prob,
     pressure_timeserie = static_timeserie
   )
-  save(geopressureviz, file = "~/geopressureviz.RData")
-
-  appDir <- system.file("geopressureviz", package = "GeoPressureR")
-  shiny::runApp(appDir, launch.browser = getOption("browser"))
 
+  # Pressure
+  path_modified_timeserie <- geopressure_ts_path(path_modified, pam$pressure)
+
+  path_modified_ts_bind <- do.call("rbind", path_modified_timeserie) %>%
+    filter(!is.na(sta_id))
+
+  # To make the labeling easier, you can replace pam$pressure$obs by the difference between
+  # pam$pressure$obs and pressure_timeserie$pressure0 to see the anomalies.
+  # Because trainset_read does not read the actual value of obs, but simply the label, it won't
+  # impact your code
+  pam_diff <- pam
+  pam_diff$pressure <- pam_diff$pressure %>%
+    left_join(path_modified_ts_bind %>% dplyr::select(c("date","pressure0")), by="date") %>%
+    rename(obs_ref = pressure0)
+
+  trainset_write(pam_diff, "data/1_pressure/labels/", filename = paste0(pam$id, "_act_pres"))
+
+  # We can automatically extract some outlier based on s value
+  # pam_diff_pressure <- pam_diff_pressure %>%
+  #  mutate( isoutlier = ifelse(!isoutlier&abs(diff)>4*gpr$prob_map_s&sta_id>0, TRUE, isoutlier))
+  # message(sum(pam_diff_pressure$isoutlier)-sum(pam$pressure$isoutlier), " new outlier automatically added")
+  twl_path <- left_join(twl, path_modified) %>%
+    mutate(
+      twilight = twilight(twilight,
+                          lon = lon, lat = lat, rise = rise, zenith = 96)
+    ) %>%
+    filter(!is.na(twilight))
+
+  write.csv(
+    df <- rbind(
+      data.frame(
+        series = ifelse(twl$rise, "Rise", "Set"),
+        timestamp = strftime(twl$twilight, "%Y-%m-%dT00:00:00Z", tz = "UTC"),
+        value = (as.numeric(format(twl$twilight, "%H")) * 60 + as.numeric(format(twl$twilight, "%M"))
+                 + gpr$shift_k / 60 + 60 * 12) %% (60 * 24),
+        label = ifelse(is.na(twl$delete), "", ifelse(twl$delete, "Delete", ""))
+      ),
+      data.frame(
+        series = ifelse(twl_path$rise, "Set_ref", "Rise_ref"),
+        timestamp = strftime(twl_path$twilight, "%Y-%m-%dT00:00:00Z", tz = "UTC"),
+        value = (as.numeric(format(twl_path$twilight, "%H")) * 60 + as.numeric(format(twl_path$twilight, "%M"))
+                 + gpr$shift_k / 60 + 60 * 12) %% (60 * 24),
+        label = ""
+      )
+    ),
+    paste0("data/2_light/labels/", gpr$gdl_id, "_light.csv"),
+    row.names = FALSE
+  )
+}
 
+if (debug) {
   # Check 1
   static_prob_n <- lapply(static_prob, function(x) {
     probt <- raster::as.matrix(x)
diff --git a/analysis/4-basic-graph.R b/analysis/4-basic-graph.R
@@ -25,12 +25,6 @@ grl$p <- grl$ps * flight_prob(grl$gs, method = "gamma", shape = 7, scale = 7, lo
 
 
 
-
-
-
-
-
-
 # Marginal map ----
 static_prob_marginal <- graph_marginal(grl)
 
@@ -58,72 +52,15 @@ path_sim <- graph_simulation(grl, nj = nj)
 
 
 if (debug) {
-
-  # Rapid visual check
-  sta_duration <- unlist(lapply(static_prob_marginal, function(x) {
-    as.numeric(difftime(metadata(x)$temporal_extent[2], metadata(x)$temporal_extent[1], units = "days"))
-  }))
-
-  m <- leaflet(width = "100%") %>%
-    addProviderTiles(providers$Stamen.TerrainBackground) %>%
-    addFullscreenControl()
-  for (i in seq_len(nj)) {
-    m <- m %>%
-      addPolylines(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = 0.7, weight = 1, color = "#808080") %>%
-      addCircles(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = 1, weight = 1, color = "#000")
-  }
-
-  m %>%
-    addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
-    addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#000", weight = sta_duration^(0.3) * 10)
-
-
-  # Light comparison
   load(paste0("data/2_light/", gdl, "_light_prob.Rdata"))
-  raw_geolight <- pam$light %>%
-    transmute(
-      Date = date,
-      Light = obs
-    )
-  lightImage(
-    tagdata = raw_geolight,
-    offset = gpr$shift_k / 60 / 60
-  )
-  tsimagePoints(twl$twilight,
-    offset = -gpr$shift_k / 60 / 60, pch = 16, cex = 1.2,
-    col = ifelse(twl$deleted, "grey20", ifelse(twl$rise, "firebrick", "cornflowerblue"))
-  )
-  for (ts in shortest_path_timeserie) {
-    twl_fl <- twl %>%
-      filter(twilight > ts$date[1] & twilight < tail(ts$date, 1))
-    tsimageDeploymentLines(twl_fl$twilight,
-      lon = ts$lon[1], ts$lat[1],
-      offset = gpr$shift_k / 60 / 60, lwd = 3, col = adjustcolor("orange", alpha.f = 0.5)
-    )
-  }
-
-  # In depth analysis with GeoPressureViz
-  load(paste0("data/1_pressure/", gdl, "_pressure_prob.Rdata"))
-  sta_marginal <- unlist(lapply(static_prob_marginal, function(x) raster::metadata(x)$sta_id))
-  sta_pres <- unlist(lapply(pressure_prob, function(x) raster::metadata(x)$sta_id))
-  sta_light <- unlist(lapply(light_prob, function(x) raster::metadata(x)$sta_id))
-  pressure_prob <- pressure_prob[sta_pres %in% sta_marginal]
-  light_prob <- light_prob[sta_light %in% sta_marginal]
-
-
-  geopressureviz <- list(
+  geopressureviz(
     pam = pam,
     static_prob = static_prob,
     static_prob_marginal = static_prob_marginal,
     pressure_prob = pressure_prob,
     light_prob = light_prob,
     pressure_timeserie = shortest_path_timeserie
   )
-  save(geopressureviz, file = "~/geopressureviz.RData")
-
-  shiny::runApp(system.file("geopressureviz", package = "GeoPressureR"),
-    launch.browser = getOption("browser")
-  )
 }
 
 
diff --git a/analysis/5-3-wind-graph_analyse.R b/analysis/5-3-wind-graph_analyse.R
@@ -56,32 +56,14 @@ path_sim <- graph_simulation(grl, nj = nj)
 
 
 if (debug) {
-
-  # In depth analysis with GeoPressureViz
-  load(paste0("data/1_pressure/", gdl, "_pressure_prob.Rdata"))
-  load(paste0("data/2_light/", gdl, "_light_prob.Rdata"))
-  load(paste0("data/3_static/", gdl, "_static_prob.Rdata"))
-
-  sta_marginal <- unlist(lapply(static_prob_marginal, function(x) raster::metadata(x)$sta_id))
-  sta_pres <- unlist(lapply(pressure_prob, function(x) raster::metadata(x)$sta_id))
-  sta_light <- unlist(lapply(light_prob, function(x) raster::metadata(x)$sta_id))
-  pressure_prob <- pressure_prob[sta_pres %in% sta_marginal]
-  light_prob <- light_prob[sta_light %in% sta_marginal]
-
-
-  geopressureviz <- list(
+  geopressureviz(
     pam = pam,
     static_prob = static_prob,
     static_prob_marginal = static_prob_marginal,
     pressure_prob = pressure_prob,
     light_prob = light_prob,
     pressure_timeserie = shortest_path_timeserie
   )
-  save(geopressureviz, file = "~/geopressureviz.RData")
-
-  shiny::runApp(system.file("geopressureviz", package = "GeoPressureR"),
-    launch.browser = getOption("browser")
-  )
 }
 
 
diff --git a/analysis/99-combined.R b/analysis/99-combined.R
