Get various climate and terrain data files for Korea, do some things in R

korean-aster-raster-data.R

# we have downloaded the files from NASA Land Processes Distributed Active
# Archive Center (LP DAAC) at https://earthexplorer.usgs.gov/
# and unzipped the downloaded file
# which contains many zip files...

# read in the Zipped files
setwd("C:/Users/bmarwick/Downloads/030432578418876/DEM-ASTER-Korea")
library(fs)
zip_files <- dir_ls(regexp = "[.]zip$")

# unzip them
library(tidyverse)
map(zip_files, unzip)

# read in the TIF files
library(rgdal)
library(raster)

raster_files <- dir_ls(regexp = "dem.tif$")

# import rasters into a list
many_rasters <- map(raster_files, raster)

# make a mosaic of the many rasters in a list 
names(many_rasters) <- NULL
many_rasters$fun <- mean
mos <- do.call(mosaic, many_rasters)

# draw a plot
plot(mos)

# save combined raster
writeRaster(mos, "ASTER-raster-Korea-raster")

korean-earthenv-raster-data.R

# files from http://www.earthenv.org/DEM
# ee <- dir_ls('C:/Users/bmarwick/Downloads/EarthEnv-DEM90_N35E125', regexp = "bil$")
ee <- list.files('C:/Users/bmarwick/Downloads/EarthEnv-DEM90_N35E125', 
								 pattern = "bil$",
								 full.names = TRUE)

ee_raster <- raster(ee)

plot(ee_raster)

korean-site-ages-calibrated.R

