makeDF_aroundKatsuura.R

#############################################################
#
# makeDf_araoundKatuura.R
#
# UNDER CONSTRUCTION !
#
#############################################################

pacman::p_load(
  tidyverse,
  sf,
  rmapshaper,
  reshape2,
  here,
  ncdf4)


JP_geo <- "https://github.com/igproj-fusion/R-gis/raw/main/japan.geojson"
JP <- read_sf(JP_geo) |> 
  filter(PREF != "Okinawa") |> 
  ms_dissolve()


kat_lon <- 140 + 18.7 / 60
kat_lat <- 35 + 9 / 60
theta <- seq(-pi, pi/2, length = 100)
lon <- cos(theta) * 5 + kat_lon
lat <- sin(theta) * 5 + kat_lat
polygon_matrix = cbind(
  x = c(kat_lon, lon, kat_lon),
  y = c(kat_lat, lat, kat_lat)
)
circle <- st_sfc(st_polygon(list(polygon_matrix))) |>  
  st_set_crs(st_crs(JP))

#########################################################



# Get filenames
# All the daily .nc files we downloaded:
all_fnames <- fs::dir_ls(here("raw"), 
                         recurse = TRUE, glob = "*.nc")

#
# lon = (Var2 - 1) * 0.25 + 0.125
#     = Var2 / 4 -  0.125
# Var2 = (lon - 0.125) / 0.25 + 1
#     = lon * 4 + 0.5
#
# lat = (Var1 - 1) * 0.25 + 0.125 - 90
#     = Var1 / 4 - 90.125
# Var1 = (lat + 90 - 0.125) / 0.25 + 0.125 
#     = lat * 4 + 360.5
#
#

Var2_center = kat_lon * 4 + 0.5
Var1_center = kat_lat  * 4 + 360.5




############################


nc4data <- nc_open(all_fnames[1])
sst <- as.vector(ncvar_get(nc4data, "sst")) 
DF <- melt(matrix(sst, ncol = 1440, byrow =　T)) |> 
  filter(value != is.na(value)) |> 
  filter(Var1 >= Var1_center - 20 & Var1 <= Var1_center + 20) |> 
  filter(Var2 >= Var2_center - 20 & Var2 <= Var2_center + 20) |> 
  mutate(lon = Var2 / 4 -  0.125) |> 
  mutate(lat = Var1 / 4 - 90.125) |> 
  select(lon, lat, sst = value) |> 
  rowwise() |> 
  mutate(point = list(st_point(c(lon, lat)))) |> 
  filter(st_intersects(point, circle, sparse = FALSE)[1,1])


ggplot() + 
  geom_point(data = DF, aes(x = lon, y = lat, color = sst), size = 2) + 
  geom_sf(data = JP, fill = "white") +
  geom_point(aes(kat_lon, kat_lat), color = "red", size = 1.5) +
  xlim(c(125, 150)) + ylim(c(25, 47)) +
  labs(x = "", y = "")