Bayesian Spatiotemporal Stochastic Partial Differential Equation for High-Resolution Earthquake Magnitude Mapping: Application to Sumatra Island, Indonesia
2024-11-23
Authors:
Dila Fitriani Azuri, I Gede Nyoman Mindra Jaya, Yudi Rosandi
Abstract:
Indonesia is highly prone to earthquakes, situated at the intersection
of three major tectonic plates: the Eurasian, Indo-Australian, and
Pacific plates. Sumatra Island, in particular, is one of the most
seismically active regions, lying within the megathrust zone where large
earthquakes frequently occur. The consequences of these earthquakes
often include significant loss of life and damage to infrastructure.
Therefore, proactive measures are essential to mitigate the devastating
effects of such disasters. One of the key preventive strategies is the
development of an Early Warning System (EWS), which can help predict
potential earthquakes based on seismic activity in specific locations.
However, not all areas have recorded earthquake events, making it
necessary to estimate seismic activity in unmeasured regions. This can
be achieved using the spatial interpolation technique, combined with the
Gaussian Field (GF) approach. Despite its accuracy, the use of GF
becomes computationally challenging due to the ‘big n problem,’ where
large datasets make computations more complex. To address this, the
Gaussian Markov Random Field (GMRF) method is applied. The GMRF
simplifies computations by breaking down continuous spatial data into
discrete components using Stochastic Partial Differential Equations
(SPDE). SPDE improves spatial modelling efficiency via triangulation and
GMRF transformation. The model results indicate that the optimal
parameters are a 10 km range and a standard deviation of 0.1. By
identifying zones with earthquake potential through SPDE spatiotemporal
modeling and high-resolution mapping, it provides essential insights for
policymakers to establish seismic safety regulations and design public
awareness campaigns, ultimately enhancing community resilience and
improving emergency response planning.
Data Sources:
There are several data sources for this research:
- Earthquake data: https://www.isc.ac.uk/isc-ehb/search/catalogue/
- Indonesia maps: https://gadm.org/download_country36.html
- Indonesia’s Active Faults: https://geoaccess.id/dataspasial/
- District and Sub-district Level: https://www.indonesia-geospasial.com/2020/04/download-shapefile-shp-batas.html
Steps:
1. Install and Read the Library
#set working directory
setwd("/Users/dila.azuri/Documents/Magister/Tesis/SPDE/Syntax")
#read libraries
library("devtools")## Loading required package: usethis## Warning: package 'usethis' was built under R version 4.3.3library(ggplot2)
library(sf)## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUElibrary(gstat)
library(sp)
library(INLA) #INLA## Loading required package: Matrix## This is INLA_23.09.09 built 2023-10-16 17:38:01 UTC.
## - See www.r-inla.org/contact-us for how to get help.library(writexl) #save to xls
library(dplyr) #for pivot in statistics descriptive##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(plotly) # for surface plot##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(akima) #for smoothness interpolation result
library(cowplot) #gabungin plot
library(patchwork) #label title di ggplot## Warning: package 'patchwork' was built under R version 4.3.3##
## Attaching package: 'patchwork'## The following object is masked from 'package:cowplot':
##
## align_plots2. Input Data
2.1 Input Data Earthquake
#====================== Earthquake data
Earthquakes<-read.csv("INLA - Gempa.csv",sep=";")
head(Earthquakes)## EVENTID TYPE AUTHOR DATE Tahun Bulan TIME LAT LON
## 1 610793744 ke ISC-EHB 2000-01-06 2000 1 00:56:22.27 2,28 98,18
## 2 610794108 ke ISC-EHB 2000-02-21 2000 2 02:57:22.42 -5,826 106,228
## 3 610794223 ke ISC-EHB 2000-03-07 2000 3 18:49:15.42 3,623 95,963
## 4 610794245 ke ISC-EHB 2000-03-10 2000 3 21:32:11.69 4,746 95,908
## 5 610794275 ke ISC-EHB 2000-03-14 2000 3 22:10:24.95 -4,839 103,32
## 6 610794483 ke ISC-EHB 2000-04-03 2000 4 19:02:30.95 2,507 97,153
## DEPTH DEPFIX DEPQUAL AUTHOR.1 TYPE.1 MAG
## 1 60 TRUE L3 ISC mb 5
## 2 140 TRUE L3 ISC mb 4,5
## 3 40 TRUE L3 ISC mb 4,6
## 4 15,5 TRUE L2 ISC mb 5
## 5 50,9 NA L2 ISC mb 5
## 6 30 TRUE L3 ISC mb 4,7#======= Choose Data: 2020
Earthquakes20=Earthquakes[Earthquakes$Tahun==2020,]
head(Earthquakes20)## EVENTID TYPE AUTHOR DATE Tahun Bulan TIME LAT LON
## 3496 626319012 se ISC-EHB 2020-01-07 2020 1 06:05:20.15 2,435 96,373
## 3497 626319040 se ISC-EHB 2020-01-07 2020 1 22:48:03.27 -1,566 99,37
## 3498 626319145 se ISC-EHB 2020-01-12 2020 1 10:45:18.14 2,98 97,834
## 3499 626319317 se ISC-EHB 2020-01-21 2020 1 11:03:07.82 0,668 99,819
## 3500 626319329 se ISC-EHB 2020-01-22 2020 1 03:58:45.43 2,213 96,265
## 3501 626319449 se ISC-EHB 2020-01-27 2020 1 03:32:36.59 -5,145 102,408
## DEPTH DEPFIX DEPQUAL AUTHOR.1 TYPE.1 MAG
## 3496 21,6 NA L1 ISC mb 6
## 3497 28,5 NA L2 ISC mb 4,5
## 3498 66,8 NA L1 ISC mb 4,8
## 3499 157,2 NA L1 ISC mb 4,5
## 3500 23,1 NA L2 ISC mb 4,7
## 3501 41,5 NA L1 ISC mb 4,7#======= Data Preprocessing
Earthquakes20$LAT <- gsub(",", ".", Earthquakes20$LAT)
Earthquakes20$LAT <- as.numeric(Earthquakes20$LAT)
Earthquakes20$LON <- gsub(",", ".", Earthquakes20$LON)
Earthquakes20$LON <- as.numeric(Earthquakes20$LON)
Earthquakes20$MAG <- gsub(",", ".", Earthquakes20$MAG)
Earthquakes20$MAG <- as.numeric(Earthquakes20$MAG)
colnames(Earthquakes20)[colnames(Earthquakes20) == "Bulan"] <- "Month"
head(Earthquakes20)## EVENTID TYPE AUTHOR DATE Tahun Month TIME LAT LON
## 3496 626319012 se ISC-EHB 2020-01-07 2020 1 06:05:20.15 2.435 96.373
## 3497 626319040 se ISC-EHB 2020-01-07 2020 1 22:48:03.27 -1.566 99.370
## 3498 626319145 se ISC-EHB 2020-01-12 2020 1 10:45:18.14 2.980 97.834
## 3499 626319317 se ISC-EHB 2020-01-21 2020 1 11:03:07.82 0.668 99.819
## 3500 626319329 se ISC-EHB 2020-01-22 2020 1 03:58:45.43 2.213 96.265
## 3501 626319449 se ISC-EHB 2020-01-27 2020 1 03:32:36.59 -5.145 102.408
## DEPTH DEPFIX DEPQUAL AUTHOR.1 TYPE.1 MAG
## 3496 21,6 NA L1 ISC mb 6.0
## 3497 28,5 NA L2 ISC mb 4.5
## 3498 66,8 NA L1 ISC mb 4.8
## 3499 157,2 NA L1 ISC mb 4.5
## 3500 23,1 NA L2 ISC mb 4.7
## 3501 41,5 NA L1 ISC mb 4.72.1.1 Exploratory Data Analysis
2.1.1.1 Statistics Descriptive
# Descriptive statistics: Monthly
mag_stats_by_month <- Earthquakes20 %>%
group_by(Month) %>%
summarise(
count = n(),
mean_mag = mean(MAG, na.rm = TRUE),
median_mag = median(MAG, na.rm = TRUE),
sd_mag = sd(MAG, na.rm = TRUE),
min_mag = min(MAG, na.rm = TRUE),
max_mag = max(MAG, na.rm = TRUE)
)
mag_stats_by_month## # A tibble: 12 × 7
## Month count mean_mag median_mag sd_mag min_mag max_mag
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 8 4.82 4.7 0.512 4.4 6
## 2 2 8 4.72 4.7 0.358 4.2 5.3
## 3 3 6 4.77 4.9 0.528 4.1 5.5
## 4 4 7 4.7 4.7 0.365 4.1 5.2
## 5 5 3 4.67 4.7 0.351 4.3 5
## 6 6 5 4.8 4.7 0.447 4.3 5.4
## 7 7 9 4.57 4.5 0.354 4.1 5.3
## 8 8 13 4.89 4.9 0.685 4 6.2
## 9 9 11 4.6 4.6 0.361 4.1 5.2
## 10 10 19 4.94 4.9 0.358 4.4 5.6
## 11 11 13 4.86 4.8 0.411 4.3 5.8
## 12 12 7 4.67 4.5 0.309 4.4 5.1#Descriptive statistics: overall
mean(Earthquakes20$MAG)
sd(Earthquakes20$MAG)
min(Earthquakes20$MAG)
max(Earthquakes20$MAG)2.1.1.2 Boxplot
# Outlier Detection: Boxplot
ggplot(Earthquakes20, aes(x = as.factor(Month), y = MAG, fill = as.factor(Month))) +
geom_boxplot(color = "black") +
labs(title = "Boxplot of Earthquake Magnitude by Month",
x = "Month",
y = "Magnitude (Mw)") +
scale_x_discrete(labels = c("Jan", "Feb", "Mar", "Apr", "May", "June",
"July", "Aug", "Sept", "Oct", "Nov", "Dec")) +
scale_fill_manual(values = c("1" = "cadetblue3", "2" = "#377EB8", "3" = "cornsilk1",
"4" = "#984EA3", "5" = "coral", "6" = "deeppink3",
"7" = "bisque", "8" = "#F781BF", "9" = "cornsilk4",
"10" = "aliceblue", "11" = "brown", "12" = "darkseagreen")) +
theme_minimal()+
theme(
panel.background = element_rect(fill = "white", color = NA),
panel.border = element_rect(color = "white", fill = NA, linewidth = 1),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid")
)
2.2 Input Sumatra Maps
The maps consist of all provinces in Indonesia. For this study, the focus is on Sumatra Island, which includes the provinces of Kepulauan Riau, Sumatera Selatan, Lampung, Sumatera Utara, Bangka Belitung, Jambi, Riau, Bengkulu, Sumatera Barat, and Aceh.
Indonesia<-readRDS('gadm36_IDN_1_sp.rds')
Indonesia$NAME_1
Sumatra<-Indonesia[Indonesia$NAME_1 %in%c("Kepulauan Riau","Sumatera Selatan","Lampung","Sumatera Utara","Bangka Belitung","Jambi","Riau","Bengkulu","Sumatera Barat","Aceh"),]
Sumatra_sf <- st_as_sf(Sumatra)2.3 Input Active Faults
Active faults are also incorporated into the mapping to determine whether the earthquake-prone areas are located near these faults.
fault = st_read(file.choose())
sumatera_fault <- st_bbox(c(xmin = 95, ymin = -6, xmax = 106, ymax = 6), crs = st_crs(fault))
fault_sumatera <- st_crop(fault, sumatera_fault)2.4 Epicenters Visualization
The map of the study area illustrates the distribution of relocated earthquake epicenters with magnitudes exceeding 4.0 in 2020. This map is presented as Figure 1 in the paper, with the blue line representing the fault.
2.4.1 Sumatra Maps
Month_names <- c( "1"="January",
"2"="February",
"3"="March",
"4"="April",
"5"="May",
"6"="June",
"7"="July",
"8"="August",
"9"="September",
"10"="October",
"11"="November",
"12"="December")
ggplot() +
geom_sf(data = Sumatra_sf) +
geom_sf(data = fault_sumatera, color = "blue", size = 1) +
geom_point(data = as.data.frame(Earthquakes20), aes(x = LON, y = LAT, color =MAG)) +
scale_color_gradient(low = "yellow", high = "red")+
labs(title = "Distribution of relocated earthquake epicenters", x = "Longitude", y = "Latitude") +
theme_bw()
2.4.2 Monthly
ggplot() +
geom_sf(data = Sumatra_sf) + # Menggunakan data = Sumatra_sf
geom_sf(data = patahan_sumatera, color = "bisque", size = 0.6) +
geom_point(data = as.data.frame(Earthquakes20), aes(x = LON, y = LAT, color =MAG),size=0.8) + # Mempetakan koordinat Earthquakes
facet_wrap(~Month, labeller = as_labeller(Month_names)) + # Membagi peta berdasarkan bulan
labs(title = "Peta Sumatra", x = "Longitude", y = "Latitude") +
scale_color_gradient(low = "yellow", high = "red")+
theme_bw() +
theme(axis.text = element_blank())
3. Data Preparation
In data analysis using a Bayesian Spatiotemporal model, it’s essential to have datasets with consistent dimensions across different time points. However, earthquake data often do not occur at the same locations across different times. Therefore, the initial step in conducting this analysis involves creating a spatiotemporal dataset with uniform dimensions.
To achieve this, we define a unique set of locations and repeat this set for each time period, such as monthly. For locations that are not recorded in the initial dataset or are not found at a given time, we assign a value of “NA.” Below is the detailed process for this procedure.
# Create a unique list of all Longitude and Latitude combinations
locations <- unique(Earthquakes20[, c("LON", "LAT")])
# Create a grid of all combinations of Longitude, Latitude, and Time
all_combinations <- merge(locations, data.frame(Month = unique(Earthquakes20$Month)))
head(all_combinations)
# Merge with original Earthquakes20
final_Earthquakes20 <- merge(all_combinations, Earthquakes20,
by = c("LON", "LAT", "Month"), all.x = TRUE)
head(final_Earthquakes20)
# Fill NA for Locations Missing Data at Specific Times
final_Earthquakes20$MAG[is.na(final_Earthquakes20$MAG)] <- NA
final_Earthquakes20$DEPTH[is.na(final_Earthquakes20$DEPTH)] <- NA
final_Earthquakes20$Year[is.na(final_Earthquakes20$Year)] <- 2020
# Construct final dataframe
Earthquakes20<-data.frame(LON=final_Earthquakes20$LON,LAT=final_Earthquakes20$LAT,
Month=final_Earthquakes20$Month, MAG=final_Earthquakes20$MAG,
DEPTH=final_Earthquakes20$DEPTH)
head(Earthquakes20)
ggplot() +
geom_sf(data = Sumatra_sf) + # Menggunakan data = Sumatra_sf
geom_sf(data = patahan_sumatera, color = "coral", size = 1) +
geom_point(data = as.data.frame(Earthquakes20), aes(x = LON, y = LAT, color =MAG),size=0.7) + # Mempetakan koordinat Earthquakes
facet_wrap(~Month, labeller = as_labeller(Month_names)) + # Membagi peta berdasarkan bulan
labs(title = "Peta Sumatra", x = "Longitude", y = "Latitude") +
scale_color_gradient(low = "yellow", high = "red", na.value = "darkseagreen")+
theme_bw() +
theme(axis.text = element_blank())
4. Grid
Once we have datasets with consistent spatial and temporal dimensions, the next step is to define the grid area for analysis. A challenge with earthquake data is the absence of predefined geographic boundaries. Therefore, we need to establish our own prediction area boundaries.
To address this, we create a grid that encompasses all earthquake points recorded during the observed time period, such as the year 2020. This approach ensures that our grid covers the entire area where earthquakes have occurred, providing a comprehensive framework for analysis. It is important to note that in spatial data analysis, coordinates must be transformed from longitude and latitude to Universal Transverse Mercator (UTM). This is necessary because the calculation of Euclidean distances relies on UTM coordinates, not on longitude and latitude. The transformation process is outlined in detail in the following steps.
An important consideration is determining the grid interval, as this will influence the accuracy of clustering calculations. Typically, finding the optimal interval involves a trial-and-error process until the desired clustering results are achieved. As a general guideline, if observation points are sparse, it is advisable to use a wider grid interval to capture the broader spatial distribution of the data.
4.1 Several Interval Grid
# Geographic coordinate data (longitude and latitude)
polygon_coords <- data.frame(
longitude = c(94, 94, 98, 108, 102, 94),
latitude = c(2, 6, 6, -7, -7, 2)
)
polygon_sf <- st_as_sf(polygon_coords, coords = c("longitude", "latitude"), crs = 4326)
# Convert to UTM (e.g., UTM Zone 48N)
polygon_utm <- st_transform(polygon_sf, crs = 32648)
utm_coords <- as.data.frame(st_coordinates(polygon_utm))
# Membuat data frame untuk utm_coords jika belum ada
utm_coords_df <- data.frame(longitude = utm_coords$X, latitude = utm_coords$Y)4.1.1 10km x 10km
#====================== 4.1 GRID: 10km x 10km
# Define the interval between grids. Update if needed
interval1 <- 10000
grid_points1 <- expand.grid(
longitude = seq(min(utm_coords$X), max(utm_coords$X), by = interval1),
latitude = seq(min(utm_coords$Y), max(utm_coords$Y), by = interval1)
)
inside1 <- point.in.polygon(grid_points1$longitude, grid_points1$latitude,
utm_coords$X, utm_coords$Y)
inside_points1 <- grid_points1[inside1 == 1, ]
# create grid plot
ggplot() +
geom_path(data = utm_coords_df, aes(x = longitude, y = latitude), color = "red") +
geom_point(data = inside_points1, aes(x = longitude, y = latitude), color = "blue", size = 0.5) +
labs(x = "Longitude", y = "Latitude", title = "Interval Grid: 10000") +
theme_minimal()
4.1.2 20km x 20km
#====================== 4.2 GRID: 20km x 20km
# Define the interval between grids. Update if needed
interval2 <- 20000
grid_points2 <- expand.grid(
longitude = seq(min(utm_coords$X), max(utm_coords$X), by = interval2),
latitude = seq(min(utm_coords$Y), max(utm_coords$Y), by = interval2)
)
inside2 <- point.in.polygon(grid_points2$longitude, grid_points2$latitude,
utm_coords$X, utm_coords$Y)
inside_points2 <- grid_points2[inside2 == 1, ]
# create grid plot
ggplot() +
geom_path(data = utm_coords_df, aes(x = longitude, y = latitude), color = "red") +
geom_point(data = inside_points2, aes(x = longitude, y = latitude), color = "blue", size = 0.5) +
labs(x = "Longitude", y = "Latitude", title = "Interval Grid: 20000") +
theme_minimal()
4.1.3 30km x 30km
#====================== 4.3 GRID v1: 30km x 30km
# Define the interval between grids. Update if needed
interval3 <- 30000
grid_points3 <- expand.grid(
longitude = seq(min(utm_coords$X), max(utm_coords$X), by = interval3),
latitude = seq(min(utm_coords$Y), max(utm_coords$Y), by = interval3)
)
inside3 <- point.in.polygon(grid_points3$longitude, grid_points3$latitude,
utm_coords$X, utm_coords$Y)
inside_points3 <- grid_points3[inside3 == 1, ]
# Membuat plot dengan ggplot2
ggplot() +
geom_path(data = utm_coords_df, aes(x = longitude, y = latitude), color = "red") +
geom_point(data = inside_points3, aes(x = longitude, y = latitude), color = "blue", size = 0.5) +
labs(x = "Longitude", y = "Latitude", title = "Interval Grid: 30000") +
theme_minimal()
4.1.4 40km x 40km
#====================== 4.4 GRID: 40km x 40km
# Define the interval between grids. Update if needed
interval4 <- 40000
grid_points4 <- expand.grid(
longitude = seq(min(utm_coords$X), max(utm_coords$X), by = interval4),
latitude = seq(min(utm_coords$Y), max(utm_coords$Y), by = interval4)
)
inside4 <- point.in.polygon(grid_points4$longitude, grid_points4$latitude,
utm_coords$X, utm_coords$Y)
inside_points4 <- grid_points4[inside4 == 1, ]
# create grid plot
ggplot() +
geom_path(data = utm_coords_df, aes(x = longitude, y = latitude), color = "red") +
geom_point(data = inside_points4, aes(x = longitude, y = latitude), color = "blue", size = 0.5) +
labs(x = "Longitude", y = "Latitude", title = "Interval Grid: 40000") +
theme_minimal()
4.1.5 50km x 50km
#====================== 4.5 GRID: 50km x 50km
# Define the interval between grids. Update if needed
interval5 <- 50000
grid_points5 <- expand.grid(
longitude = seq(min(utm_coords$X), max(utm_coords$X), by = interval5),
latitude = seq(min(utm_coords$Y), max(utm_coords$Y), by = interval5)
)
inside5 <- point.in.polygon(grid_points5$longitude, grid_points5$latitude,
utm_coords$X, utm_coords$Y)
inside_points5 <- grid_points5[inside5 == 1, ]
# create grid plot
ggplot() +
geom_path(data = utm_coords_df, aes(x = longitude, y = latitude), color = "red") +
geom_point(data = inside_points5, aes(x = longitude, y = latitude), color = "blue", size = 0.5) +
labs(x = "Longitude", y = "Latitude", title = "Interval Grid: 50000") +
theme_minimal()
4.2 Grid align with the observation
Another crucial step is to ensure that the grid aligns with the locations of the observation data. Here’s how we can verify this alignment.
4.2.1 10km x 10km
## Square grid with constant interval
coordinates(inside_points1)<-~longitude+latitude
gridded(inside_points1) = TRUE
## Get the UTM grid location
coordsGrid.UTM1<-inside_points1 ## Get the UTM grid location
UTM<-CRS("+proj=utm +zone=48 ellps=WGS84")
Sumatra_UTM <- spTransform(Sumatra, UTM)
coordsGrid.UTM_df1 <- as.data.frame(coordsGrid.UTM1)
Sumatra_UTM_sf <- st_as_sf(Sumatra_UTM)
ggplot() +
geom_sf(data = Sumatra_UTM_sf, fill = "white", color = "black") +
geom_point(data = coordsGrid.UTM_df1, aes(x = longitude, y = latitude), color = "blue", size = 0.5) +
theme_minimal() +
ggtitle("Sumatra Maps with Grid UTM")
4.2.2 20km x 20km
## Square grid with constant interval
coordinates(inside_points2)<-~longitude+latitude
gridded(inside_points2) = TRUE
coordsGrid.UTM2<-inside_points2 ## Get the UTM grid location
coordsGrid.UTM_df2 <- as.data.frame(coordsGrid.UTM2)
ggplot() +
geom_sf(data = Sumatra_UTM_sf, fill = "white", color = "black") +
geom_point(data = coordsGrid.UTM_df2, aes(x = longitude, y = latitude), color = "blue", size = 0.5) +
theme_minimal() +
ggtitle("Sumatra Maps with Grid UTM")
4.2.3 30km x 30km
## Square grid with constant interval
coordinates(inside_points3)<-~longitude+latitude
gridded(inside_points3) = TRUE
coordsGrid.UTM3<-inside_points3 ## Get the UTM grid location
coordsGrid.UTM_df3 <- as.data.frame(coordsGrid.UTM3)
ggplot() +
geom_sf(data = Sumatra_UTM_sf, fill = "white", color = "black") +
geom_point(data = coordsGrid.UTM_df3, aes(x = longitude, y = latitude), color = "blue", size = 0.5) +
theme_minimal() +
ggtitle("Sumatra Maps with Grid UTM")
4.2.4 40km x 40km
## Square grid with constant interval
coordinates(inside_points4)<-~longitude+latitude
gridded(inside_points4) = TRUE
coordsGrid.UTM4<-inside_points4 ## Get the UTM grid location
coordsGrid.UTM_df4 <- as.data.frame(coordsGrid.UTM4)
ggplot() +
geom_sf(data = Sumatra_UTM_sf, fill = "white", color = "black") +
geom_point(data = coordsGrid.UTM_df4, aes(x = longitude, y = latitude), color = "blue", size = 0.5) +
theme_minimal() +
ggtitle("Sumatra Maps with Grid UTM")
4.2.5 50km x 50km
## Square grid with constant interval
coordinates(inside_points5)<-~longitude+latitude
gridded(inside_points5) = TRUE
coordsGrid.UTM5<-inside_points5 ## Get the UTM grid location
coordsGrid.UTM_df5 <- as.data.frame(coordsGrid.UTM5)
ggplot() +
geom_sf(data = Sumatra_UTM_sf, fill = "white", color = "black") +
geom_point(data = coordsGrid.UTM_df5, aes(x = longitude, y = latitude), color = "blue", size = 0.5) +
theme_minimal() +
ggtitle("Sumatra Maps with Grid UTM")
5. Prepare the data in the UTM system
st_crs(Sumatra_sf)## Coordinate Reference System:
## User input: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## wkt:
## BOUNDCRS[
## SOURCECRS[
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Geocentric translations (geog2D domain)",
## ID["EPSG",9603]],
## PARAMETER["X-axis translation",0,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",0,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",0,
## ID["EPSG",8607]]]]Sumatra_sf_utm1 <- st_transform(Sumatra_sf, crs = 32648)
## Observed
Observed<-data.frame(X=Earthquakes20$LON,Y=Earthquakes20$LAT)
coordinates(Observed) <- ~X+Y # Asumsikan kolom X dan Y adalah kolom lon dan lat
proj4string(Observed) <- CRS("+proj=longlat +datum=WGS84") # Set CRS ke WGS84 (EPSG: 4326)
Observed_utm <- spTransform(Observed, CRS("+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"))
Observed_utm1<-coordinates(Observed_utm)
colnames(Observed_utm1)<-c("X","Y")
Earthquakes20$UTM_x<-Observed_utm1[,1]
Earthquakes20$UTM_y<-Observed_utm1[,2]
Earthquakes20$IT<-Earthquakes20$Month
head(Earthquakes20)6. Construct Model
6.1 Define the Triangulation Mesh
In this study, the mesh utilized is designed as shown. The construction of the mesh may involve a trial-and-error approach. Further details about the methodology can be found in this referenced journal: https://www.tandfonline.com/doi/full/10.1080/03610926.2018.1536209
Earthquakes_mesh <- inla.mesh.2d(
loc=cbind(Earthquakes20$UTM_x,Earthquakes20$UTM_y),
offset=c(50000, 200000),
max.edge=c(100000, 1000000))
plot(Earthquakes_mesh,col="gray80", main="Earthquakes Mesh of Sumatra")
points(Observed_utm1,col="red")
plot(Sumatra_UTM,add=T, border="blue")
6.2 Define the model for observed and grids areas
coordinates.allyear<-as.matrix(data.frame(x=Earthquakes20$UTM_x,y=Earthquakes20$UTM_y))
A_est <- inla.spde.make.A(mesh=Earthquakes_mesh,loc=coordinates.allyear,
group=Earthquakes20$Month, n.group=12)
dim(A_est)## [1] 1308 62886.3 Define the INLA parameters required for computation
#============ INLA
control <- list(
results = list(return.marginals.random = TRUE, return.marginals.predictor=TRUE),
compute = list(hyperpar=TRUE, return.marginals.predictor=TRUE, return.marginals=TRUE, dic=TRUE, mlik = TRUE, cpo = TRUE,
po = TRUE, waic=TRUE))6.4 Define prior distribution of Penalized Complexity
6.4.1 For 10km x 10km
Grid.CoordUTM1<-as.data.frame(coordsGrid.UTM1)
m1<-nrow(Grid.CoordUTM1)
Grid.CoordUTM1<-as.matrix(Grid.CoordUTM1)
GroupTime1<-rep(c(1:12),each=m1)
All.Grid1<-kronecker(matrix(1, 12, 1), Grid.CoordUTM1)6.4.1.1 Range: 10km and Standard Deviation 0.1
prior.median.sd1001 = 0.1; prior.median.range10 = 10000
Earthquakes_spde1001 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range10, .5),
prior.sigma = c(prior.median.sd1001, .5))
s_index1001 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde1001$n.spde,
n.group=12)
stack_est1001 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index1001) , tag="est")
A_pred1 <- inla.spde.make.A(mesh=Earthquakes_mesh,
loc=as.matrix(All.Grid1),
group=GroupTime1, #selected day for prediction
n.group=12)
Cov1<-data.frame(UTM_x=All.Grid1[,1]/1000,UTM_y=All.Grid1[,2]/1000,DEPTH=NA,Time=GroupTime1)
stack_pred1001 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred1),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid1)),
Time=GroupTime1, UTM_x=Cov1$UTM_x,
UTM_y=Cov1$UTM_y, DEPTH=Cov1$DEPTH), s_index1001),
tag="pred")
stack1001 <- inla.stack(stack_est1001, stack_pred1001) 6.4.1.2 Range: 10km and Standard Deviation 0.5
prior.median.sd1005 = 0.5; prior.median.range10 = 10000
Earthquakes_spde1005 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range10, .5),
prior.sigma = c(prior.median.sd1005, .5))
s_index1005 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde1005$n.spde,
n.group=12)
stack_est1005 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index1005) , tag="est")
stack_pred1005 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred1),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid1)),
Time=GroupTime1, UTM_x=Cov1$UTM_x,
UTM_y=Cov1$UTM_y, DEPTH=Cov1$DEPTH), s_index1005),
tag="pred")
stack1005 <- inla.stack(stack_est1005, stack_pred1005)6.4.1.3 Range: 10km and Standard Deviation 1
prior.median.sd101 = 1; prior.median.range10 = 10000
Earthquakes_spde101 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range10, .5),
prior.sigma = c(prior.median.sd101, .5))
s_index101 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde101$n.spde,
n.group=12)
stack_est101 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index101) , tag="est")
stack_pred101 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred1),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid1)),
Time=GroupTime1, UTM_x=Cov1$UTM_x,
UTM_y=Cov1$UTM_y, DEPTH=Cov1$DEPTH), s_index101),
tag="pred")
stack101 <- inla.stack(stack_est101, stack_pred101) 6.4.1.4 Range: 10km and Standard Deviation 1.5
prior.median.sd1015 = 1.5; prior.median.range10 = 10000
Earthquakes_spde1015 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range10, .5),
prior.sigma = c(prior.median.sd1015, .5))
s_index1015 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde1015$n.spde,
n.group=12)
stack_est1015 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index1015) , tag="est")
stack_pred1015 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred1),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid1)),
Time=GroupTime1, UTM_x=Cov1$UTM_x,
UTM_y=Cov1$UTM_y, DEPTH=Cov1$DEPTH), s_index1015),
tag="pred")
stack1015 <- inla.stack(stack_est1015, stack_pred1015) 6.4.1.5 Range: 10km and Standard Deviation 2
prior.median.sd102 = 2; prior.median.range10 = 10000
Earthquakes_spde102 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range10, .5),
prior.sigma = c(prior.median.sd102, .5))
s_index102 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde102$n.spde,
n.group=12)
stack_est102 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index102) , tag="est")
stack_pred102 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred1),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid1)),
Time=GroupTime1, UTM_x=Cov1$UTM_x,
UTM_y=Cov1$UTM_y, DEPTH=Cov1$DEPTH), s_index102),
tag="pred")
stack102 <- inla.stack(stack_est102, stack_pred102) 6.5 Run INLA program
It should be noted that this calculation takes quite a long time, based on the size of grid.
6.5.1 10km x 10km
6.5.1.1 Range: 10km and Standard Deviation 0.1
formula1001 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde1001,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output1001 <- inla(formula1001,
data=inla.stack.data(stack1001, spde=Earthquakes_spde1001),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack1001), compute=TRUE))
summary(output1001)Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles
= quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials,
strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ",
" lincomb = lincomb, selection = selection, control.compute = control.compute, ", "
control.predictor = control.predictor, control.family = control.family, ", " control.inla =
control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert =
control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", "
control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso =
control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg =
inla.arg, num.threads = num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = debug, .parent.frame =
.parent.frame)" )
Time used:
Pre = 1.6, Running = 71.8, Post = 4.72, Total = 78.1
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 0.447 21.437 -44.813 0.561 45.399 0.511 0
UTM_x -0.002 0.156 -0.365 -0.001 0.358 -0.002 0
UTM_y -0.001 0.109 -0.251 -0.001 0.245 -0.001 0
Random effects:
Name Model
spatial.field SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 521.16 4922.7 4.24e-01 45.42 3571.10 0.15
Range for spatial.field 52442.31 6825.159 4.03e+04 52001.18 67125.71 51125.46
Stdev for spatial.field 7.80 0.688 6.55e+00 7.77 9.26 7.69
Deviance Information Criterion (DIC) ...............: 7.83
Deviance Information Criterion (DIC, saturated) ....: 230.95
Effective number of parameters .....................: 121.95
Watanabe-Akaike information criterion (WAIC) ...: 107.98
Effective number of parameters .................: 191.47
Marginal log-Likelihood: 1406.62
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')6.5.1.2 Range: 10km and Standard Deviation 0.5
formula1005 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde1005,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output1005 <- inla(formula1005,
data=inla.stack.data(stack1005, spde=Earthquakes_spde1005),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack1005), compute=TRUE))
summary(output1005)Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles
= quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials,
strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ",
" lincomb = lincomb, selection = selection, control.compute = control.compute, ", "
control.predictor = control.predictor, control.family = control.family, ", " control.inla =
control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert =
control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", "
control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso =
control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg =
inla.arg, num.threads = num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = debug, .parent.frame =
.parent.frame)" )
Time used:
Pre = 1.06, Running = 85.3, Post = 5.57, Total = 91.9
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 0.381 23.228 -48.436 0.509 48.916 0.457 0
UTM_x -0.002 0.188 -0.427 -0.001 0.420 -0.002 0
UTM_y -0.001 0.141 -0.320 -0.001 0.315 -0.001 0
Random effects:
Name Model
spatial.field SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 617.62 8466.65 4.56e-01 33.66 3825.71 3.33e-01
Range for spatial.field 44484.79 8869.88 2.96e+04 43633.15 64330.83 4.20e+04
Stdev for spatial.field 15.18 2.12 1.15e+01 15.02 19.81 1.47e+01
Deviance Information Criterion (DIC) ……………: 24.62
Deviance Information Criterion (DIC, saturated) ….: 198.32
Effective number of parameters …………………: 89.32
Watanabe-Akaike information criterion (WAIC) …: 151.24
Effective number of parameters ……………..: 194.68
Marginal log-Likelihood: 1461.83
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')6.5.1.3 Range: 10km and Standard Deviation 1
formula101 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde101,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output101 <- inla(formula101,
data=inla.stack.data(stack101, spde=Earthquakes_spde101),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack101), compute=TRUE))
summary(output101)Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles
= quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials,
strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ",
" lincomb = lincomb, selection = selection, control.compute = control.compute, ", "
control.predictor = control.predictor, control.family = control.family, ", " control.inla =
control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert =
control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", "
control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso =
control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg =
inla.arg, num.threads = num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = debug, .parent.frame =
.parent.frame)" )
Time used:
Pre = 1.04, Running = 84.1, Post = 7.75, Total = 92.9
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 0.378 23.224 -48.368 0.494 48.834 0.461 0
UTM_x -0.002 0.190 -0.433 -0.001 0.426 -0.002 0
UTM_y -0.001 0.145 -0.329 -0.001 0.325 -0.001 0
Random effects:
Name Model
spatial.field SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 308.69 2696.16 6.23e-01 31.27 2099.73 6.42e-01
Range for spatial.field 41752.89 9993.45 2.55e+04 40626.13 64561.24 3.85e+04
Stdev for spatial.field 20.10 3.44 1.43e+01 19.78 27.75 1.91e+01
Deviance Information Criterion (DIC) ...............: 31.60
Deviance Information Criterion (DIC, saturated) ....: 202.69
Effective number of parameters .....................: 93.69
Watanabe-Akaike information criterion (WAIC) ...: 145.26
Effective number of parameters .................: 182.24
Marginal log-Likelihood: 1473.14
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')6.5.1.4 Range: 10km and Standard Deviation 1.5
formula1015 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde1015,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output1015 <- inla(formula1015,
data=inla.stack.data(stack1015, spde=Earthquakes_spde1015),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack1015), compute=TRUE))
summary(output1015)Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles
= quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials,
strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ",
" lincomb = lincomb, selection = selection, control.compute = control.compute, ", "
control.predictor = control.predictor, control.family = control.family, ", " control.inla =
control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert =
control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", "
control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso =
control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg =
inla.arg, num.threads = num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = debug, .parent.frame =
.parent.frame)" )
Time used:
Pre = 1.03, Running = 86.9, Post = 8.69, Total = 96.6
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 0.376 23.269 -48.488 0.496 48.950 0.459 0
UTM_x -0.002 0.193 -0.439 -0.001 0.432 -0.002 0
UTM_y -0.001 0.148 -0.336 -0.001 0.331 -0.001 0
Random effects:
Name Model
spatial.field SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 356.77 3493.36 5.54e-01 31.23 2390.21 5.11e-01
Range for spatial.field 40305.75 10764.65 2.31e+04 38991.97 65128.75 3.65e+04
Stdev for spatial.field 23.63 4.55 1.61e+01 23.16 33.92 2.22e+01
Deviance Information Criterion (DIC) ...............: 31.91
Deviance Information Criterion (DIC, saturated) ....: 202.41
Effective number of parameters .....................: 93.42
Watanabe-Akaike information criterion (WAIC) ...: 150.38
Effective number of parameters .................: 186.65
Marginal log-Likelihood: 1477.87
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')6.5.1.5 Range: 10km and Standard Deviation 2
formula102 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde102,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output102 <- inla(formula102,
data=inla.stack.data(stack102, spde=Earthquakes_spde102),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack102), compute=TRUE))
summary(output102)Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles
= quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials,
strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ",
" lincomb = lincomb, selection = selection, control.compute = control.compute, ", "
control.predictor = control.predictor, control.family = control.family, ", " control.inla =
control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert =
control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", "
control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso =
control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg =
inla.arg, num.threads = num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = debug, .parent.frame =
.parent.frame)" )
Time used:
Pre = 1.09, Running = 75.4, Post = 9.86, Total = 86.4
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 0.374 23.343 -48.616 0.491 49.071 0.458 0
UTM_x -0.002 0.194 -0.441 -0.001 0.434 -0.002 0
UTM_y -0.001 0.149 -0.338 -0.001 0.333 -0.001 0
Random effects:
Name Model
spatial.field SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 313.11 2821.83 6.15e-01 30.52 2118.86 6.49e-01
Range for spatial.field 39281.81 11361.33 2.14e+04 37806.95 65704.47 3.51e+04
Stdev for spatial.field 26.51 5.57 1.74e+01 25.89 39.23 2.46e+01
Deviance Information Criterion (DIC) ...............: 33.01
Deviance Information Criterion (DIC, saturated) ....: 200.59
Effective number of parameters .....................: 91.58
Watanabe-Akaike information criterion (WAIC) ...: 150.07
Effective number of parameters .................: 183.50
Marginal log-Likelihood: 1480.51
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')6.5.2 20km x 20km, 30km x 30km, 40km x 40km, and 50km x 50km
Do the same steps like 10km x 10km for 20km, 30km, 40km, and 50km with following syntax
#=========================== 7.2 20km * 20km
Grid.CoordUTM2<-as.data.frame(coordsGrid.UTM2)
m2<-nrow(Grid.CoordUTM2)
Grid.CoordUTM2<-as.matrix(Grid.CoordUTM2)
GroupTime2<-rep(c(1:12),each=m2)
All.Grid2<-kronecker(matrix(1, 12, 1), Grid.CoordUTM2)
#============ 7.2.1 PENALIZED COMPLEXITY : 20km, 0.1
prior.median.sd2001 = 0.1; prior.median.range20 = 20000
Earthquakes_spde2001 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range20, .5),
prior.sigma = c(prior.median.sd2001, .5))
s_index2001 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde2001$n.spde,
n.group=12)
stack_est2001 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index2001) , tag="est")
A_pred2 <- inla.spde.make.A(mesh=Earthquakes_mesh,
loc=as.matrix(All.Grid2),
group=GroupTime2, #selected day for prediction
n.group=12)
Cov2<-data.frame(UTM_x=All.Grid2[,1]/1000,UTM_y=All.Grid2[,2]/1000,DEPTH=NA,Time=GroupTime2)
stack_pred2001 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred2),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid2)),
Time=GroupTime2, UTM_x=Cov2$UTM_x,
UTM_y=Cov2$UTM_y, DEPTH=Cov2$DEPTH), s_index2001),
tag="pred")
stack2001 <- inla.stack(stack_est2001, stack_pred2001)
formula2001 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde2001,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output2001 <- inla(formula2001,
data=inla.stack.data(stack2001, spde=Earthquakes_spde2001),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack2001), compute=TRUE))
summary(output2001)
#============ 7.2.2 PENALIZED COMPLEXITY : 20km, 0.5
prior.median.sd2005 = 0.5; prior.median.range20 = 20000
Earthquakes_spde2005 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range20, .5),
prior.sigma = c(prior.median.sd2005, .5))
s_index2005 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde2005$n.spde,
n.group=12)
stack_est2005 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index2005) , tag="est")
stack_pred2005 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred2),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid2)),
Time=GroupTime2, UTM_x=Cov2$UTM_x,
UTM_y=Cov2$UTM_y, DEPTH=Cov2$DEPTH), s_index2005),
tag="pred")
stack2005 <- inla.stack(stack_est2005, stack_pred2005)
formula2005 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde2005,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output2005 <- inla(formula2005,
data=inla.stack.data(stack2005, spde=Earthquakes_spde2005),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack2005), compute=TRUE))
summary(output2005)
#============ 7.2.3 PENALIZED COMPLEXITY : 20km, 1
prior.median.sd201 = 1; prior.median.range20 = 20000
Earthquakes_spde201 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range20, .5),
prior.sigma = c(prior.median.sd201, .5))
s_index201 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde201$n.spde,
n.group=12)
stack_est201 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index201) , tag="est")
stack_pred201 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred2),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid2)),
Time=GroupTime2, UTM_x=Cov2$UTM_x,
UTM_y=Cov2$UTM_y, DEPTH=Cov2$DEPTH), s_index201),
tag="pred")
stack201 <- inla.stack(stack_est201, stack_pred201)
formula201 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde201,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output201 <- inla(formula201,
data=inla.stack.data(stack201, spde=Earthquakes_spde201),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack201), compute=TRUE))
summary(output201)
#============ 7.2.4 PENALIZED COMPLEXITY : 20km, 1.5
prior.median.sd2015 = 1.5; prior.median.range20 = 20000
Earthquakes_spde2015 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range20, .5),
prior.sigma = c(prior.median.sd2015, .5))
s_index2015 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde2015$n.spde,
n.group=12)
stack_est2015 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index2015) , tag="est")
stack_pred2015 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred2),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid2)),
Time=GroupTime2, UTM_x=Cov2$UTM_x,
UTM_y=Cov2$UTM_y, DEPTH=Cov2$DEPTH), s_index2015),
tag="pred")
stack2015 <- inla.stack(stack_est2015, stack_pred2015)
formula2015 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde2015,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output2015 <- inla(formula2015,
data=inla.stack.data(stack2015, spde=Earthquakes_spde2015),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack2015), compute=TRUE))
summary(output2015)
#============ 7.2.5 PENALIZED COMPLEXITY : 20km, 2
prior.median.sd202 = 2; prior.median.range20 = 20000
Earthquakes_spde202 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range20, .5),
prior.sigma = c(prior.median.sd202, .5))
s_index202 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde202$n.spde,
n.group=12)
stack_est202 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index202) , tag="est")
stack_pred202 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred2),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid2)),
Time=GroupTime2, UTM_x=Cov2$UTM_x,
UTM_y=Cov2$UTM_y, DEPTH=Cov2$DEPTH), s_index202),
tag="pred")
stack202 <- inla.stack(stack_est202, stack_pred202)
formula202 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde202,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output202 <- inla(formula202,
data=inla.stack.data(stack202, spde=Earthquakes_spde202),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack202), compute=TRUE))
summary(output202)
#=========================== 7.3 30km * 30km
Grid.CoordUTM3<-as.data.frame(coordsGrid.UTM3)
m3<-nrow(Grid.CoordUTM3)
Grid.CoordUTM3<-as.matrix(Grid.CoordUTM3)
GroupTime3<-rep(c(1:12),each=m3)
All.Grid3<-kronecker(matrix(1, 12, 1), Grid.CoordUTM3)
#============ 7.3.1 PENALIZED COMPLEXITY : 30km, 0.1
prior.median.sd3001 = 0.1; prior.median.range30 = 30000
Earthquakes_spde3001 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range30, .5),
prior.sigma = c(prior.median.sd3001, .5))
s_index3001 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde3001$n.spde,
n.group=12)
stack_est3001 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index3001) , tag="est")
A_pred3 <- inla.spde.make.A(mesh=Earthquakes_mesh,
loc=as.matrix(All.Grid3),
group=GroupTime3, #selected day for prediction
n.group=12)
Cov3<-data.frame(UTM_x=All.Grid3[,1]/1000,UTM_y=All.Grid3[,2]/1000,DEPTH=NA,Time=GroupTime3)
stack_pred3001 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred3),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid3)),
Time=GroupTime3, UTM_x=Cov3$UTM_x,
UTM_y=Cov3$UTM_y, DEPTH=Cov3$DEPTH), s_index3001),
tag="pred")
stack3001 <- inla.stack(stack_est3001, stack_pred3001)
formula3001 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde3001,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output3001 <- inla(formula3001,
data=inla.stack.data(stack3001, spde=Earthquakes_spde3001),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack3001), compute=TRUE))
summary(output3001)
#============ 7.3.2 PENALIZED COMPLEXITY : 30km, 0.5
prior.median.sd3005 = 0.5; prior.median.range30 = 30000
Earthquakes_spde3005 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range30, .5),
prior.sigma = c(prior.median.sd3005, .5))
s_index3005 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde3005$n.spde,
n.group=12)
stack_est3005 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index3005) , tag="est")
stack_pred3005 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred3),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid3)),
Time=GroupTime3, UTM_x=Cov3$UTM_x,
UTM_y=Cov3$UTM_y, DEPTH=Cov3$DEPTH), s_index3005),
tag="pred")
stack3005 <- inla.stack(stack_est3005, stack_pred3005)
formula3005 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde3005,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output3005 <- inla(formula3005,
data=inla.stack.data(stack3005, spde=Earthquakes_spde3005),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack3005), compute=TRUE))
summary(output3005)
#============ 7.3.3 PENALIZED COMPLEXITY : 30km, 1
prior.median.sd301 = 1; prior.median.range30 = 30000
Earthquakes_spde301 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range30, .5),
prior.sigma = c(prior.median.sd301, .5))
s_index301 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde301$n.spde,
n.group=12)
stack_est301 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index301) , tag="est")
stack_pred301 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred3),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid3)),
Time=GroupTime3, UTM_x=Cov3$UTM_x,
UTM_y=Cov3$UTM_y, DEPTH=Cov3$DEPTH), s_index301),
tag="pred")
stack301 <- inla.stack(stack_est301, stack_pred301)
formula301 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde301,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output301 <- inla(formula301,
data=inla.stack.data(stack301, spde=Earthquakes_spde301),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack301), compute=TRUE))
summary(output301)
#============ 7.3.4 PENALIZED COMPLEXITY : 30km, 1.5
prior.median.sd3015 = 1.5; prior.median.range30 = 30000
Earthquakes_spde3015 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range30, .5),
prior.sigma = c(prior.median.sd3015, .5))
s_index3015 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde3015$n.spde,
n.group=12)
stack_est3015 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index3015) , tag="est")
stack_pred3015 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred3),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid3)),
Time=GroupTime3, UTM_x=Cov3$UTM_x,
UTM_y=Cov3$UTM_y, DEPTH=Cov3$DEPTH), s_index3015),
tag="pred")
stack3015 <- inla.stack(stack_est3015, stack_pred3015)
formula3015 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde3015,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output3015 <- inla(formula3015,
data=inla.stack.data(stack3015, spde=Earthquakes_spde3015),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack3015), compute=TRUE))
summary(output3015)
#============ 7.3.5 PENALIZED COMPLEXITY : 30km, 2
prior.median.sd302 = 2; prior.median.range30 = 30000
Earthquakes_spde302 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range30, .5),
prior.sigma = c(prior.median.sd302, .5))
s_index302 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde302$n.spde,
n.group=12)
stack_est302 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index302) , tag="est")
stack_pred302 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred3),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid3)),
Time=GroupTime3, UTM_x=Cov3$UTM_x,
UTM_y=Cov3$UTM_y, DEPTH=Cov3$DEPTH), s_index302),
tag="pred")
stack302 <- inla.stack(stack_est302, stack_pred302)
formula302 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde302,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output302 <- inla(formula302,
data=inla.stack.data(stack302, spde=Earthquakes_spde302),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack302), compute=TRUE))
summary(output302)
#=========================== 7.4 40km * 40km
Grid.CoordUTM4<-as.data.frame(coordsGrid.UTM4)
m4<-nrow(Grid.CoordUTM4)
Grid.CoordUTM4<-as.matrix(Grid.CoordUTM4)
GroupTime4<-rep(c(1:12),each=m4)
All.Grid4<-kronecker(matrix(1, 12, 1), Grid.CoordUTM4)
#============ 7.4.1 PENALIZED COMPLEXITY : 40km, 0.1
prior.median.sd4001 = 0.1; prior.median.range40 = 40000
Earthquakes_spde4001 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range40, .5),
prior.sigma = c(prior.median.sd4001, .5))
s_index4001 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde4001$n.spde,
n.group=12)
stack_est4001 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index4001) , tag="est")
A_pred4 <- inla.spde.make.A(mesh=Earthquakes_mesh,
loc=as.matrix(All.Grid4),
group=GroupTime4, #selected day for prediction
n.group=12)
Cov4<-data.frame(UTM_x=All.Grid4[,1]/1000,UTM_y=All.Grid4[,2]/1000,DEPTH=NA,Time=GroupTime4)
stack_pred4001 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred4),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid4)),
Time=GroupTime4, UTM_x=Cov4$UTM_x,
UTM_y=Cov4$UTM_y, DEPTH=Cov4$DEPTH), s_index4001),
tag="pred")
stack4001 <- inla.stack(stack_est4001, stack_pred4001)
formula4001 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde4001,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output4001 <- inla(formula4001,
data=inla.stack.data(stack4001, spde=Earthquakes_spde4001),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack4001), compute=TRUE))
summary(output4001)
#============ 7.4.2 PENALIZED COMPLEXITY : 40km, 0.5
prior.median.sd4005 = 0.5; prior.median.range40 = 40000
Earthquakes_spde4005 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range40, .5),
prior.sigma = c(prior.median.sd4005, .5))
s_index4005 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde4005$n.spde,
n.group=12)
stack_est4005 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index4005) , tag="est")
stack_pred4005 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred4),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid4)),
Time=GroupTime4, UTM_x=Cov4$UTM_x,
UTM_y=Cov4$UTM_y, DEPTH=Cov4$DEPTH), s_index4005),
tag="pred")
stack4005 <- inla.stack(stack_est4005, stack_pred4005)
formula4005 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde4005,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output4005 <- inla(formula4005,
data=inla.stack.data(stack4005, spde=Earthquakes_spde4005),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack4005), compute=TRUE))
summary(output4005)
#============ 7.4.3 PENALIZED COMPLEXITY : 40km, 1
prior.median.sd401 = 1; prior.median.range40 = 40000
Earthquakes_spde401 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range40, .5),
prior.sigma = c(prior.median.sd401, .5))
s_index401 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde401$n.spde,
n.group=12)
stack_est401 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index401) , tag="est")
stack_pred401 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred4),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid4)),
Time=GroupTime4, UTM_x=Cov4$UTM_x,
UTM_y=Cov4$UTM_y, DEPTH=Cov4$DEPTH), s_index401),
tag="pred")
stack401 <- inla.stack(stack_est401, stack_pred401)
formula401 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde401,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output401 <- inla(formula401,
data=inla.stack.data(stack401, spde=Earthquakes_spde401),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack401), compute=TRUE))
summary(output401)
#============ 7.4.4 PENALIZED COMPLEXITY : 40km, 1.5
prior.median.sd4015 = 1.5; prior.median.range40 = 40000
Earthquakes_spde4015 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range40, .5),
prior.sigma = c(prior.median.sd4015, .5))
s_index4015 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde4015$n.spde,
n.group=12)
stack_est4015 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index4015) , tag="est")
stack_pred4015 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred4),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid4)),
Time=GroupTime4, UTM_x=Cov4$UTM_x,
UTM_y=Cov4$UTM_y, DEPTH=Cov4$DEPTH), s_index4015),
tag="pred")
stack4015 <- inla.stack(stack_est4015, stack_pred4015)
formula4015 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde4015,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output4015 <- inla(formula4015,
data=inla.stack.data(stack4015, spde=Earthquakes_spde4015),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack4015), compute=TRUE))
summary(output4015)
#============ 7.4.5 PENALIZED COMPLEXITY : 40km, 2
prior.median.sd402 = 2; prior.median.range40 = 40000
Earthquakes_spde402 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range40, .5),
prior.sigma = c(prior.median.sd402, .5))
s_index402 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde402$n.spde,
n.group=12)
stack_est402 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index402) , tag="est")
stack_pred402 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred4),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid4)),
Time=GroupTime4, UTM_x=Cov4$UTM_x,
UTM_y=Cov4$UTM_y, DEPTH=Cov4$DEPTH), s_index402),
tag="pred")
stack402 <- inla.stack(stack_est402, stack_pred402)
formula402 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde402,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output402 <- inla(formula402,
data=inla.stack.data(stack402, spde=Earthquakes_spde402),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack402), compute=TRUE))
summary(output402)
#=========================== 7.5 50km * 50km
Grid.CoordUTM5<-as.data.frame(coordsGrid.UTM5)
m5<-nrow(Grid.CoordUTM5)
Grid.CoordUTM5<-as.matrix(Grid.CoordUTM5)
GroupTime5<-rep(c(1:12),each=m5)
All.Grid5<-kronecker(matrix(1, 12, 1), Grid.CoordUTM5)
#============ 7.5.1 PENALIZED COMPLEXITY : 50km, 0.1
prior.median.sd5001 = 0.1; prior.median.range50 = 50000
Earthquakes_spde5001 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range50, .5),
prior.sigma = c(prior.median.sd5001, .5))
s_index5001 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde5001$n.spde,
n.group=12)
stack_est5001 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index5001) , tag="est")
A_pred5 <- inla.spde.make.A(mesh=Earthquakes_mesh,
loc=as.matrix(All.Grid5),
group=GroupTime5, #selected day for prediction
n.group=12)
Cov5<-data.frame(UTM_x=All.Grid5[,1]/1000,UTM_y=All.Grid5[,2]/1000,DEPTH=NA,Time=GroupTime5)
stack_pred5001 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred5),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid5)),
Time=GroupTime5, UTM_x=Cov5$UTM_x,
UTM_y=Cov5$UTM_y, DEPTH=Cov5$DEPTH), s_index5001),
tag="pred")
stack5001 <- inla.stack(stack_est5001, stack_pred5001)
formula5001 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde5001,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output5001 <- inla(formula5001,
data=inla.stack.data(stack5001, spde=Earthquakes_spde5001),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack5001), compute=TRUE))
summary(output5001)
#============ 7.5.2 PENALIZED COMPLEXITY : 50km, 0.5
prior.median.sd5005 = 0.5; prior.median.range50 = 50000
Earthquakes_spde5005 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range50, .5),
prior.sigma = c(prior.median.sd5005, .5))
s_index5005 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde5005$n.spde,
n.group=12)
stack_est5005 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index5005) , tag="est")
stack_pred5005 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred5),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid5)),
Time=GroupTime5, UTM_x=Cov5$UTM_x,
UTM_y=Cov5$UTM_y, DEPTH=Cov5$DEPTH), s_index5005),
tag="pred")
stack5005 <- inla.stack(stack_est5005, stack_pred5005)
formula5005 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde5005,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output5005 <- inla(formula5005,
data=inla.stack.data(stack5005, spde=Earthquakes_spde5005),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack5005), compute=TRUE))
summary(output5005)
#============ 7.5.3 PENALIZED COMPLEXITY : 50km, 1
prior.median.sd501 = 1; prior.median.range50 = 50000
Earthquakes_spde501 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range50, .5),
prior.sigma = c(prior.median.sd501, .5))
s_index501 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde501$n.spde,
n.group=12)
stack_est501 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index501) , tag="est")
stack_pred501 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred5),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid5)),
Time=GroupTime5, UTM_x=Cov5$UTM_x,
UTM_y=Cov5$UTM_y, DEPTH=Cov5$DEPTH), s_index501),
tag="pred")
stack501 <- inla.stack(stack_est501, stack_pred501)
formula501 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde501,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output501 <- inla(formula501,
data=inla.stack.data(stack501, spde=Earthquakes_spde501),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack501), compute=TRUE))
summary(output501)
#============ 7.5.4 PENALIZED COMPLEXITY : 50km, 1.5
prior.median.sd5015 = 1.5; prior.median.range50 = 50000
Earthquakes_spde5015 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range50, .5),
prior.sigma = c(prior.median.sd5015, .5))
s_index5015 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde5015$n.spde,
n.group=12)
stack_est5015 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index5015) , tag="est")
stack_pred5015 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred5),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid5)),
Time=GroupTime5, UTM_x=Cov5$UTM_x,
UTM_y=Cov5$UTM_y, DEPTH=Cov5$DEPTH), s_index5015),
tag="pred")
stack5015 <- inla.stack(stack_est5015, stack_pred5015)
formula5015 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde5015,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output5015 <- inla(formula5015,
data=inla.stack.data(stack5015, spde=Earthquakes_spde5015),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack5015), compute=TRUE))
summary(output5015)
#============ 7.5.5 PENALIZED COMPLEXITY : 50km, 2
prior.median.sd502 = 2; prior.median.range50 = 50000
Earthquakes_spde502 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range50, .5),
prior.sigma = c(prior.median.sd502, .5))
s_index502 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde502$n.spde,
n.group=12)
stack_est502 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index502) , tag="est")
stack_pred502 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred5),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid5)),
Time=GroupTime5, UTM_x=Cov5$UTM_x,
UTM_y=Cov5$UTM_y, DEPTH=Cov5$DEPTH), s_index502),
tag="pred")
stack502 <- inla.stack(stack_est502, stack_pred502)
formula502 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde502,group=spatial.field.group,
control.group=list(model="rw1"))
# ATTENTION: the run is computationally intensive!
output502 <- inla(formula502,
data=inla.stack.data(stack502, spde=Earthquakes_spde502),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack502), compute=TRUE))
summary(output502)7. Model Selection
7.1 Table DIC, WAIC, logCPO
model_outputs <- list(
output1001, output1005, output101, output1015, output102,
output2001, output2005, output201, output2015, output202,
output3001, output3005, output301, output3015, output302,
output4001, output4005, output401, output4015, output402,
output5001, output5005, output501, output5015, output502
)
model_names <- c(
"output1001", "output1005", "output101", "output1015", "output102",
"output2001", "output2005", "output201", "output2015", "output202",
"output3001", "output3005", "output301", "output3015", "output302",
"output4001", "output4005", "output401", "output4015", "output402",
"output5001", "output5005", "output501", "output5015", "output502"
)
model_metrics <- lapply(model_outputs, function(output) {
list(
DIC = output$dic$dic,
WAIC = output$waic$waic,
LogCPO = sum(log(output$cpo$cpo), na.rm = TRUE)
)
})
model_selection <- data.frame(
model = model_names,
DIC = round(sapply(model_metrics, function(x) x$DIC), 2),
WAIC = round(sapply(model_metrics, function(x) x$WAIC), 2),
LogCPO = round(sapply(model_metrics, function(x) x$LogCPO), 2)
)
print(model_selection)7.2 Plot DIC
range_km <- c(10,20,30,40,50,
10,20,30,40,50,
10,20,30,40,50,
10,20,30,40,50,
10,20,30,40,50)
std <- c(0.1, 0.5, 1, 1.5, 2,
0.1, 0.5, 1, 1.5, 2,
0.1, 0.5, 1, 1.5, 2,
0.1, 0.5, 1, 1.5, 2,
0.1, 0.5, 1, 1.5, 2)
DIC <- c(7.84, 24.47, 30.93, 32.30, 33.05,
8.30, 27.40, 30.44, 31.97, 32.76,
7.93, 26.42, 30.98, 32.36, 32.63,
3.12, 30.20, 30.12, 31.94, 32.55,
6.21, 29.65, 30.76, 31.71, 32.71)
z <- matrix(DIC, nrow = length(unique(std)), byrow = TRUE)
# Create the surface plot
plot_ly(x = ~std, y = ~range_km, z = ~z, type = "surface",
colorscale = list(list(0, "cornsilk"), list(0.5, "bisque"), list(1, "darkorange")),
colorbar = list(title = "DIC", orientation = "h", x = 0.5, y = -0.2)) %>%
layout(title = "Surface Plot of DIC",
scene = list(xaxis = list(title = "Standard Deviation"),
yaxis = list(title = "Range (km)"),
zaxis = list(title = "DIC")))
7.3 Plot WAIC
# Define the values based on your image
WAIC <- c(107.93, 152.24, 146.73, 150.01, 150.21,
114.73, 139.88, 147.94, 150.78, 150.72,
108.88, 142.59, 146.52, 148.25, 150.99,
113.84, 134.07, 148.45, 149.69, 151.03,
109.31, 133.67, 145.94, 149.60, 150.55)
z1 <- matrix(WAIC, nrow = length(unique(std)), byrow = TRUE)
# Create the surface plot
plot_ly(x = ~std, y = ~range_km, z = ~z1, type = "surface",
colorscale = list(list(0, "cornsilk"), list(0.5, "bisque"), list(1, "darkorange4")),
colorbar = list(title = "WAIC", orientation = "h", x = 0.5, y = -0.2)) %>%
layout(title = "Surface Plot of WAIC",
scene = list(xaxis = list(title = "Standard Deviation"),
yaxis = list(title = "Range (km)"),
zaxis = list(title = "WAIC")))
7.4 Plot logCPO
# Define the values based on your image
CPO <- c(-155.68, -177.73, -174.83, -176.24, -176.36,
-160.16, -171.65, -175.36, -176.60, -176.59,
-155.96, -172.91, -174.69, -175.50, -176.70,
-159.37, -168.27, -175.58, -176.12, -176.73,
-156.85, -167.09, -174.47, -176.09, -176.46)
z2 <- matrix(CPO, nrow = length(unique(std)), byrow = TRUE)
quartz()
# Create the surface plot
plot_ly(x = ~std, y = ~range_km, z = ~z2, type = "surface",
colorscale = list(list(0, "cornsilk"), list(0.5, "bisque"), list(1, "darkorange4")),
colorbar = list(title = "logCPO", orientation = "h", x = 0.5, y = -0.2)) %>%
layout(title = "Surface Plot of logCPO",
scene = list(xaxis = list(title = "Standard Deviation"),
yaxis = list(title = "Range (km)"),
zaxis = list(title = "logCPO")))
Chosen Model
#summary(output1001)
plot(output1001)
8. Calculate the exceedance probability to serve as the foundation for spatiotemporal clustering
To calculate the exceedance probability, we need to establish the threshold for what constitutes an earthquake of interest, such as those with a magnitude of 5 Mw or higher.
MAG<-log(5)
ExProb <- unlist(lapply(output1001$marginals.fitted.values, function(X) {
if (all(is.finite(X))) {
return(1 - inla.pmarginal(MAG, X))
} else {
return(NA) # or any other appropriate handling
}
}))9. Visualization
9.1 Prepare output visualization in map
index_pred <- inla.stack.index(stack1001,"pred")$data
index_est <- inla.stack.index(stack1001,"est")$data
Earthquakes_Pred<-data.frame(x=All.Grid1[,1],y=All.Grid1[,2],Time=GroupTime1)
Earthquakes_Est<-data.frame(x=Earthquakes20$UTM_x,y=Earthquakes20$UTM_y, Time=Earthquakes20$Month)
Earthquakes_Pred$pred_mean <- exp(output1001$summary.fitted.values[index_pred, "mean"])
Earthquakes_Pred$pred_ll <- exp(output1001$summary.fitted.values[index_pred, "0.025quant"])
Earthquakes_Pred$pred_ul <- exp(output1001$summary.fitted.values[index_pred, "0.975quant"])Earthquakes_Est$pred_mean <- exp(output1001$summary.fitted.values[index_est, "mean"])
Earthquakes_Est$pred_ll <- exp(output1001$summary.fitted.values[index_est, "0.025quant"])
Earthquakes_Est$pred_ul <- exp(output1001$summary.fitted.values[index_est, "0.975quant"])
summary(Earthquakes_Est$pred_mean)Earthquakes_Values<-rbind(Earthquakes_Est,Earthquakes_Pred)
Earthquakes_Pred$pred_prob <- ExProb[index_pred]
write.csv(Earthquakes_Pred,"Earthquakes_Pred.csv")
write.csv(Earthquakes_Est,"Earthquakes_Est.csv")
Earthquakes_Pred<-read.csv("Earthquakes_Pred.csv", sep=",")Summary estimation result
summary(Earthquakes_Est$pred_mean)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.998 4.400 4.700 4.775 5.000 6.200Summary prediction result
summary(Earthquakes_Est$pred_mean)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.998 4.400 4.700 4.775 5.000 6.2009.2 High-Resolution of Earthquake Interpolation
9.2.1 Overall: By Month
Earthquakes_Pred$title <- "By Month"
all = ggplot() +
geom_tile(data = Earthquakes_Pred, aes(x = x, y = y, fill = pred_mean)) +
geom_contour(data = Earthquakes_Pred, aes(x = x, y = y, z = pred_mean, colour = stat(level)), size = 0.2, colour = "red", alpha = 0.5) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red")) +
facet_wrap(~ Time, labeller = as_labeller(Month_names), ncol = 4) +
theme_bw() +
ylab("") +
xlab("") +
geom_sf(data = Sumatra_sf_utm1, fill = NA, size=0.01) +
theme(legend.position = "bottom", text = element_text(size = 19)) +
guides(fill = guide_colorbar(title.position = "left", title.vjust = 1,
frame.colour = "black",
barwidth = 20,
barheight = 1.5)) +
theme(panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
axis.text.x = element_blank(), #remove x axis labels
axis.ticks.x = element_blank(), #remove x axis ticks
axis.text.y = element_blank(), #remove y axis labels
axis.ticks.y = element_blank() #remove y axis ticks
) +
labs(fill = "Earthquakes \nMagnitude")
all
9.2.2 Overall: By Province
For this visualization we choose one periode, because there are no significant different secara temporal.
Earthquakes_Pred_dec <- Earthquakes_Pred[Earthquakes_Pred$Time == 12, ]
Earthquakes_Pred_dec$title = "By Province"
ggplot() +
geom_tile(data = Earthquakes_Pred_dec, aes(x = x, y = y, fill = pred_mean)) +
geom_contour(data = Earthquakes_Pred_dec, aes(x = x, y = y, z = pred_mean, colour = stat(level)),
size = 0.2, colour = "red", alpha = 0.5) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red")) +
theme_bw() +
ylab("Latitude") +
xlab("Longitude") +
geom_sf(data = Sumatra_sf_utm1, fill = NA, size = 0.01) +
geom_sf_text(data = Sumatra_sf_utm1, aes(label = NAME_1), size = 3, color = "black", fontface = "bold") +
geom_sf(data = patahan_sumatera, color = "blue", size = 1) +
theme(legend.position = "bottom",
text = element_text(size = 12)) +
theme(
panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
plot.title = element_text(hjust = 0.5, size = 12), # Menyesuaikan ukuran dan posisi judul
panel.grid = element_blank() # Menghilangkan grid
)+
guides(fill = guide_colorbar(title.position = "left", title.vjust = 1,
frame.colour = "black",
barwidth = 20,
barheight = 1.5))+
facet_grid(. ~ title)
9.2.3 District Level
max_pred = max(Earthquakes_Pred$pred_mean)
kab = st_read(file.choose())
kab <- st_transform(kab, crs = st_crs(Sumatra))
kab_valid <- st_make_valid(kab)9.2.3.1 Bengkulu Province
Sumatra_Bengkulu <- Sumatra_sf_utm1 %>% filter(NAME_1 == 'Bengkulu')
Sumatra_Bengkulu_wgs84 <- st_transform(Sumatra_Bengkulu, crs = st_crs(kab))
bbox_Sumatra_Bengkulu <- st_bbox(Sumatra_Bengkulu_wgs84)
# Crop kab based on bounding box from Sumatra_Bengkulu
kab_cropped <- st_crop(kab_valid, bbox_Sumatra_Bengkulu)
kab_bengkulu<-kab_cropped[kab_cropped$KAB_KOTA %in%c("MUKOMUKO","LEBONG","BENGKULU UTARA", "REJANG LEBONG", "BENGKULU TENGAH", "KEPAHIANG","KOTA BENGKULU", "SELUMA", "BENGKULU SELATAN", "KAUR"),]
kab_bengkulu_sf <- st_as_sf(kab_bengkulu)
kab_bengkulu_sf_utm1 <- st_transform(kab_bengkulu_sf, crs = 32648)
Earthquakes_Pred_Bengkulu <- Earthquakes_Pred_dec %>%
filter(x >= min(st_bbox(kab_bengkulu_sf_utm1)$xmin) & x <= max(st_bbox(kab_bengkulu_sf_utm1)$xmax) &
y >= min(st_bbox(kab_bengkulu_sf_utm1)$ymin) & y <= max(st_bbox(kab_bengkulu_sf_utm1)$ymax))
min(Earthquakes_Pred_Bengkulu$pred_mean)
max(Earthquakes_Pred_Bengkulu$pred_mean)
#==== plot with smoothness
Earthquakes_Pred_Bengkulu$x <- jitter(Earthquakes_Pred_Bengkulu$x, amount = 0.0001)
Earthquakes_Pred_Bengkulu$y <- jitter(Earthquakes_Pred_Bengkulu$y, amount = 0.0001)
# Interpolation after jitter
interpolated_bengkulu <- with(Earthquakes_Pred_Bengkulu,
interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200))
# Dataframe format
interpolated_df_beng <- as.data.frame(expand.grid(x = interpolated_bengkulu$x, y = interpolated_bengkulu$y))
interpolated_df_beng$pred_mean <- as.vector(interpolated_bengkulu$z)
interpolated_df_beng <- na.omit(interpolated_df_beng) # Menghapus nilai NA
#Plot
interpolated_df_beng$title = "Bengkulu"
bengkulu = ggplot() +
geom_raster(data = interpolated_df_beng, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) +
geom_sf(data = kab_bengkulu_sf_utm1, fill = NA, color = "black", linewidth = 0.3) +
geom_sf_text(data = kab_bengkulu, aes(label = KAB_KOTA), size = 1.5, color = "black", fontface = "bold") +
coord_sf(crs = 32648) +
theme_bw() +
theme(
panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
) +
labs(
x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude"
) +
facet_grid(. ~ title)
bengkulu
9.2.3.2 Sumatra Barat Province
Sumatra_Barat <- Sumatra_sf_utm1 %>% filter(NAME_1 == 'Sumatera Barat')
Sumatra_Barat_wgs84 <- st_transform(Sumatra_Barat, crs = st_crs(kab))
bbox_Sumatra_Barat <- st_bbox(Sumatra_Barat_wgs84)
kab_cropped <- st_crop(kab_valid, bbox_Sumatra_Barat)
kab_sumbar<-kab_cropped[kab_cropped$KAB_KOTA %in%c("AGAM","DHARMASRAYA", "KEPULAUAN MENTAWAI", "LIMA PULUH KOTA","PADANG PARIAMAN", "PASAMAN", "PASAMAN BARAT", "PESISIR SELATAN","SIJUNJUNG", "SOLOK", "SOLOK SELATAN", "TANAH DATAR", "KOTA BUKITTINGGI", "KOTA PADANG", "KOTA PADANG PANJANG", "KOTA PARIAMAN","KOTA PAYAKUMBUH", "KOTA SAWAHLUNTO", "KOTA SOLOK"),]
kab_sumbar_sf <- st_as_sf(kab_sumbar)
kab_sumbar_sf_utm1 <- st_transform(kab_sumbar_sf, crs = 32648)
Earthquakes_Pred_Sumbar <- Earthquakes_Pred_dec %>%
filter(x >= min(st_bbox(kab_sumbar_sf_utm1)$xmin) & x <= max(st_bbox(kab_sumbar_sf_utm1)$xmax) &
y >= min(st_bbox(kab_sumbar_sf_utm1)$ymin) & y <= max(st_bbox(kab_sumbar_sf_utm1)$ymax))
min(Earthquakes_Pred_Sumbar$pred_mean)
max(Earthquakes_Pred_Sumbar$pred_mean)
Earthquakes_Pred_Sumbar$x <- jitter(Earthquakes_Pred_Sumbar$x, amount = 0.0001)
Earthquakes_Pred_Sumbar$y <- jitter(Earthquakes_Pred_Sumbar$y, amount = 0.0001)
interpolated_sumbar <- with(Earthquakes_Pred_Sumbar,
interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200))
interpolated_df_sumbar <- as.data.frame(expand.grid(x = interpolated_sumbar$x, y = interpolated_sumbar$y))
interpolated_df_sumbar$pred_mean <- as.vector(interpolated_sumbar$z)
interpolated_df_sumbar <- na.omit(interpolated_df_sumbar) # Menghapus nilai NA
#Plot
interpolated_df_sumbar$title = "Sumatra Barat"
sumbar = ggplot() +
geom_raster(data = interpolated_df_sumbar, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) +
geom_sf(data = kab_sumbar_sf_utm1, fill = NA, color = "black", linewidth = 0.3) +
geom_sf_text(data = kab_sumbar, aes(label = KAB_KOTA), size = 1.5, color = "black", fontface = "bold") +
coord_sf(crs = 32648) +
theme_bw() +
theme(
panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
) +
labs(
x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude"
) +
facet_grid(. ~ title)
sumbar
9.2.3.3 Sumatra Utara Province
#==== SUMATERA UTARA
Sumatra_Utara <- Sumatra_sf_utm1 %>% filter(NAME_1 == 'Sumatera Utara')
Sumatra_Utara_wgs84 <- st_transform(Sumatra_Utara, crs = st_crs(kab))
bbox_Sumatra_Utara <- st_bbox(Sumatra_Utara_wgs84)
kab_cropped <- st_crop(kab_valid, bbox_Sumatra_Utara)
kab_sumut<-kab_cropped[kab_cropped$KAB_KOTA %in%c("ASAHAN", "BATU BARA", "DAIRI", "DELI SERDANG", "HUMBANG HASUNDUTAN","KARO", "LABUHANBATU", "LABUHANBATU SELATAN", "LABUHANBATU UTARA", "LANGKAT","MANDAILING NATAL", "NIAS", "NIAS BARAT","NIAS SELATAN", "NIAS UTARA", "PADANG LAWAS","PADANG LAWAS UTARA","PAKPAK BHARAT", "SAMOSIR","SERDANG BEDAGAI", "SIMALUNGUN","TAPANULI SELATAN","TAPANULI TENGAH","TAPANULI UTARA","TOBA SAMOSIR","KOTA BINJAI","KOTA GUNUNGSITOLI", "KOTA MEDAN", "KOTA PADANG SIDIMPUAN","KOTA PEMATANGSIANTAR","KOTA SIBOLGA","KOTA TANJUNG BALAI",
"KOTA TEBING TINGGI"),]
kab_sumut_sf <- st_as_sf(kab_sumut)
kab_sumut_sf_utm1 <- st_transform(kab_sumut_sf, crs = 32648)
Earthquakes_Pred_Sumut <- Earthquakes_Pred_dec %>%
filter(x >= min(st_bbox(kab_sumut_sf_utm1)$xmin) & x <= max(st_bbox(kab_sumut_sf_utm1)$xmax) &
y >= min(st_bbox(kab_sumut_sf_utm1)$ymin) & y <= max(st_bbox(kab_sumut_sf_utm1)$ymax))
#min(Earthquakes_Pred_Sumut$pred_mean)
#max(Earthquakes_Pred_Sumut$pred_mean)
Earthquakes_Pred_Sumut$x <- jitter(Earthquakes_Pred_Sumut$x, amount = 0.0001)
Earthquakes_Pred_Sumut$y <- jitter(Earthquakes_Pred_Sumut$y, amount = 0.0001)
interpolated_sumut <- with(Earthquakes_Pred_Sumut,
interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200))
interpolated_df_sumut <- as.data.frame(expand.grid(x = interpolated_sumut$x, y = interpolated_sumut$y))
interpolated_df_sumut$pred_mean <- as.vector(interpolated_sumut$z)
interpolated_df_sumut <- na.omit(interpolated_df_sumut) # Menghapus nilai NA
#Plot
interpolated_df_sumut$title = "Sumatra Utara"
sumut = ggplot() +
geom_raster(data = interpolated_df_sumut, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) +
geom_sf(data = kab_sumut_sf_utm1, fill = NA, color = "black", linewidth = 0.3) +
geom_sf_text(data = kab_sumut, aes(label = KAB_KOTA), size = 1.5, color = "black", fontface = "bold") +
coord_sf(crs = 32648) +
theme_bw() +
theme(
panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
) +
labs(
x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude"
) +
facet_grid(. ~ title)
sumut
9.2.4 Sub-district Level
kec = st_read(file.choose())
kec <- st_transform(kec, crs = st_crs(Sumatra))
kec_valid <- st_make_valid(kec)9.2.4.1 Bengkulu Province
Mukomuko, bengkulu utara, and lebong
kab_bengkulu_wgs84 <- st_transform(kab_bengkulu, crs = st_crs(kec)) bbox_kab_bengkulu <- st_bbox(kab_bengkulu_wgs84) kec_cropped1 <- st_crop(kec_valid, bbox_kab_bengkulu) kec_mbul<-kec_cropped1[kec_cropped1$KECAMATAN %in%c("V KOTO", "XIV KOTO","AIR DIKIT","AIR MAJUNTO","AIR RAMI", "IPUH","KOTA MUKOMUKO","LUBUK PINANG","MALIN DEMAN","PENARIK","PONDOK SUGUH","SELAGAN RAYA","SUNGAI RUMBAI","TERAMANG JAYA", "TERAS TERUNJAM", "AIR BESI","AIR NAPAL", "AIR PADANG", "ARMA JAYA","BATIK NAU", "GIRI MULYA","HULU PALIK","KERKAP","KETAHUN", "KOTA ARGA MAKMUR","LAIS","MARGA SAKTI SEBELAT","NAPAL PUTIH","PADANG JAYA","PINANG RAYA","PUTRI HIJAU","TANJUNG AGUNG PALIK","ULOK KUPAI", "AMEN","BINGIN KUNING", "LEBONG ATAS", "LEBONG SAKTI","LEBONG SELATAN","LEBONG TENGAH","LEBONG UTARA","PINANG BELAPIS","RIMBO PENGADANG","TOPOS","TUBEI","URAM JAYA"),] kec_mbul_sf <- st_as_sf(kec_mbul) kec_mbul_sf_utm1 <- st_transform(kec_mbul_sf, crs = 32648) Earthquakes_Pred_MBUL <- Earthquakes_Pred_dec %>% filter(x >= min(st_bbox(kec_mbul_sf_utm1)$xmin) & x <= max(st_bbox(kec_mbul_sf_utm1)$xmax) & y >= min(st_bbox(kec_mbul_sf_utm1)$ymin) & y <= max(st_bbox(kec_mbul_sf_utm1)$ymax)) min(Earthquakes_Pred_MBUL$pred_mean) max(Earthquakes_Pred_MBUL$pred_mean) # Menambahkan jitter pada koordinat x dan y sebelum interpolasi Earthquakes_Pred_MBUL$x <- jitter(Earthquakes_Pred_MBUL$x, amount = 0.0001) Earthquakes_Pred_MBUL$y <- jitter(Earthquakes_Pred_MBUL$y, amount = 0.0001) interpolated_mbul <- with(Earthquakes_Pred_MBUL, interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200)) interpolated_df_mbul <- as.data.frame(expand.grid(x = interpolated_mbul$x, y = interpolated_mbul$y)) interpolated_df_mbul$pred_mean <- as.vector(interpolated_mbul$z) interpolated_df_mbul <- na.omit(interpolated_df_mbul) # Menghapus nilai NA #Plot interpolated_df_mbul$title = "Mukomuko, Bengkulu Utara, & Lebong" mbul = ggplot() + geom_raster(data = interpolated_df_mbul, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) + scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) + geom_sf(data = kec_mbul_sf_utm1, fill = NA, color = "black", linewidth = 0.3) + geom_sf_text(data = kec_mbul, aes(label = KECAMATAN), size = 1.5, color = "black", fontface = "bold") + coord_sf(crs = 32648) + theme_bw() + theme( panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7), panel.grid.major = element_line(color = "grey95", linetype = "solid"), panel.grid.minor = element_line(color = "grey95", linetype = "solid"), ) + labs( x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude" ) + facet_grid(. ~ title) mbul
Bengkulu Selatan and Kaur
#== 2. Bengkulu Selatan and Kaur kec_bsk<-kec_cropped1[kec_cropped1$KECAMATAN %in%c("KAUR SELATAN","KAUR UTARA","KELAM TENGAH", "KAUR TENGAH","KINAL","LUNGKANG KULE","LUAS", "MAJE","MUARA SAHUNG","NASAL","PADANG GUCI HILIR", "PADANG GUCI HULU","SEMIDANG GUMAY","TANJUNG KEMUNING","TETAP", "KEDURANG","SEGINIM","PINO","MANNA","KOTA MANNA","PINO RAYA", "KEDURANG ILIR", "AIR NIPIS","ULU MANNA","BUNGA MAS", "PASAR MANNA"),] kec_bsk_sf <- st_as_sf(kec_bsk) kec_bsk_sf_utm1 <- st_transform(kec_bsk_sf, crs = 32648) Earthquakes_Pred_BSK <- Earthquakes_Pred_dec %>% filter(x >= min(st_bbox(kec_bsk_sf_utm1)$xmin) & x <= max(st_bbox(kec_bsk_sf_utm1)$xmax) & y >= min(st_bbox(kec_bsk_sf_utm1)$ymin) & y <= max(st_bbox(kec_bsk_sf_utm1)$ymax)) min(Earthquakes_Pred_BSK$pred_mean) max(Earthquakes_Pred_BSK$pred_mean) Earthquakes_Pred_BSK$x <- jitter(Earthquakes_Pred_BSK$x, amount = 0.0001) Earthquakes_Pred_BSK$y <- jitter(Earthquakes_Pred_BSK$y, amount = 0.0001) interpolated_BSK <- with(Earthquakes_Pred_BSK, interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200)) interpolated_df_BSK <- as.data.frame(expand.grid(x = interpolated_BSK$x, y = interpolated_BSK$y)) interpolated_df_BSK$pred_mean <- as.vector(interpolated_BSK$z) interpolated_df_BSK <- na.omit(interpolated_df_BSK) # Menghapus nilai NA #Plot interpolated_df_BSK$title = "Bengkulu Selatan and Kaur" bsk = ggplot() + geom_raster(data = interpolated_df_BSK, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) + scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) + geom_sf(data = kec_bsk_sf_utm1, fill = NA, color = "black", linewidth = 0.3) + geom_sf_text(data = kec_bsk, aes(label = KECAMATAN), size = 1.5, color = "black", fontface = "bold") + coord_sf(crs = 32648) + theme_bw() + theme( panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7), panel.grid.major = element_line(color = "grey95", linetype = "solid"), panel.grid.minor = element_line(color = "grey95", linetype = "solid"), ) + labs( x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude" ) + facet_grid(. ~ title) bsk
9.2.4.2 Sumatra Barat Province
Kepulauan Mentawai
#=== 1. Kab Mentawai kab_sumbar_wgs84 <- st_transform(kab_sumbar, crs = st_crs(kec)) bbox_kab_sumbar <- st_bbox(kab_sumbar_wgs84) kec_cropped <- st_crop(kec_valid, bbox_kab_sumbar) kec_mentawai<-kec_cropped[kec_cropped$KECAMATAN %in%c("PAGAI SELATAN","PAGAI UTARA","SIBERUT BARAT","SIBERUT BARAT DAYA","SIBERUT SELATAN", "SIBERUT UTARA","SIBERUT TENGAH","SIKAKAP","SIPORA SELATAN","SIPORA UTARA"),] kec_mentawai_sf <- st_as_sf(kec_mentawai) kec_mentawai_sf_utm1 <- st_transform(kec_mentawai_sf, crs = 32648) Earthquakes_Pred_Mentawai <- Earthquakes_Pred_dec %>% filter(x >= min(st_bbox(kec_mentawai_sf_utm1)$xmin) & x <= max(st_bbox(kec_mentawai_sf_utm1)$xmax) & y >= min(st_bbox(kec_mentawai_sf_utm1)$ymin) & y <= max(st_bbox(kec_mentawai_sf_utm1)$ymax)) #min(Earthquakes_Pred_Mentawai$pred_mean) #max(Earthquakes_Pred_Mentawai$pred_mean) Earthquakes_Pred_Mentawai$x <- jitter(Earthquakes_Pred_Mentawai$x, amount = 0.0001) Earthquakes_Pred_Mentawai$y <- jitter(Earthquakes_Pred_Mentawai$y, amount = 0.0001) interpolated_mentawai <- with(Earthquakes_Pred_Mentawai, interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200)) interpolated_df_mentawai <- as.data.frame(expand.grid(x = interpolated_mentawai$x, y = interpolated_mentawai$y)) interpolated_df_mentawai$pred_mean <- as.vector(interpolated_mentawai$z) interpolated_df_mentawai <- na.omit(interpolated_df_mentawai) interpolated_df_mentawai$title = "Kep. Mentawai" mentawai = ggplot() + geom_raster(data = interpolated_df_mentawai, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) + scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) + geom_sf(data = kec_mentawai_sf_utm1, fill = NA, color = "black", linewidth = 0.3) + geom_sf_text(data = kec_mentawai, aes(label = KECAMATAN), size = 1.5, color = "black", fontface = "bold") + coord_sf(crs = 32648) + theme_bw() + theme( panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7), panel.grid.major = element_line(color = "grey95", linetype = "solid"), panel.grid.minor = element_line(color = "grey95", linetype = "solid"), ) + labs( x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude" ) + facet_grid(. ~ title) mentawai
Pesisir Selatan
#=== 2. Pesisir Selatan kec_pesisir_selatan<-kec_cropped[kec_cropped$KECAMATAN %in%c("IV JURAI","IV NAGARI BAYANG UTARA","AIRPURA","BASA AMPEK BALAI TAPAN","BATANG KAPAS","BAYANG","KOTO XI TARUSAN","LINGGO SARI BAGANTI","LENGAYANG","LUNANG","PANCUNG SOAL","RANAH AMPEK HULU TAPAN","RANAH PESISIR","SILAUT","SUTERA"),] kec_pesisir_selatan_sf <- st_as_sf(kec_pesisir_selatan) kec_pesisir_selatan_sf_utm1 <- st_transform(kec_pesisir_selatan_sf, crs = 32648) Earthquakes_Pred_Pesisir_Selatan <- Earthquakes_Pred_dec %>% filter(x >= min(st_bbox(kec_pesisir_selatan_sf_utm1)$xmin) & x <= max(st_bbox(kec_pesisir_selatan_sf_utm1)$xmax) & y >= min(st_bbox(kec_pesisir_selatan_sf_utm1)$ymin) & y <= max(st_bbox(kec_pesisir_selatan_sf_utm1)$ymax)) #min(Earthquakes_Pred_Pesisir_Selatan$pred_mean) #max(Earthquakes_Pred_Pesisir_Selatan$pred_mean) Earthquakes_Pred_Pesisir_Selatan$x <- jitter(Earthquakes_Pred_Pesisir_Selatan$x, amount = 0.0001) Earthquakes_Pred_Pesisir_Selatan$y <- jitter(Earthquakes_Pred_Pesisir_Selatan$y, amount = 0.0001) interpolated_PS <- with(Earthquakes_Pred_Pesisir_Selatan, interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200)) interpolated_df_PS <- as.data.frame(expand.grid(x = interpolated_PS$x, y = interpolated_PS$y)) interpolated_df_PS$pred_mean <- as.vector(interpolated_PS$z) interpolated_df_PS <- na.omit(interpolated_df_PS) # Menghapus nilai NA interpolated_df_PS$title = "Pesisir Selatan" ps = ggplot() + geom_raster(data = interpolated_df_PS, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) + scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) + geom_sf(data = kec_pesisir_selatan_sf_utm1, fill = NA, color = "black", linewidth = 0.3) + geom_sf_text(data = kec_pesisir_selatan, aes(label = KECAMATAN), size = 1.5, color = "black", fontface = "bold") + coord_sf(crs = 32648) + theme_bw() + theme( panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7), panel.grid.major = element_line(color = "grey95", linetype = "solid"), panel.grid.minor = element_line(color = "grey95", linetype = "solid"), ) + labs( x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude" ) + facet_grid(. ~ title) ps
9.2.4.3 Sumatra Utara Province
Padang Lawas and Mandailing Natal
#=== 1. Padang lawas and Mandailing natal kab_sumut_wgs84 <- st_transform(kab_sumut, crs = st_crs(kec)) bbox_kab_sumut <- st_bbox(kab_sumut_wgs84) kec_cropped2 <- st_crop(kec_valid, bbox_kab_sumut) kec_plm<-kec_cropped2[kec_cropped2$KECAMATAN %in%c("AEK NABARA BARUMUN","BARUMUN","BARUMUN BARAT","BARUMUN BARU","BARUMUN SELATAN","BARUMUN TENGAH","BATANG LUBU SUTAM","HURISTAK","HUTA RAJA TINGGI","LUBUK BARUMUN","SIHAPAS BARUMUN","SOSA","SOSA JULU","SOSA TIMUR","SOSOPAN","ULU BARUMUN","ULU SOSA", "BATAHAN","BATANG NATAL","BUKIT MALINTANG","HUTA BARGOT","KOTANOPAN","LEMBAH SORIK MARAPI","LINGGA BAYU","MUARA BATANG GADIS","MUARA SIPONGI","NAGA JUANG","NATAL","PEKANTAN","PENYABUNGAN BARAT","PENYABUNGAN KOTA","PENYABUNGAN SELATAN","PENYABUNGAN TIMUR","PENYABUNGAN UTARA","PUNCAK SORIK MARAPI","RANTO BAEK","SIABU","SINUNUKAN","TAMBANGAN","ULU PUNGKUT"),] kec_plm_sf <- st_as_sf(kec_plm) kec_plm_sf_utm1 <- st_transform(kec_plm_sf, crs = 32648) Earthquakes_Pred_PLM <- Earthquakes_Pred_dec %>% filter(x >= min(st_bbox(kec_plm_sf_utm1)$xmin) & x <= max(st_bbox(kec_plm_sf_utm1)$xmax) & y >= min(st_bbox(kec_plm_sf_utm1)$ymin) & y <= max(st_bbox(kec_plm_sf_utm1)$ymax)) #min(Earthquakes_Pred_PLM$pred_mean) #max(Earthquakes_Pred_PLM$pred_mean) Earthquakes_Pred_PLM$x <- jitter(Earthquakes_Pred_PLM$x, amount = 0.0001) Earthquakes_Pred_PLM$y <- jitter(Earthquakes_Pred_PLM$y, amount = 0.0001) interpolated_PLM <- with(Earthquakes_Pred_PLM, interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200)) interpolated_df_PLM <- as.data.frame(expand.grid(x = interpolated_PLM$x, y = interpolated_PLM$y)) interpolated_df_PLM$pred_mean <- as.vector(interpolated_PLM$z) interpolated_df_PLM <- na.omit(interpolated_df_PLM) # Menghapus nilai NA interpolated_df_PLM$title = "Padang Lawas & Mandailing Natal" plm = ggplot() + geom_raster(data = interpolated_df_PLM, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) + scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) + geom_sf(data = kec_plm_sf_utm1, fill = NA, color = "black", linewidth = 0.3) + geom_sf_text(data = kec_plm, aes(label = KECAMATAN), size = 1.5, color = "black", fontface = "bold") + coord_sf(crs = 32648) + theme_bw() + theme( panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7), panel.grid.major = element_line(color = "grey95", linetype = "solid"), panel.grid.minor = element_line(color = "grey95", linetype = "solid"), ) + labs( x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude" ) + facet_grid(. ~ title) plm
Nias
#=== 2. Nias kec_nias<-kec_cropped2[kec_cropped2$KECAMATAN %in%c("AFULU","ALASA","ALASA TALUMUZOI","LAHEWA","LAHEWA TIMUR","LOTU","NAMOHALU","SAWO","SITOLU ORI","TUGALA OYO","TUHEMBERUA","GUNUNGSITOLI","GUNUNGSITOLI ALO'OA","GUNUNGSITOLI BARAT","GUNUNGSITOLI IDANOI","GUNUNGSITOLI SELATAN","GUNUNGSITOLI UTARA","BAWOLATO","BOTOMUZOI","GIDO","HILIDUHO","HILISERANGKAI","IDANOGAWO","MA'U","SOGAE'ADU","SOMOLO-MOLO","ULUGAWO","LAHOMI","LOLOFITU MOI", "MANDREHE","MANDREHE BARAT","MANDREHE UTARA","MORO'O","SORIMBU","UKU MORO'O"),] kec_nias_sf <- st_as_sf(kec_nias) kec_nias_sf_utm1 <- st_transform(kec_nias_sf, crs = 32648) Earthquakes_Pred_Nias <- Earthquakes_Pred_dec %>% filter(x >= min(st_bbox(kec_nias_sf_utm1)$xmin) & x <= max(st_bbox(kec_nias_sf_utm1)$xmax) & y >= min(st_bbox(kec_nias_sf_utm1)$ymin) & y <= max(st_bbox(kec_nias_sf_utm1)$ymax)) #min(Earthquakes_Pred_Nias$pred_mean) #max(Earthquakes_Pred_Nias$pred_mean) Earthquakes_Pred_Nias$x <- jitter(Earthquakes_Pred_Nias$x, amount = 0.0001) Earthquakes_Pred_Nias$y <- jitter(Earthquakes_Pred_Nias$y, amount = 0.0001) interpolated_Nias <- with(Earthquakes_Pred_Nias, interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200)) interpolated_df_Nias <- as.data.frame(expand.grid(x = interpolated_Nias$x, y = interpolated_Nias$y)) interpolated_df_Nias$pred_mean <- as.vector(interpolated_Nias$z) interpolated_df_Nias<- na.omit(interpolated_df_Nias) # Menghapus nilai NA interpolated_df_Nias$title = "Nias" nias = ggplot() + geom_raster(data = interpolated_df_Nias, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) + scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) + geom_sf(data = kec_nias_sf_utm1, fill = NA, color = "black", linewidth = 0.3) + geom_sf_text(data = kec_nias, aes(label = KECAMATAN), size = 1.5, color = "black", fontface = "bold") + coord_sf(crs = 32648) + theme_bw() + theme( panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7), panel.grid.major = element_line(color = "grey95", linetype = "solid"), panel.grid.minor = element_line(color = "grey95", linetype = "solid"), ) + labs( x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude" ) + facet_grid(. ~ title) nias
Nias Selatan
#=== 3. Nias Selatan kec_nias_sel<-kec_cropped2[kec_cropped2$KECAMATAN %in%c("AMANDRAYA","ARAMO","BORONADU","FANAYAMA","GOMO","HIBALA","HILIMEGAI","HILISALAWA'AHE","HURUNA","IDANOTAE","LAHUSA","LOLOMATUA","LOLOWAU","LUAHAGUNDRE MANIAMOLO","MANIAMOLO","MAZINO","MAZO","O'O'U","ONOHAZUMBA","ONOLALU","PULAU-PULAU BATU","PULAU-PULAU BARAT","PULAU-PULAU BATU TIMUR", "PULAU-PULAU BATU UTARA","SIDUA'ORI","SIMUK", "SOMAMBAWA","SUSUA","TANAH MASA","TOMA","ULUNOYO", "ULU IDANOTAE","UMBUNASI","ULUSUSUA"),] kec_niassel_sf <- st_as_sf(kec_nias_sel) kec_niassel_sf_utm1 <- st_transform(kec_niassel_sf, crs = 32648) Earthquakes_Pred_Nias_Sel <- Earthquakes_Pred_dec %>% filter(x >= min(st_bbox(kec_niassel_sf_utm1)$xmin) & x <= max(st_bbox(kec_niassel_sf_utm1)$xmax) & y >= min(st_bbox(kec_niassel_sf_utm1)$ymin) & y <= max(st_bbox(kec_niassel_sf_utm1)$ymax)) Earthquakes_Pred_Nias_Sel$x <- jitter(Earthquakes_Pred_Nias_Sel$x, amount = 0.0001) Earthquakes_Pred_Nias_Sel$y <- jitter(Earthquakes_Pred_Nias_Sel$y, amount = 0.0001) interpolated_Nias_Sel <- with(Earthquakes_Pred_Nias_Sel, interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200)) interpolated_df_Nias_Sel <- as.data.frame(expand.grid(x = interpolated_Nias_Sel$x, y = interpolated_Nias_Sel$y)) interpolated_df_Nias_Sel$pred_mean <- as.vector(interpolated_Nias_Sel$z) interpolated_df_Nias_Sel<- na.omit(interpolated_df_Nias_Sel) # Menghapus nilai NA interpolated_df_Nias_Sel$title = "Nias Selatan" nias_sel = ggplot() + geom_raster(data = interpolated_df_Nias_Sel, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) + scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) + geom_sf(data = kec_niassel_sf_utm1, fill = NA, color = "black", linewidth = 0.3) + geom_sf_text(data = kec_nias_sel, aes(label = KECAMATAN), size = 1.5, color = "black", fontface = "bold") + coord_sf(crs = 32648) + theme_bw() + theme( panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7), panel.grid.major = element_line(color = "grey95", linetype = "solid"), panel.grid.minor = element_line(color = "grey95", linetype = "solid"), ) + labs( x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude" ) + facet_grid(. ~ title) nias_sel
9.3 High-Resolution of Exceedance Probability
9.3.1 Overall: By Month
ggplot() +
geom_tile(data = Earthquakes_Pred, aes(x = x, y = y, fill = pred_prob)) +
geom_contour(data = Earthquakes_Pred, aes(x = x, y = y, z = pred_mean, colour = after_stat(level)),
linewidth = 0.2, colour = "red", alpha = 0.5) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red")) +
facet_wrap(~ Time, labeller = as_labeller(Month_names), ncol = 4) +
theme_bw() +
ylab("") +
xlab("") +
geom_sf(data = Sumatra_sf_utm1, fill = NA, size = 0.01) +
theme(legend.position = "bottom", text = element_text(size = 19)) +
guides(fill = guide_colorbar(title.position = "left", title.vjust = 1,
frame.colour = "black",
barwidth = 20,
barheight = 1.5)) +
theme(axis.text.x = element_blank(), #remove x axis labels
axis.ticks.x = element_blank(), #remove x axis ticks
axis.text.y = element_blank(), #remove y axis labels
axis.ticks.y = element_blank() #remove y axis ticks
) +
labs(fill = "Earthquakes \nExceedance Probability") +
ggtitle("Earthquakes Exceedance Probability")
9.3.2 Overall: By Province
max_prob = max(Earthquakes_Pred$pred_prob); max_prob
ggplot() +
geom_tile(data = Earthquakes_Pred_dec, aes(x = x, y = y, fill = pred_prob)) +
geom_contour(data = Earthquakes_Pred_dec, aes(x = x, y = y, z = pred_mean, colour = stat(level)),
size = 0.2, colour = "red", alpha = 0.5) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red"),limits = c(0, max_prob)) +
theme_bw() +
ylab("Latitude") +
xlab("Longitude") +
geom_sf(data = Sumatra_sf_utm1, fill = NA, size = 0.01) +
geom_sf_text(data = Sumatra_sf_utm1, aes(label = NAME_1), size = 3, color = "black", fontface = "bold") +
#geom_sf(data = patahan_sumatera, color = "blue", size = 1) +
theme(legend.position = "bottom",
text = element_text(size = 12)) +
theme(
panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
plot.title = element_text(hjust = 0.5, size = 12), # Menyesuaikan ukuran dan posisi judul
panel.grid = element_blank() # Menghilangkan grid
)+
guides(fill = guide_colorbar(title.position = "left", title.vjust = 1,
frame.colour = "black",
barwidth = 20,
barheight = 1.5))+
facet_grid(. ~ title)
9.4 Interpolation
9.4.1 Aceh
Sumatra_Aceh <- Sumatra_sf_utm1 %>% filter(NAME_1 == 'Aceh')
Sumatra_Aceh_wgs84 <- st_transform(Sumatra_Aceh, crs = st_crs(kab))
bbox_Sumatra_Aceh <- st_bbox(Sumatra_Aceh_wgs84)
kab_cropped <- st_crop(kab_valid, bbox_Sumatra_Aceh)
kab_aceh<-kab_cropped[kab_cropped$KAB_KOTA %in%c("ACEH BARAT","ACEH BARAT DAYA","ACEH BESAR",
"ACEH JAYA","ACEH SELATAN","ACEH SINGKIL","ACEH TAMIANG","ACEH TENGAH","ACEH TENGGARA","ACEH TIMUR","ACEH UTARA","BENER MERIAH","BIREUEN","GAYO LUES","NAGAN RAYA","PIDIE","PIDIE JAYA","SIMEULUE","BANDA ACEH","KOTA LANGSA","KOTA LHOKSEUMAWE","KOTA SABANG","KOTA SUBULUSSALAM"),]
kab_aceh_sf <- st_as_sf(kab_aceh)
kab_aceh_sf_utm1 <- st_transform(kab_aceh_sf, crs = 32648)
Earthquakes_Pred_Aceh <- Earthquakes_Pred_dec %>%
filter(x >= min(st_bbox(kab_aceh_sf_utm1)$xmin) & x <= max(st_bbox(kab_aceh_sf_utm1)$xmax) &
y >= min(st_bbox(kab_aceh_sf_utm1)$ymin) & y <= max(st_bbox(kab_aceh_sf_utm1)$ymax))
min(Earthquakes_Pred_Aceh$pred_mean)
max(Earthquakes_Pred_Aceh$pred_mean)
Earthquakes_Pred_Aceh$x <- jitter(Earthquakes_Pred_Aceh$x, amount = 0.0001)
Earthquakes_Pred_Aceh$y <- jitter(Earthquakes_Pred_Aceh$y, amount = 0.0001)
interpolated_aceh <- with(Earthquakes_Pred_Aceh,
interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200))
interpolated_df_aceh <- as.data.frame(expand.grid(x = interpolated_aceh$x, y = interpolated_aceh$y))
interpolated_df_aceh$pred_mean <- as.vector(interpolated_aceh$z)
interpolated_df_aceh <- na.omit(interpolated_df_aceh) # Menghapus nilai NA
interpolated_df_aceh$title = "Aceh"
aceh = ggplot() +
geom_raster(data = interpolated_df_aceh, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) +
geom_sf(data = kab_aceh_sf_utm1, fill = NA, color = "black", linewidth = 0.3) +
geom_sf_text(data = kab_aceh, aes(label = KAB_KOTA), size = 2.5, color = "black", fontface = "bold") +
coord_sf(crs = 32648) +
theme_bw() +
theme(
panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
) +
labs(
x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude"
) +
facet_grid(. ~ title)
aceh
9.4.2 Riau
Sumatra_Riau <- Sumatra_sf_utm1 %>% filter(NAME_1 == 'Riau')
Sumatra_Riau_wgs84 <- st_transform(Sumatra_Riau, crs = st_crs(kab))
bbox_Sumatra_Riau <- st_bbox(Sumatra_Riau_wgs84)
# Crop kab berdasarkan bounding box dari Sumatra_Bengkulu
kab_cropped <- st_crop(kab_valid, bbox_Sumatra_Riau)
kab_riau<-kab_cropped[kab_cropped$KAB_KOTA %in%c("BENGKALIS","INDRAGIRI HILIR","INDRAGIRI HULU","KAMPAR","KEPULAUAN MERANTI", "KUANTAN SINGINGI","PELALAWAN","ROKAN HILIR","SIAK","KOTA DUMAI","KOTA PEKANBARU"),]
kab_riau_sf <- st_as_sf(kab_riau)
kab_riau_sf_utm1 <- st_transform(kab_riau_sf, crs = 32648)
Earthquakes_Pred_Riau <- Earthquakes_Pred_dec %>%
filter(x >= min(st_bbox(kab_riau_sf_utm1)$xmin) & x <= max(st_bbox(kab_riau_sf_utm1)$xmax) &
y >= min(st_bbox(kab_riau_sf_utm1)$ymin) & y <= max(st_bbox(kab_riau_sf_utm1)$ymax))
#min earthquake
min(Earthquakes_Pred_Riau$pred_mean)
max(Earthquakes_Pred_Riau$pred_mean)
#==== plot dengan ada lang lot dan smoothness
# Menambahkan jitter pada koordinat x dan y sebelum interpolasi
Earthquakes_Pred_Riau$x <- jitter(Earthquakes_Pred_Riau$x, amount = 0.0001)
Earthquakes_Pred_Riau$y <- jitter(Earthquakes_Pred_Riau$y, amount = 0.0001)
# Melakukan interpolasi setelah jitter
interpolated_riau <- with(Earthquakes_Pred_Riau,
interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200))
# Mengubah hasil interpolasi menjadi format data frame yang bisa di-plot dengan ggplot
interpolated_df_riau <- as.data.frame(expand.grid(x = interpolated_riau$x, y = interpolated_riau$y))
interpolated_df_riau$pred_mean <- as.vector(interpolated_riau$z)
interpolated_df_riau <- na.omit(interpolated_df_riau) # Menghapus nilai NA
#Plot
interpolated_df_riau$title = "Riau"
riau = ggplot() +
geom_raster(data = interpolated_df_riau, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) +
geom_sf(data = kab_riau_sf_utm1, fill = NA, color = "black", linewidth = 0.3) +
geom_sf_text(data = kab_riau, aes(label = KAB_KOTA), size = 2.5, color = "black", fontface = "bold") +
coord_sf(crs = 32648) +
theme_bw() +
theme(
panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
) +
labs(
x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude"
) +
facet_grid(. ~ title)
riau
9.4.3 Kepulauan Riau
Sumatra_kRiau <- Sumatra_sf_utm1 %>% filter(NAME_1 == 'Kepulauan Riau')
Sumatra_kRiau_wgs84 <- st_transform(Sumatra_kRiau, crs = st_crs(kab))
bbox_Sumatra_kRiau <- st_bbox(Sumatra_kRiau_wgs84)
# Crop kab berdasarkan bounding box dari Sumatra_Bengkulu
kab_cropped <- st_crop(kab_valid, bbox_Sumatra_kRiau)
kab_kriau<-kab_cropped[kab_cropped$KAB_KOTA %in%c("BINTAN","KARIMUN","KEPULAUAN ANAMBAS","LINGGA","NATUNA","KOTA BATAM","KOTA TANJUNGPINANG"),]
kab_kriau_sf <- st_as_sf(kab_kriau)
kab_kriau_sf_utm1 <- st_transform(kab_kriau_sf, crs = 32648)
Earthquakes_Pred_kRiau <- Earthquakes_Pred_dec %>%
filter(x >= min(st_bbox(kab_kriau_sf_utm1)$xmin) & x <= max(st_bbox(kab_kriau_sf_utm1)$xmax) &
y >= min(st_bbox(kab_kriau_sf_utm1)$ymin) & y <= max(st_bbox(kab_kriau_sf_utm1)$ymax))
#min earthquake
min(Earthquakes_Pred_kRiau$pred_mean)
max(Earthquakes_Pred_kRiau$pred_mean)
#==== plot dengan ada lang lot dan smoothness
# Menambahkan jitter pada koordinat x dan y sebelum interpolasi
Earthquakes_Pred_kRiau$x <- jitter(Earthquakes_Pred_kRiau$x, amount = 0.0001)
Earthquakes_Pred_kRiau$y <- jitter(Earthquakes_Pred_kRiau$y, amount = 0.0001)
# Melakukan interpolasi setelah jitter
interpolated_kriau <- with(Earthquakes_Pred_kRiau,
interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200))
# Mengubah hasil interpolasi menjadi format data frame yang bisa di-plot dengan ggplot
interpolated_df_kriau <- as.data.frame(expand.grid(x = interpolated_kriau$x, y = interpolated_kriau$y))
interpolated_df_kriau$pred_mean <- as.vector(interpolated_kriau$z)
interpolated_df_kriau <- na.omit(interpolated_df_kriau) # Menghapus nilai NA
#Plot
interpolated_df_kriau$title = "Kepulauan Riau"
kriau = ggplot() +
geom_raster(data = interpolated_df_kriau, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) +
geom_sf(data = kab_kriau_sf_utm1, fill = NA, color = "black", linewidth = 0.3) +
geom_sf_text(data = kab_kriau, aes(label = KAB_KOTA), size =1.5, color = "black", fontface = "bold") +
coord_sf(crs = 32648) +
theme_bw() +
theme(
panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
) +
labs(
x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude"
) +
facet_grid(. ~ title)
kriau
9.4.4 Jambi
Sumatra_jambi <- Sumatra_sf_utm1 %>% filter(NAME_1 == 'Jambi')
Sumatra_jambi_wgs84 <- st_transform(Sumatra_jambi, crs = st_crs(kab))
bbox_Sumatra_jambi <- st_bbox(Sumatra_jambi_wgs84)
kab_cropped <- st_crop(kab_valid, bbox_Sumatra_jambi)
kab_jambi<-kab_cropped[kab_cropped$KAB_KOTA %in%c("BATANGHARI","BUNGO","KERINCI","MERANGIN","MUARO JAMBI", "SAROLANGUN","TANJUNG JABUNG BARAT","TANJUNG JABUNG TIMUR","TEBO","KOTA JAMBI","KOTA SUNGAI PENUH"),]
kab_jambi_sf <- st_as_sf(kab_jambi)
kab_jambi_sf_utm1 <- st_transform(kab_jambi_sf, crs = 32648)
Earthquakes_Pred_jambi <- Earthquakes_Pred_dec %>%
filter(x >= min(st_bbox(kab_jambi_sf_utm1)$xmin) & x <= max(st_bbox(kab_jambi_sf_utm1)$xmax) &
y >= min(st_bbox(kab_jambi_sf_utm1)$ymin) & y <= max(st_bbox(kab_jambi_sf_utm1)$ymax))
#min earthquake
min(Earthquakes_Pred_jambi$pred_mean)
max(Earthquakes_Pred_jambi$pred_mean)
#==== plot dengan ada lang lot dan smoothness
# Menambahkan jitter pada koordinat x dan y sebelum interpolasi
Earthquakes_Pred_jambi$x <- jitter(Earthquakes_Pred_jambi$x, amount = 0.0001)
Earthquakes_Pred_jambi$y <- jitter(Earthquakes_Pred_jambi$y, amount = 0.0001)
# Melakukan interpolasi setelah jitter
interpolated_jambi <- with(Earthquakes_Pred_jambi,
interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200))
# Mengubah hasil interpolasi menjadi format data frame yang bisa di-plot dengan ggplot
interpolated_df_jambi <- as.data.frame(expand.grid(x = interpolated_jambi$x, y = interpolated_jambi$y))
interpolated_df_jambi$pred_mean <- as.vector(interpolated_jambi$z)
interpolated_df_jambi <- na.omit(interpolated_df_jambi) # Menghapus nilai NA
#Plot
interpolated_df_jambi$title = "Jambi"
jambi = ggplot() +
geom_raster(data = interpolated_df_jambi, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) +
geom_sf(data = kab_jambi_sf_utm1, fill = NA, color = "black", linewidth = 0.3) +
geom_sf_text(data = kab_jambi, aes(label = KAB_KOTA), size =1.5, color = "black", fontface = "bold") +
coord_sf(crs = 32648) +
theme_bw() +
theme(
panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
) +
labs(
x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude"
) +
facet_grid(. ~ title)
jambi
9.4.5 Lampung
Sumatra_lampung <- Sumatra_sf_utm1 %>% filter(NAME_1 == 'Lampung')
Sumatra_lampung_wgs84 <- st_transform(Sumatra_lampung, crs = st_crs(kab))
bbox_Sumatra_lampung <- st_bbox(Sumatra_lampung_wgs84)
kab_cropped <- st_crop(kab_valid, bbox_Sumatra_lampung)
kab_lampung<-kab_cropped[kab_cropped$KAB_KOTA %in%c("LAMPUNG BARAT","LAMPUNG SELATAN","LAMPUNG TENGAH","LAMPUNG TIMUR","LAMPUNG UTARA","MESUJI","PESAWARAN","PESISIR BARAT","PRINGSEWU","TANGGAMUS","TULANG BAWANG", "TULANG BAWANG BARAT","WAY KANAN","KOTA BANDAR LAMPUNG","KOTA METRO"),]
kab_lampung_sf <- st_as_sf(kab_lampung)
kab_lampung_sf_utm1 <- st_transform(kab_lampung_sf, crs = 32648)
Earthquakes_Pred_lampung <- Earthquakes_Pred_dec %>%
filter(x >= min(st_bbox(kab_lampung_sf_utm1)$xmin) & x <= max(st_bbox(kab_lampung_sf_utm1)$xmax) &
y >= min(st_bbox(kab_lampung_sf_utm1)$ymin) & y <= max(st_bbox(kab_lampung_sf_utm1)$ymax))
#min earthquake
min(Earthquakes_Pred_lampung$pred_mean)
max(Earthquakes_Pred_lampung$pred_mean)
#==== plot dengan ada lang lot dan smoothness
# Menambahkan jitter pada koordinat x dan y sebelum interpolasi
Earthquakes_Pred_lampung$x <- jitter(Earthquakes_Pred_lampung$x, amount = 0.0001)
Earthquakes_Pred_lampung$y <- jitter(Earthquakes_Pred_lampung$y, amount = 0.0001)
# Melakukan interpolasi setelah jitter
interpolated_lampung <- with(Earthquakes_Pred_lampung,
interp(x = x, y = y, z = pred_mean, nx = 20, ny = 20))
# Mengubah hasil interpolasi menjadi format data frame yang bisa di-plot dengan ggplot
interpolated_df_lampung <- as.data.frame(expand.grid(x = interpolated_lampung$x, y = interpolated_lampung$y))
interpolated_df_lampung$pred_mean <- as.vector(interpolated_lampung$z)
interpolated_df_lampung<- na.omit(interpolated_df_lampung) # Menghapus nilai NA
#Plot
interpolated_df_lampung$title = "Lampung"
lampung = ggplot() +
geom_raster(data = interpolated_df_lampung, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) +
geom_sf(data = kab_lampung_sf_utm1, fill = NA, color = "black", linewidth = 0.3) +
geom_sf_text(data = kab_lampung, aes(label = KAB_KOTA), size =1.5, color = "black", fontface = "bold") +
coord_sf(crs = 32648) +
theme_bw() +
theme(
panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
) +
labs(
x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude"
) +
facet_grid(. ~ title)
lampung
9.4.6 Sumatra Selatan
Sumatra_sumsel <- Sumatra_sf_utm1 %>% filter(NAME_1 == 'Sumatera Selatan')
Sumatra_sumsel_wgs84 <- st_transform(Sumatra_sumsel, crs = st_crs(kab))
bbox_Sumatra_sumsel <- st_bbox(Sumatra_sumsel_wgs84)
kab_cropped <- st_crop(kab_valid, bbox_Sumatra_sumsel)
kab_sumsel<-kab_cropped[kab_cropped$KAB_KOTA %in%c("BANYUASIN","EMPAT LAWANG","LAHAT","MUARA ENIM","MUSI BANYUASIN","MUSI RAWAS","MUSI RAWAS UTARA","OGAN ILIR","OGAN KOMERING ILIR", "OGAN KOMERING ULU","OGAN KOMERING ULU SELATAN","OGAN KOMERING ULU TIMUR","PENUKAL ABAB LEMATANG ILIR","KOTA LUBUK LINGGAU","KOTA PAGARALAM","KOTA PALEMBANG","KOTA PRABUMULIH"),]
kab_sumsel_sf <- st_as_sf(kab_sumsel)
kab_sumsel_sf_utm1 <- st_transform(kab_sumsel_sf, crs = 32648)
Earthquakes_Pred_sumsel <- Earthquakes_Pred_dec %>%
filter(x >= min(st_bbox(kab_sumsel_sf_utm1)$xmin) & x <= max(st_bbox(kab_sumsel_sf_utm1)$xmax) &
y >= min(st_bbox(kab_sumsel_sf_utm1)$ymin) & y <= max(st_bbox(kab_sumsel_sf_utm1)$ymax))
#min earthquake
min(Earthquakes_Pred_sumsel$pred_mean)
max(Earthquakes_Pred_sumsel$pred_mean)
#==== plot dengan ada lang lot dan smoothness
# Menambahkan jitter pada koordinat x dan y sebelum interpolasi
Earthquakes_Pred_sumsel$x <- jitter(Earthquakes_Pred_sumsel$x, amount = 0.0001)
Earthquakes_Pred_sumsel$y <- jitter(Earthquakes_Pred_sumsel$y, amount = 0.0001)
# Melakukan interpolasi setelah jitter
interpolated_sumsel <- with(Earthquakes_Pred_sumsel,
interp(x = x, y = y, z = pred_mean, nx = 200, ny = 200))
# Mengubah hasil interpolasi menjadi format data frame yang bisa di-plot dengan ggplot
interpolated_df_sumsel <- as.data.frame(expand.grid(x = interpolated_sumsel$x, y = interpolated_sumsel$y))
interpolated_df_sumsel$pred_mean <- as.vector(interpolated_sumsel$z)
interpolated_df_sumsel <- na.omit(interpolated_df_sumsel) # Menghapus nilai NA
#Plot
interpolated_df_sumsel$title = "Sumatra Selatan"
sumsel = ggplot() +
geom_raster(data = interpolated_df_sumsel, aes(x = x, y = y, fill = pred_mean), interpolate = TRUE) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red"), limits = c(0, max_pred)) +
geom_sf(data = kab_sumsel_sf_utm1, fill = NA, color = "black", linewidth = 0.3) +
geom_sf_text(data = kab_sumsel, aes(label = KAB_KOTA), size =1.5, color = "black", fontface = "bold") +
coord_sf(crs = 32648) +
theme_bw() +
theme(
panel.background = element_rect(fill = "white", color = "black", linewidth = 0.7),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
) +
labs(
x = "Longitude", y = "Latitude", fill = "Earthquake Magnitude"
) +
facet_grid(. ~ title)
sumsel