############### PART5: Extraction of temporal profile #############

#### Extract time series from raster stack

### Add dates to mkae it a date object?

#-72,42.24 #near Spencer MA

#48.281724, -71.401815 near Saguenay Quebec

#39.805293, -89.872626 near Springfield Illinois

#32.009676, -106.990266 near El Paso Us Texas

#39.955529, -105.263851: near Boulder Colorado

#45.384051, -121.859146 : near lost Creek near Mount Hood Oregon

#53.283685, -113.341702: leduc Canada near Edmonton

#39.009052, -77.026875: Silver Spring Sligo Creek, MD

#36.627806, -119.928901 : Raisin City, Fresno CA

#36.677139, -115.330122: Las Vegas

#35.696143, -78.396011: near Raleigh NC

list_lat_long <- list(

c( -72, 42.24),

c( -71.401815, 48.281724),

c( -89.872626, 39.805293),

c( -106.990266, 32.009676),

c( -105.263851, 39.955529), 

c( -121.859146,45.384051),

c( -113.341702,53.283685),

c( -77.026875,39.009052),

c( -119.928901,36.627806),

c(  -115.330122,36.677139),

c( -78.396011,35.696143)) 

query_for_station_lat_long <- function(point_val,df_points_spdf,step_x=0.25,step_y=0.25){

  #make this function better using gbuffer later!!!, works for now 

  #Improve: circular query + random within it

  y_val <- point_val[2]

  x_val <- point_val[1]

  y_val_tmp <- y_val + step_y

  if(x_val>0){

    x_val_tmp <- x_val - step_x

  }else{

    x_val_tmp <- x_val + step_x

  }

  test_df <- subset(df_points_spdf,(y_val < df_points_spdf$y & df_points_spdf$y < y_val_tmp))

  test_df2 <- subset(test_df,(x_val < test_df$x & test_df$x < x_val_tmp))

  #plot(test_df2)

 if(nrow(test_df2)>0){

   df_selected <- test_df2[1,]

 }else{

   df_selected <- NULL

 }

 return(df_selected)

}

test_day_query2 <- lapply(list_lat_long,FUN=query_for_station_lat_long,df_points_spdf=df_point_day,step_x=1,step_y=1)

#test_day_query <-query_for_station_lat_long(c(-72,42),df_points_spdf=df_point_day,step_x=1,step_y=0.25)

df_stations_selected <- do.call(rbind,test_day_query2)

proj4string(df_stations_selected) <- proj_str

#debug(query_for_station_lat_long)

##Next use the day results and select a mean station, quartile and min and max?

#list_id_data_v <- unique(data_stations_var_pred_data_v$id)

#list_id_data_s <- unique(data_stations_var_pred_data_s$id)

#Started at 4pm: on sept 9

list_id_data_v <- df_stations_selected$id

list_id_data_s <- df_stations_selected$id

write.table(data_v_subset,"data_v_subset_test.txt")

write.table(data_s_subset,"data_s_subset_test.txt")

##finish at 16.10pm on 09/09

#16.30 pm on 09/09

#data_stations_var_pred <- aggregate(id2 + date ~ x + y + dailyTmax + mod1 + res_mod1 ,data = data_stations, FUN=mean ) #+ mod1 + res_mod1 , data = data_stations, min)

dim(data_stations_var_pred)

#md <- melt(mydata, id=(c("id", "time")),)

md <- melt(data_stations, id=(c("id", "date")),measure.vars=c("x","y","dailyTmax","mod1","res_mod1"))

#formula_str <- "id + date ~ x + y + dailyTmax + mod1 + res_mod1"

data_stations_var_pred <- cast(md, id + date ~ variable, fun.aggregate = mean, 

  na.rm = TRUE)

write.table(data_stations_var_pred,

            file=file.path(out_dir,paste0("data_stations_var_pred_tmp_",out_suffix,".txt",

                                                                 sep=",")))

write.table(data_stations_var_pred_training_testing,"data_stations_var_pred_training_testing.txt")

md <- melt(data_stations, id=(c("id", "date")),measure.vars=c("training","testing"))

data_stations_training_testing <- cast(md, id + date ~ variable, fun.aggregate = sum, 

  na.rm = TRUE)

write.table(data_stations_training_testing,

            file=file.path(out_dir,paste0("data_stations_training_testing_",out_suffix,".txt",

                                                                 sep=",")))

#data_stations_var_pred <- aggregate(id2 ~ x + y + date + dailyTmax + mod1 + res_mod1 ,data = data_stations, mean ) #+ mod1 + res_mod1 , data = data_stations, min)

#data_stations$id2 <- as.numeric(data_stations$id)

#data_stations$date <- as.character(data_stations$date)

#[1] 57154     7

unique(data_stations_var_pred$id)

dim(data_stations_training_testing)

#[1] 57154     4

#mask_pred <- FALSE

#matching <- FALSE #to be added after mask_pred option

#list_param_pre_process <- list(raster_name_lf,python_bin,infile_mask,scaling,mask_pred,NA_flag_val,out_suffix,out_dir) 

#names(list_param_pre_process) <- c("lf","python_bin","infile_mask","scaling","mask_pred","NA_flag_val","out_suffix","out_dir") 

#lf_mosaic_scaled <- mclapply(1:length(raster_name_lf),FUN=pre_process_raster_mosaic_fun,list_param=list_param_pre_process,mc.preschedule=FALSE,mc.cores = num_cores)                         

##Now grab year year 1992 or matching year...maybe put this in a data.frame??

if(is.null(r_mosaic_fname)){

  pattern_str <-"*.tif"

  lf_mosaic_list <- list.files(path=in_dir_mosaic,pattern=pattern_str,recursive=F,full.names=T)

  r_mosaic <- stack(lf_mosaic_list,quick=T)

  save(r_mosaic,file="r_mosaic.RData")

}else{

  r_mosaic <- load_obj(r_mosaic_fname) #load raster stack of images

}

date_end <- "19991231"

df_points_day <- data_stations_var_pred_data_s

  #15.17 (on 09/08)

  #22.27 (on 09/08)

max_idst <- 0.009*5 #5km in degree

min_dist <- 0    #minimum distance to start with

###this needs be modified

## target number

target_max_nb <- 1 #this is not actually used yet in the current implementation

target_min_nb <- 5 #this is the target number of stations we would like

max_dist <- 1000 # the maximum distance used for pruning ie removes stations that are closer than 1000m 

min_dist <- 0    #minimum distance to start with

step_dist <- 1000 #iteration step to remove the stations

target_range_nb <- c(target_min_nb,target_max_nb) #target range of number of stations

#First increase distance till 5km

#then use random sampling...to get the extact target?

#First test with maximum distance of 100m

test1 <- sub_sampling_by_dist(target_range_nb,dist=min_dist,max_dist=max_dist,step=step_dist,data_in=df_points_day)