-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathIsotropic_variograms.R
More file actions
107 lines (84 loc) · 4.51 KB
/
Copy pathIsotropic_variograms.R
File metadata and controls
107 lines (84 loc) · 4.51 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
# Load necessary libraries
library(sf)
library(sp) # For converting to sp objects
library(gstat) # For variogram computation and modeling
library(ggplot2) # For additional plotting
# Define file paths
base_dir <- normalizePath("../") # Adjust base directory path as needed
# Update with your points shapefile path
# observation_shapefile <- file.path(base_dir, "kriging/San_Antonio/shapefiles_discharge", "discharge_sn_20160603_04_00.shp")
# observation_shapefile <- file.path(base_dir, "kriging/San_Antonio/shapefiles_discharge", "discharge_sn_20150522_12_00.shp")
# observation_shapefile <- file.path(base_dir, "kriging/San_Antonio/shapefiles_discharge", "discharge_sn_20150524_09_00.shp")
# observation_shapefile <- file.path(base_dir, "kriging/Guadalupe/guad_nexus/shapefiles_discharge", "discharge_gd_20181017_21_00.shp")
# observation_shapefile <- file.path(base_dir, "kriging/Guadalupe/guad_nexus/shapefiles_discharge", "discharge_gd_20181016_23_00.shp") #guadalupe_storm
# observation_shapefile <- file.path(base_dir, "kriging/Llano/llano_nexus/shapefiles_discharge", "discharge_ln_20181016_21_00.shp") #llano_storm_paper
observation_shapefile <- file.path(base_dir, "kriging/Llano/llano_nexus/shapefiles_discharge", "discharge_ln_20181017_22_00.shp") #20180926_06_00 nonstorm
# Read the shapefile (POINTS shapefile)
observation_points <- st_read(observation_shapefile)
# Convert sf object to sp object (gstat requires Spatial format)
observation_points_sp <- as(observation_points, "Spatial")
# Choose an attribute for the variogram calculation
# Replace 'your_attribute' with the actual attribute name in the shapefile
observation_points_sp$obs <- observation_points_sp$normalized
# Compute empirical variogram using gstat
# empirical_variogram <- variogram(obs ~ 1, observation_points_sp, cutoff = 75000, width = 5000) # Lag width of 3 km
empirical_variogram <- variogram(obs ~ 1, observation_points_sp, cutoff = 100000, width = 5000) # Lag width of 3 km
# Fit a theoretical variogram model (Exponential model with automatic fitting)
fitted_variogram <- fit.variogram(empirical_variogram, model = vgm(model = "Exp", nugget = 0.01), fit.sills = TRUE)
# Print the fitted variogram parameters
print(fitted_variogram)
# Extract the range from the fitted variogram
range_value <- fitted_variogram[fitted_variogram$model == "Exp", "range"]
# Convert empirical variogram to a data frame
variogram_df <- as.data.frame(empirical_variogram)
# Extract fitted variogram points for plotting
fitted_values <- variogramLine(fitted_variogram, maxdist = max(empirical_variogram$dist))
fitted_df <- data.frame(dist = fitted_values$dist, gamma = fitted_values$gamma)
library(ggplot2)
library(grid)
# Get max x and y values in the same units used in the plot
max_dist <- max(c(variogram_df$dist, fitted_df$dist)) # meters
# max_gamma <- max(c(variogram_df$gamma, fitted_df$gamma)) * 10000 for sn nonstorm, 100 for sn storm
max_gamma <- max(c(variogram_df$gamma, fitted_df$gamma)) *10
# Position for annotation
label_x <- max_dist * 0.01 / 1000 # km
label_y <- max_gamma * 0.98
variogram_plot <- ggplot() +
geom_point(data = variogram_df, aes(x = dist / 1000, y = gamma*100, color = "Empirical"), size = 3, alpha = 0.6) +
geom_line(data = fitted_df, aes(x = dist / 1000, y = gamma*100, color = "Fitted"), linewidth = 1.2) +
# Annotated panel letter
# annotate("text", x = label_x, y = label_y,
# label = "b", hjust = 0, vjust = 1.2, size = 5) +
labs(
x = "Distance (km)",
y = expression(Semivariance~(~10^-2)),
color = NULL
) +
theme_classic() +
theme(
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text = element_text(size = 16), # Match second plot
panel.border = element_rect(color = "black", fill = NA, linewidth = 1.2),
# legend.position = "none", # Match second plot
legend.position = c(0.88, 0.18),
legend.justification = c(1, 0),
# legend.background = element_rect(fill = "white", color = "black"),
legend.text = element_text(size = 14),
plot.margin = margin(t = 20, r = 10, b = 10, l = 10)
) +
scale_color_manual(
values = c("Empirical" = "darkgreen", "Fitted" = "red")
)
print(variogram_plot)
# Save using same size as Figure4b
# ggsave(
# # filename = file.path(base_dir, "Results_paper/paper_final/Figure3b_iso_test.pdf"),
# filename = file.path(base_dir, "Results_paper/paper_revision/llano_isotropic_receding560.pdf"),
# plot = variogram_plot,
# height = 3.07,
# width = 3.93,
# units = "in",
# dpi = 1000,
# scale = 1
# )