dates <- read.table(text = "41500±1500
40690±1500
39280±470
38500±1000
37300±400
36070±380
31700±900
31200±900
31000±600
30690±3000
32000±800
30800±100
27500±300
29800±1500
29900±100
29760±300")

library(tidyverse)
dates_1 <- 
dates %>% 
	separate(V1, into = c('age', 'error')) %>% 
	mutate_all(as.numeric)

library(Bchron)
ages1 = BchronCalibrate(ages=dates_1$age,
												ageSds=dates_1$error,
												calCurves=rep('intcal13', nrow(dates_1)))

plot(ages1)

# if we want at  want a credible interval (i.e. a single contiguous range) 			
# First create age samples for each date
age_samples = SampleAges(ages1)
# Now summarise them with quantile - this gives a 95% credible interval
apply(age_samples, 2, quantile, prob=c(0.025,0.975))


ages1Dens = BchronDensity(ages=dates_1$age,
												 ageSds=dates_1$error,
												 calCurves=rep('intcal13', nrow(dates_1)))

#  output the possible start/end dates of phases:

summary(ages1Dens, prob = 0.95)

plot(ages1Dens, 
									xlab='Age (cal years BP)', 
									xaxp=c(0, 50000, 16))

plot(ages1Dens$ageGrid,
		 ages1Dens$densities,
		 type='l')

korean_site_age_densities <- 
ggplot(data_frame(ageGrid = ages1Dens$ageGrid/1000,
									densities = ages1Dens$densities[,1]), 
			 aes(ageGrid, 
			 		densities)) +
	geom_line() +
	scale_x_continuous(limits = c(20, 70),
										 name = "x 1000 years ago") 

library(gsloid)
# subset the MIS data for the last 250 ka years
mis_last_70ka <- LR04_MISboundaries[LR04_MISboundaries$LR04_Age_ka_start %in% 20:70, ]

delta_oxy <- 
ggplot() +
	annotate("rect", 
					 xmin = mis_last_70ka$LR04_Age_ka_end, 
					 xmax = mis_last_70ka$LR04_Age_ka_start, 
					 ymin = -Inf, 
					 ymax = Inf,
					 alpha = .2,
					 fill = c("grey70", "white")) +
	annotate("text", 
					 label = mis_last_70ka$label_MIS, 
					 x =     mis_last_70ka$LR04_Age_ka_mid,
					 y = 		seq(3.1, 2.8, length.out = 2),
					 size = 3) +
	geom_line(data = lisiecki2005,                # add d18O line
						aes(Time,
								d18O)) +  
	scale_x_continuous(limits = c(20, 70),
										 name = "x 1000 years ago") +
	scale_y_reverse(name = bquote(delta^18*O)) +
	theme_bw()

library(cowplot)
plot_grid(delta_oxy, korean_site_age_densities, ncol = 1, align = "v", axis = 'b')

plot-TraCE-data-after-Barton-et-al.R

# We follow the code from Barton et al

library(ncdf4) # import GCM data
library(rgdal) # read GCM data
library(raster) # process GCM data
library(rasterVis) # plotting GCM data
library(tidyverse) # data management and plotting
library(magrittr) # pipes for code readability
library(EMD) # calculate trends in the data
library(dismo) # for latitudinally weighted samples
library(mgcv) # fit GAM for downscaling

# this is Korea
bbox <- extent(c(33.5, 39, 124, 130))
plot(bbox)

                               # approx locations only!
samp_pts <- matrix(c( 129, 35, # Busan
										  127, 36, # Daejeon
										  127, 37, # Seoul
										  128, 37), # Pyeongchang
									 ncol = 2, byrow = T)

# see a map of Korea
library(ggmap)
get_map(location = c(left = 124, 
										 bottom =33.5, 
										 right = 130, 
										 top = 39)) %>% 
	ggmap() + 
	labs(x = 'Longitude', y = 'Latitude') + 
	geom_point(data = as_data_frame(samp_pts), 
						 aes(x = V1, y = V2), 
						 size = 3, 
						 color = 'red')


# Create a matrix with the coordinates for the locations of interest 
# just some random points in Korea...

samp_pts <- matrix(c( 129, 35, # Busan
										  127, 36, # Daejeon
										  127, 37, # Seoul
										  128, 37), # Pyeongchang
									 ncol = 2, byrow = T)

plot(samp_pts)

# Now pull all the TraCE data into one data frame, with one row per year, and one column per variable/location combination. First rbind the two sets of TraCE data and transpose the results, turning the 6 rows into 6 columns. Add a column for the Year (in ka BP), and use to select only the entries earlier than 6,000 BP.

setwd("C:/Users/bmarwick/Downloads/Korean-palaeoclimate-data")

trace_dat <- rbind(
	brick('TraCE-21k/trace.01-36.22000BP.cam2.TREFHT.22000BP_decavg_400BCE.nc') %>%
		raster::extract(samp_pts) %>%
		subtract(273.15), # convert from kelvin to C
	brick('TraCE-21k/trace.01-36.22000BP.cam2.PRECT.22000BP_decavg_400BCE.nc') %>%
		raster::extract(samp_pts) %>% # extract data at these coordinates
		multiply_by(3.154e+10)) %>% # convert to mm/year
	t  %>% # transpose
	as.data.frame %>%
	set_colnames(c('tmp,Busan', 'tmp,Daejeon', 'tmp,Seoul', 'tmp,Pyeongchang',
								 'prc,Busan', 'prc,Daejeon', 'prc,Seoul', 'prc,Pyeongchang')) %>%
	rownames_to_column('Year') %>%
	mutate(Year = as.numeric(substring(Year, 3))) %>%
	filter(Year > 6) # get all the decades up to 6ka BP


# Trend Analysis

# Now organize the temperature and precipitation data to make plotting easier using functions from tidyr, and use the EMD package to calculate trend lines using the empirical mode decomposition approach.

trace_plot <- 
	trace_dat %>%
	gather(key, value, - Year) %>%
	separate(key, c('Variable', 'Region'), ',') %>%
	mutate(Region = factor(Region, levels = c('Busan', 'Daejeon',
																						'Seoul', 'Pyeongchang')),
				 Variable = ifelse(Variable == 'tmp', 'Temperature (°C)', 'Precipitation (mm)'))

emd_res <- function(x) emd(x, boundary = 'wave')$residue

trace_emd <- 
	trace_dat %>%
	mutate_at(vars(-Year), emd_res) %>%
	gather(key, value, - Year) %>%
	separate(key, c('Variable', 'Region'), ',') %>%
	mutate(Region = factor(Region, levels = c('Busan', 'Daejeon',
																						'Seoul', 'Pyeongchang')),
				 Variable = ifelse(Variable == 'tmp',
				 									'Temperature (°C)', 'Precipitation (mm)'))
                          
 # now plot
 
 ggplot(data = trace_plot, 
			 aes(x = Year, y = value)) +
	facet_grid(Variable ~ Region, 
						 switch = 'y', 
						 scale = 'free_y') +
	geom_vline(xintercept = c(22, 19, 14, 10, 6), 
						 lty = 2) +
	geom_point(aes(color = Variable), 
						 alpha = .3) +
	geom_line(data = trace_emd, 
						size = 1.2, 
						color = "black", 
						alpha = .8) +
	scale_x_reverse(breaks = seq(6,22,4)) +
	labs(x = '1,000 Years BP', y = '') +
	guides(color = "none") +
	theme_bw(base_size = 20) +
	theme(strip.background = element_blank())