SWAT-plus_LOCA1-SWATplus_single_2012HUC.Rmd

```{r}


  library(package = "stringr")
  library(package = "forcats")
  library(package = "readr")
  library(package = "tidyverse")
  library(package = "lubridate")
  library(package = "RCurl")
  library(package = "elevatr")
  library(package = "tidync")



file="./LOCA2_MASKS.nc"

HUC08_MRB_Code_Frame = tidync(file) %>%
  hyper_tibble()  %>%
  filter(!is.na(LOCA2_MASK)) %>%
  rename(elevation = surface_altitude)


```

```{r}

thredds_root= "http://kyrill.ias.sdsmt.edu:8080/thredds/fileServer/LOCA2/Specific_Regional_Aggregate_Sets/USGS_HUC08_Basins/R_Daily_Files/"

                
get_metadata = "http://kyrill.ias.sdsmt.edu:8080/thredds/fileServer/LOCA2/Specific_Regional_Aggregate_Sets/USGS_HUC08_Basins/USGS_HUC08_LUT.RData"

load(url(get_metadata), verbose=TRUE)


hucs_to_process = c("10120202")

hucs_to_process = c("10120110")

HUC08_MRB_LUT = USGS_HUC08_LUT %>% filter(huc08 %in% hucs_to_process)





```

```{r}



```

```{r}
ggplot(data = HUC08_MRB_Code_Frame) +
  aes(x = lon,
      y = lat,
      color = elevation,
      fill  = elevation) + 
  geom_point()


ggplot(data = HUC08_MRB_Code_Frame) +
  aes(x = lon,
      y = lat,
      color = LOCA2_HUC08,
      fill  = LOCA2_HUC08) + 
  geom_point( show.legend = FALSE)


 for (basin in c("10120202","10120110"))
  {   
       spatial = HUC08_MRB_Code_Frame %>% 
      filter(LOCA2_HUC08 == basin) 

    plot = ggplot(data = spatial) +
      aes(x = lon,
          y = lat,
          color = elevation,
          fill  = elevation) + 
      geom_point()
    
    print(plot)
    
       
    spatial = HUC08_MRB_Code_Frame %>% 
      filter(LOCA2_HUC08 == basin) %>%
      group_by(LOCA2_HUC08) %>% 
      summarize(max_lat = max(lat),
                min_lat = min(lat),
                max_lon = max(lon),
                min_lon = min(lon),
                lat = mean(lat),
                lon = mean(lon),
                elev = mean(elevation))
    
    print(spatial)
    
    line_1 = str_c("ID",
                  "NAME",
                  "LAT",
                  "LONG",
                  "ELEVATION",
                  sep=',')
    
    line_2 = str_c(1,
                  "loca_pcp",
                  sprintf("%8.5f", spatial$lat),
                  sprintf("%9.5f", spatial$lon),
                  sprintf("%9.5f", spatial$elev),
                  sep=',')
 }


        
        line_1 = "ID,NAME,LAT,LONG,ELEVATION"
        line_2 = "1,loca_pcp,44.40000,-103.47000,1005.79999"


```



```{r}


  Output_Scenarios = c("SSP2-4.5",
                       "SSP3-7.0",
                       "SSP5-8.5")


  for (basin in hucs_to_process)
  {

    metadata = HUC08_MRB_LUT %>% filter(huc08 == basin)
    


    
    spatial = HUC08_MRB_Code_Frame %>% 
      filter(LOCA2_HUC08 == basin) %>%
      group_by(LOCA2_HUC08) %>% 
      summarize(max_lat = max(lat),
                min_lat = min(lat),
                max_lon = max(lon),
                min_lon = min(lon),
                lat = mean(lat),
                lon = mean(lon),
                elev = mean(elevation))
    
    


    File_URL = str_c(thredds_root,
                     "LOCA2_V1_HUC08_",
                     basin,
                     ".RData",
                     sep = "")

    my.connection = url(description = File_URL)
      load(file = my.connection)
      close(con = my.connection)
      remove( my.connection)
      
    loca2_daily = loca2_daily %>%
      mutate(Ensemble = str_c(Model, ".", Member, sep = ""))

    ensembles =  unique(loca2_daily$Ensemble)
    print(ensembles)
    print(length(ensembles))
    for (ensemble in unique(loca2_daily$Ensemble))
    {

      for (scenario in Output_Scenarios)
      {
        
        out_dir = str_c("~/Desktop/swat2012_prelim_meteo_HUC08_LOCA2_10120110/",
                     ensemble,
                     "__",
                     str_replace(string = scenario,
                                 pattern     = " ",
                                 replacement = ""),
                     "/",
                     sep="")
        

                     
        system(str_c("mkdir -pv ",out_dir))

        subset = loca2_daily %>%
                       filter(((Scenario  == "Historical") |
                               (Scenario  == scenario)) &
                              (Ensemble   == ensemble) &
                              (Percentile == "MEAN")) %>%
                       mutate(Year   = year(Time),
                              Julian = yday(Time))
        


        #
        # Write the Daily Max/Min Temperature File
        #
        
        filename_hedt = str_c(out_dir,
                             "tmp.txt",
                             sep = "")
        
        line_1 = "ID,NAME,LAT,LONG,ELEVATION"
        line_2 = str_c(1,
                  "loca_tmp",
                  sprintf("%8.5f", spatial$lat),
                  sprintf("%9.5f", spatial$lon),
                  sprintf("%9.5f", spatial$elev),
                  sep=',')

        fileConn<-file(filename_hedt)
          writeLines(c(line_1,line_2), fileConn)
          close(fileConn)



        filename_tmp = str_c(out_dir,
                             "loca_tmp.txt",
                             sep = "")
        


        line_1 = "19500101"
                      

        fileConn<-file(filename_tmp)
          writeLines(c(line_1), fileConn)
          close(fileConn)
          
          subperiod_print = subset %>%
            select("tasmax","tasmin")

          
          write_delim(x     = subperiod_print, 
                      file  =  filename_tmp, 
                      delim =       ",",
                      append=TRUE) 
     





        #
        # Write the Daily Rainfall File
        #        
        filename_hedp = str_c(out_dir,
                             "pcp.txt",
                             sep = "")
        
        line_1 = "ID,NAME,LAT,LONG,ELEVATION"
        line_2 = str_c(1,
                  "loca_pcp",
                  sprintf("%8.5f", spatial$lat),
                  sprintf("%9.5f", spatial$lon),
                  sprintf("%9.5f", spatial$elev),
                  sep=',')

        fileConn<-file(filename_hedp)
          writeLines(c(line_1,line_2), fileConn)
          close(fileConn)

        filename_pcp = str_c(out_dir,
                             "loca_pcp.txt",
                             sep = "")

        line_1 = "19500101"

                      

        fileConn<-file(filename_pcp)
          writeLines(c(line_1), fileConn)
          close(fileConn)
          
          subperiod_print = subset %>%
            select("pr")

          
          write_delim(x     = subperiod_print, 
                      file  =  filename_pcp, 
                      delim =       "",
                      append=TRUE) 



        #
        # Clean Up
        #

        remove(subset)
        remove(filename_tmp)
        remove(filename_pcp)
      } # scenario

    } # ensemble

  } # basin
  remove(ensemble)
  remove(scenario)
  remove(basin)



```