library(httr) library(jsonlite) library(dplyr) library(TTR) library(plotly) library(zoo) fetch_binance_batch <- function(symbol = "BTCUSDT", interval = "1d", start_date) { base_url <- "https://api.binance.com/api/v3/klines" start_ms <- as.numeric(as.POSIXct(start_date, tz = "UTC")) * 1000 res <- GET(base_url, query = list(symbol = symbol, interval = interval, startTime = format(start_ms, scientific = FALSE), limit = 1000)) stop_for_status(res) data_list <- fromJSON(content(res, "text", encoding = "UTF-8")) if (length(data_list) == 0) return(NULL) df <- as.data.frame(data_list, stringsAsFactors = FALSE) colnames(df) <- c("OpenTime", "Open", "High", "Low", "Close", "Volume", "CloseTime", "QuoteAssetVolume", "NumberOfTrades", "TakerBuyBaseAssetVolume", "TakerBuyQuoteAssetVolume", "Ignore") numeric_cols <- setdiff(names(df), "Ignore") df[numeric_cols] <- lapply(df[numeric_cols], as.numeric) df$OpenTime <- as.POSIXct(df$OpenTime / 1000, origin = "1970-01-01", tz = "UTC") df$CloseTime <- as.POSIXct(df$CloseTime / 1000, origin = "1970-01-01", tz = "UTC") df } fetch_binance_daily_full <- function(symbol = "BTCUSDT", interval = "1d", start_date) { all_data <- list() start_ms <- as.numeric(as.POSIXct(start_date, tz = "UTC")) * 1000 batch_num <- 1 repeat { cat(sprintf("Fetching batch %d from %s\n", batch_num, as.POSIXct(start_ms / 1000, origin = "1970-01-01", tz = "UTC"))) df <- fetch_binance_batch(symbol, interval, as.POSIXct(start_ms / 1000, origin = "1970-01-01", tz = "UTC")) if (is.null(df)) { cat("No more data.\n") break } all_data[[batch_num]] <- df last_open_ms <- as.numeric(df$OpenTime[nrow(df)]) * 1000 if (last_open_ms <= start_ms) break start_ms <- last_open_ms + 1 batch_num <- batch_num + 1 Sys.sleep(0.5) } do.call(rbind, all_data) } btc_daily <- fetch_binance_daily_full("BTCUSDT", "1d", "2017-01-01") btc_daily$SMA20 <- SMA(btc_daily$Close, n = 20) btc_daily$EMA20 <- EMA(btc_daily$Close, n = 20) btc_daily$RSI14 <- RSI(btc_daily$Close, n = 14) macd_vals <- MACD(btc_daily$Close, nFast = 12, nSlow = 26, nSig = 9) btc_daily$MACD <- macd_vals[,1] btc_daily$MACD_signal <- macd_vals[,2] btc_daily$MACD_hist <- macd_vals[,1] - macd_vals[,2] bbands <- BBands(btc_daily$Close, n = 20, sd = 2) btc_daily$BB_up <- bbands[, "up"] btc_daily$BB_dn <- bbands[, "dn"] btc_daily$BB_mavg <- bbands[, "mavg"] btc_daily$BB_pctB <- bbands[, "pctB"] adx_vals <- ADX(HLC = btc_daily[, c("High", "Low", "Close")], n = 14) btc_daily$ADX <- adx_vals[, "ADX"] btc_daily$DIp <- adx_vals[, "DIp"] btc_daily$DIn <- adx_vals[, "DIn"] stoch_vals <- stoch(HLC = btc_daily[, c("High", "Low", "Close")], nFastK = 14, nFastD = 3, nSlowD = 3) btc_daily$Stoch_K <- stoch_vals[, "fastK"] btc_daily$Stoch_D <- stoch_vals[, "fastD"] btc_daily$TrendSignal <- ifelse(lag(btc_daily$SMA20) < lag(btc_daily$EMA20) & btc_daily$SMA20 >= btc_daily$EMA20, "Golden Cross", ifelse(lag(btc_daily$SMA20) > lag(btc_daily$EMA20) & btc_daily$SMA20 <= btc_daily$EMA20, "Death Cross", NA)) btc_daily$RSI_signal <- ifelse(btc_daily$RSI14 > 70, "Overbought", ifelse(btc_daily$RSI14 < 30, "Oversold", NA)) btc_daily$VolumeSpike <- btc_daily$Volume > 2 * SMA(btc_daily$Volume, 20) btc_daily$BB_width <- btc_daily$BB_up - btc_daily$BB_dn threshold <- quantile(btc_daily$BB_width, 0.1, na.rm = TRUE) btc_daily$BB_squeeze <- btc_daily$BB_width < threshold btc_daily$log_return <- c(NA, diff(log(btc_daily$Close))) btc_daily$volatility_20d <- rollapply(btc_daily$log_return, 20, sd, fill = NA, align = "right") btc_daily$OpenTimeStr <- format(btc_daily$OpenTime, "%Y-%m-%d") text_info <- paste("Date:", btc_daily$OpenTimeStr, "<br>Close:", round(btc_daily$Close, 2), "<br>SMA20:", round(btc_daily$SMA20, 2), "<br>EMA20:", round(btc_daily$EMA20, 2), "<br>Signal:", btc_daily$TrendSignal) p1 <- plot_ly(btc_daily, x = ~OpenTimeStr, type = "candlestick", open = ~Open, close = ~Close, high = ~High, low = ~Low, name = "Candles", text = text_info, hoverinfo = "text") %>% add_lines(y = ~SMA20, name = "SMA 20", line = list(color = 'orange')) %>% add_lines(y = ~EMA20, name = "EMA 20", line = list(color = 'purple')) %>% layout(yaxis = list(title = "Price (USDT)")) p2 <- plot_ly(btc_daily, x = ~OpenTimeStr) %>% add_lines(y = ~MACD_hist, name = "MACD Hist", line = list(color = 'steelblue')) %>% add_lines(y = ~MACD_signal, name = "MACD Signal", line = list(color = 'orange')) %>% add_lines(y = ~MACD, name = "MACD", line = list(color = 'purple')) %>% layout(yaxis = list(title = "MACD")) p3 <- plot_ly(btc_daily, x = ~OpenTimeStr, y = ~RSI14, type = 'scatter', mode = 'lines', name = "RSI (14)", line = list(color = 'darkgreen')) %>% add_lines(y = rep(70, nrow(btc_daily)), name = "Overbought (70)", line = list(dash = 'dash', color = 'red')) %>% add_lines(y = rep(30, nrow(btc_daily)), name = "Oversold (30)", line = list(dash = 'dash', color = 'red')) %>% layout(yaxis = list(title = "RSI")) p4 <- plot_ly(btc_daily, x = ~OpenTimeStr, y = ~Volume, type = 'bar', name = "Volume", marker = list(color = ifelse(btc_daily$VolumeSpike, 'red', 'steelblue'))) %>% layout(yaxis = list(title = "Volume")) fig <- subplot(p1, p2, p3, p4, nrows = 4, shareX = TRUE, heights = c(0.4, 0.2, 0.2, 0.2), titleY = TRUE) %>% layout( xaxis4 = list(title = "Date", showgrid = FALSE, zeroline = FALSE), yaxis1 = list(title = "Price (USDT)"), yaxis2 = list(title = "MACD"), yaxis3 = list(title = "RSI"), yaxis4 = list(title = "Volume"), hovermode = "x unified", legend = list(orientation = "h", x = 0.1, y = 1.1), margin = list(t=40, b=60) ) fig

library(ncdf4) library(raster) library(sp) library(dplyr) library(sf) library(leaflet) library(openxlsx) # Set working directory setwd("c:/users/daniel/downloads") # Open the NetCDF files d18O <- nc_open("LGMR_d18Op_climo.nc") SAT <- nc_open("LGMR_SAT_climo.nc") # Extract dimensions and variables d18O_obsdatadates <- d18O$dim$age$vals lon <- ncvar_get(d18O, "lon") lat <- ncvar_get(d18O, "lat") d18O_time <- ncvar_get(d18O, "age") d18O_tmp_array <- ncvar_get(d18O, "d18Op") SAT_tmp_array <- ncvar_get(SAT, "sat") SAT_fillvalue <- ncatt_get(SAT, "sat", "_FillValue")$value # Fill value for missing data d18O_fillvalue <- ncatt_get(d18O, "d18Op", "_FillValue")$value # Create a data frame of coordinates from the NetCDF grid lonlat <- as.matrix(expand.grid(lon, lat)) # Define SP data frame SP <- data.frame( lon = as.numeric(c("47.55", "45.366", "41.848", "42.076")), lat = as.numeric(c("19.45", "22.873", "20.675", "20.896")) ) # Convert SP to SpatialPoints coordinates(SP) <- ~lat + lon proj4string(SP) <- CRS("+init=epsg:4326") # Prepare a data frame to store results results <- data.frame() # Loop through each time slice for (t in 1:length(d18O_obsdatadates)) { # Extract the slice for the current time step for both datasets d18O_nc_slice <- d18O_tmp_array[, , t] d18O_nc_slice[d18O_nc_slice == d18O_fillvalue] <- NA # Replace fill values with NA SAT_nc_slice <- SAT_tmp_array[, , t] SAT_nc_slice[SAT_nc_slice == SAT_fillvalue] <- NA # Replace fill values with NA # Convert the d18O slice to a data frame d18O_nc_vec <- as.vector(d18O_nc_slice) d18O_VP <- data.frame(cbind(lonlat, d18O_nc_vec)) names(d18O_VP) <- c("lon", "lat", "d18Op") d18O_VP <- na.omit(d18O_VP) # Convert the SAT slice to a data frame SAT_nc_vec <- as.vector(SAT_nc_slice) SAT_VP <- data.frame(cbind(lonlat, SAT_nc_vec)) names(SAT_VP) <- c("lon", "lat", "sat") SAT_VP <- na.omit(SAT_VP) # Combine d18O and SAT data frames combined_VP <- merge(d18O_VP, SAT_VP, by = c("lon", "lat")) # Convert combined_VP to SpatialPointsDataFrame coordinates(combined_VP) <- ~lon + lat proj4string(combined_VP) <- CRS("+init=epsg:4326") # Find the nearest neighbor for each point in SP nearest_indices <- apply(coordinates(SP), 1, function(pt) { distances <- spDistsN1(as.matrix(coordinates(combined_VP)), pt, longlat = TRUE) which.min(distances) }) # Extract the nearest values nearest_d18Op <- combined_VP@data[nearest_indices, "d18Op"] nearest_sat <- combined_VP@data[nearest_indices, "sat"] # Combine results with time and SP coordinates temp_results <- data.frame( lon = coordinates(SP)[, 1], lat = coordinates(SP)[, 2], time = d18O_obsdatadates[t], d18Op = nearest_d18Op, sat = nearest_sat ) # Append to the main results data frame results <- rbind(results, temp_results) } # Close the NetCDF files nc_close(d18O) nc_close(SAT) # Save results to a CSV file write.csv(results, "nearest_values.csv", row.names = FALSE) # Create a map plot for the first age slice d18O_first_age_slice <- d18O_tmp_array[, , 1] d18O_first_age_slice[d18O_first_age_slice == d18O_fillvalue] <- NA d18O_first_age_df <- data.frame(cbind(lonlat, as.vector(d18O_first_age_slice))) colnames(d18O_first_age_df) <- c("lon", "lat", "d18Op") d18O_first_age_df <- na.omit(d18O_first_age_df) # Use sf package for easier handling of spatial data # Convert d18O_first_age_df to sf object d18O_first_age_sf <- st_as_sf(d18O_first_age_df, coords = c("lon", "lat"), crs = 4326) sat_first_age_slice <- SAT_tmp_array[, , 1] sat_first_age_slice[sat_first_age_slice == SAT_fillvalue] <- NA sat_first_age_df <- data.frame(cbind(lonlat, as.vector(sat_first_age_slice))) colnames(sat_first_age_df) <- c("lon", "lat", "d18Op") sat_first_age_df <- na.omit(sat_first_age_df) # Use sf package for easier handling of spatial data # Convert d18O_first_age_df to sf object sat_first_age_sf <- st_as_sf(sat_first_age_df, coords = c("lon", "lat"), crs = 4326) # Convert SP to a data frame SP_df <- data.frame( lon = SP@coords[, 1], lat = SP@coords[, 2], extracted_d18Op = results$d18Op[results$time == d18O_obsdatadates[1]], # Add extracted d18Op values for the first age extracted_sat = results$sat[results$time == d18O_obsdatadates[1]] # Add extracted SAT values for the first age ) # Convert SP_df to sf object SP_sf <- st_as_sf(SP_df, coords = c("lon", "lat"), crs = 4326) # Create a base plot using leaflet leaflet() %>% addTiles() %>% addCircleMarkers( data = d18O_first_age_sf %>% st_coordinates() %>% as.data.frame(), lng = ~X, lat = ~Y, radius = 4, color = "blue", fillOpacity = 0.6, label = ~paste("d18Op:", round(d18O_first_age_df$d18Op, 2), round(sat_first_age_df$d18Op, 2)) ) %>% addCircleMarkers( data = SP_sf %>% st_coordinates() %>% as.data.frame(), lng = ~X, lat = ~Y, radius = 6, color = "red", fillOpacity = 1, label = ~paste("SP Point d18Op:", round(SP_df$extracted_d18Op, 2), "SAT:", round(SP_df$extracted_sat, 2)) ) # Save results to a CSV file write.xlsx(results, "netCDF_extract_Zoli.xlsx", rowNames = FALSE)

#Distance between the stations: 105114.58 basel_row <- select_validalo_aggregate[select_validalo_aggregate$Site == "BASEL", ] rovaniemi_row <- select_validalo_aggregate[select_validalo_aggregate$Site == "ROVANIEMI", ] svurban_row <- select_validalo_aggregate_sv[select_validalo_aggregate_sv$Site == "Sv. Urban", ] #random_integer1 = 26 min_distance = 49 random_integer <- sample(1:10000, 1) random_integer1 = 6097 set.seed(random_integer1) # Calculate the number of rows for the remaining 10% sampling desired_sample_size <- round(0.1 * nrow(select_validalo_aggregate)) -3 # Subtract 3 for Basel, Rovaniemi and Sv. Urban # Randomly select the remaining rows remaining_rows <- select_validalo_aggregate[!(select_validalo_aggregate$Site %in% c("BASEL", "ROVANIEMI","Sv. Urban")), ] random_indices <- sample(nrow(remaining_rows), size = desired_sample_size, replace = FALSE) random_selection <- remaining_rows[random_indices, ] # Combine the selected rows with Basel and Rovaniemi rows final_selection <- rbind(basel_row, rovaniemi_row,svurban_row, random_selection) final_selection_sp = final_selection coordinates(final_selection_sp) <- ~Longitude + Latitude proj4string(final_selection_sp) <- CRS("+init=epsg:4326") # Assuming final_selection_sp is your SpatialPointsDataFrame # Extract the coordinates from the SpatialPointsDataFrame coords <- coordinates(final_selection_sp) # Calculate the distance matrix in meters distance_matrix <- distm(coords, fun = distVincentyEllipsoid) # If you want to find the minimum non-zero distance min_distance <- min(distance_matrix[distance_matrix > 0]) Stations = select_validalo_aggregate_sp Validation = final_selection_sp mapview(Stations, alpha = 0.1) + mapview(Validation,zcol = "Group")

library(plotly) Prediction_map_3D = Prediction_map[!is.na( Prediction_map$Prediction), ] plot_ly(Prediction_map_3D, x = ~X, y = ~Y, z = ~Prediction, type = "scatter3d", mode = "markers", marker = list(size = 3, opacity = 0.8, color = ~Prediction, colorscale = "Viridis")) %>% layout(scene = list(xaxis = list(title = "Latitude", showgrid = TRUE, gridwidth = 10), yaxis = list(title = "Longitude", showgrid = TRUE, gridwidth = 10), zaxis = list(title = "d18O"),tickvals = rev(seq(min(Prediction_map_3D$Prediction), max(Prediction_map_3D$Prediction), length.out = 5)),camera = list(eye = list(x = 0, y = -0.1, z = 2))), showlegend = FALSE)

rf

Isoscape - 2022 December using Machine Learning

afasad

EU Isoscape using ML for the August monthly data of 2008.

data <- read.csv("data_ist.csv", header=TRUE, sep = ";" ,stringsAsFactors=FALSE) colnames(data) = c("SN","Y","X") data_sp = data coordinates(data_sp) <- ~ Y + X proj4string(data_sp) <- CRS("+init=epsg:23700") data_lat_lon = spTransform(data_sp, CRSobj = CRS("+init=epsg:4326")) data_cord_converted = data.frame(data,coordinates(data_lat_lon)) colnames(data_cord_converted) = c("SN","Y","X","Longitude","Latidue") data_mapview = data.frame(data_cord_converted, coordinates(data_lat_lon)) data_mapview$ID = row.names(data_mapview) data_mapview$category=c("Lila") colnames(data_mapview) = c("SN","Y","X","Longitude","Latidue","lon","lat","ID","category") data_mapview = data_mapview [1:132,] data2 <- read_excel("istvannak.xls") colnames(data2) = c( "Település","Hrsz./Telep","száma","Kataszteri száma", "Csőperem tengerszint feletti magassága", "Y" ,"X","Talpmélység","Szűrőzés -tól/-től [m]","Szűrőzés -ig [m]" ) data_sp2 = data2 coordinates(data_sp2) <- ~ Y + X proj4string(data_sp2) <- CRS("+init=epsg:23700") data_lat_lon2 = spTransform(data_sp2, CRSobj = CRS("+init=epsg:4326")) data_cord_converted2 = data.frame(data2,coordinates(data_lat_lon2)) colnames(data_cord_converted2) = c( "Település","Hrsz./Telep","száma","Kataszteri száma", "Csőperem tengerszint feletti magassága", "Y" ,"X","Talpmélység","Szűrőzés -tól/-től [m]","Szűrőzés -ig [m]" ) data_mapview2 = data.frame(data_cord_converted2, coordinates(data_lat_lon2)) colnames(data_mapview2) = c( "Település","Hrsz./Telep","száma","Kataszteri száma", "Csőperem tengerszint feletti magassága", "Y" ,"X","Talpmélység","Szűrőzés -tól/-től [m]","Szűrőzés -ig [m]","Longitude","Latidue","lon","lat" ) data_mapview2$category = c("Sárga") data_mapview2$SN = row.names(data_mapview2) database = data.frame(matrix(ncol = ncol(data_mapview2), nrow = nrow(data_mapview))) colnames(database) =c( "Település","Hrsz./Telep","száma","Kataszteri száma", "Csőperem tengerszint feletti magassága", "Y" ,"X","Talpmélység","Szűrőzés -tól/-től [m]","Szűrőzés -ig [m]","Longitude","Latidue","lon","lat","category","SN" ) database$SN = data_mapview$SN database$Y = data_mapview$Y database$X = data_mapview$X database$Longitude = data_mapview$Longitude database$Latidue = data_mapview$Latidue database$lon = data_mapview$lon database$lat = data_mapview$lat database$category = data_mapview$category database_sp = rbind(database, data_mapview2) coordinates(database_sp) <- ~ lon + lat proj4string(database_sp) <- CRS("+init=epsg:4326") mapview(database_sp, zcol = "category", alpha.regions = 1,label="SN")

map

EU Isoscape using ML for the January monthly data of 2008.

2 Random forest prediction of Stable oxygen isotopes from precipitation over europe, a review study for GeoMATES 2022 converence

Absolute differences of the measured and predicted values of δ18O using the different methods on Dataset 1 and Dataset 2. The boxes indicate the interquartile interval (IQI), the whiskers extend up from the top of the box to the largest value less than or equal to 1.5 times the IQI and down from the bottom of the box to the smallest value that is larger than 1.5 times the IQI. Values outside this range are considered to be outliers and are indicated by dots. The × shows the mean and the horizontal black line the median {Kovács, 2012 #579}. RK: regression kriging, RERF: Regression Enhanced Random Forest, and MRRF: multiple regression random forest. CODE: all(sapply( c("sp", "gstat","magrittr", "dplyr","xlsx","randomForestSRC"), require, character.only = T)) setwd("C:/Users/Daniel/Google Drive/Erdélyi Dániel/Doktori/UNKP/results/Figs") #change it Data <-readxl::read_excel("fig3_data.xlsx",col_names=TRUE, na="NA", sheet="Sheet1") library(tidyverse) library(hrbrthemes) library(viridis) Data$Model <- factor(Data$Model , levels = c('RK','MRRF','RERF','RF'),ordered = TRUE) Fig3 = Data %>% ggplot( aes(x=Dataset, y=Differences, fill=Model, col=Model)) + geom_boxplot(outlier.colour="black", notch=TRUE, lwd=0.5) + scale_fill_manual(values = c("#00B0F0", "#002060", "#92D050", "#FFC000")) + scale_color_manual(values = c("#005D7E", "#00B0F0", "#476D1D", "#9A7500")) + geom_jitter(color="black", size=0.5, alpha=0.3) + theme_ipsum() + theme( legend.position="right", plot.title = element_text(size=11) ) + ggtitle("Absolute Differences of the measured and predicted values of d18O") + xlab("Period") Fig3

Spatial distribution of the monitoring network of the Geostatistical- and Machine Learning modell comparison of predicting the spatial distribution of stable isotopes from precipitation over the British Isles and Germany

Top: 23 January Mid: 24 January Down: 25 January

Az adatok 10%-a random leválogatva, majd leellenőrizve, hogy nincsenek egymáshoz közelebb, mint 50 km. Script: all(sapply( c("parallel","doParallel","foreach","caret","pipeR","tidyr","glmnet","randomForestSRC","raster","magrittr","sp","mapview", "elevatr", "sf", "readxl", "GSIF", "rnaturalearth","rnaturalearthhires","rgeos","geosphere"), require, character.only = T)) setwd("C:/Users/Daniel/Google Drive/Erdélyi Dániel/Doktori/UNKP") #change it full_data_own_to_change_0123 <-readxl::read_excel("input_data/filtered/Data_0123.xlsx",col_names=TRUE, na="NA", sheet="Data_0123") #circles around points circle <- full_data_own_to_change_0123 coordinates(circle) <- ~longitude + latitude projection(circle) <- "+init=epsg:4326" circle2 <- spTransform(circle, CRSobj = CRS("+init=epsg:3857")) circle3 = data.frame( coordinates(circle2)) colnames(circle3)<-c("x", "y") circles <- cbind(circle, circle3) circles <- as.data.frame(circles) data <- circles coordinates(circles)<-~x + y proj4string(circles) <- CRS("+init=epsg:3857") circle_buffer <- gBuffer(circles, width = 50000) mapview(circle_buffer) + mapview(circles, alpha = 1) # select 10% of the data randomly, then check it on the map if they between 50km distance to each other datavalidation_set = data %>% dplyr::sample_frac(.1) coordinates(datavalidation_set) <- ~x + y projection(datavalidation_set) <- "+init=epsg:3857" all_stations <- circles mapview(circle_buffer, alpha.regions = 0, alpha =1) + mapview(datavalidation_set, zcol="site",col.regions = mapviewPalette("mapviewRasterColors"), alpha.regions = 1, alpha=1 ) + mapview(all_stations, alpha.regions = 0.2)

RERF test on the UK data jan 23. daily d18O from precipitation.