I decided to add some more parameters to my model of agricultural conditions. I downloaded data from here. This study represents summer climatic reconstruction data from AD 755 to the present in Europe. Because it is only available for download as an .nc file, I had to use R to read the file and reduce it down to only the data that I’m interested in. I had to select the grid squares corresponding to England and then select the time periods of interest
library(lattice)
library(ncdf4)
library(RColorBrewer)
library(chron)
library(ggplot2)
setwd("~/Grad Year 3/Advanced Data Structures/")
I have never used a file with a .nc extension, so I found the tutorial here extremely helpful. About halfway through using the code from the tutorial, I had to start writing my own code because the tutorial code does not account for NA values in the data set (the steps in the tutorial convert an array to a long vector, but the conversion function eliminates NAs. As a result, my dataframe created from the instructions in the tutorial did not have the correct dimensions or placement of values). I am very grateful to this tutorial, however, because I had no idea how to read a .nc file without it!
ncname <- "BHM_Eu2k2_Recon"
ncfname <- "spatial data/BHM_Eu2k2_Recon.nc"
#select the variable name from the .nc file that I am using
dname <- "TT"
ncin <- nc_open(ncfname)
#extracts longitude values from the .nc file
lon <- ncvar_get(ncin, "lon")
nlon <- dim(lon)
#extracts latitude values from the .nc file
lat <- ncvar_get(ncin, "lat", verbose = F)
nlat <- dim(lat)
#extracts the year values from the .nc file
t <- ncvar_get(ncin, "time")
tunits <- ncatt_get(ncin, "time", "units")
nt <- dim(t)
#creates a three dimenionsal array with the extracted variables
tmp.array <- ncvar_get(ncin, dname)
dlname <- ncatt_get(ncin, dname, "long_name")
dunits <- ncatt_get(ncin, dname, "units")
fillvalue <- ncatt_get(ncin, dname, "_FillValue")
#closes the connection to the .nc file
nc_close(ncin)
#adjusts the default NA setting from "fillvalue" (standard for the .nc file) to "NA" (which R can understand)
tmp.array[tmp.array == fillvalue$value] <- NA
Where I left off from the tutorial and wrote my own code:
#select only the medieval years and grid squares that include England
tmp.medieval <- tmp.array[1:700, 4:5, 1:3]
#convert to a data frame and invert the rows/columns
tmp.medieval.df <- as.data.frame.array(tmp.medieval)
tmp.medieval.new <- t(tmp.medieval.df)
#set names to be the year values selected from the array
names <- seq(755, 1454, 1)
colnames(tmp.medieval.new) <- c(names)
#the data frame currently only has the temperature values. rows represent data from a single grid square. columns represent reconstructed climate by year. need to create a lat and lon column so that each row can be reaffiliated with the correct grid square.
#This can be done with repeating vectors (I checked the tables to ensure that the spatial data lined up with the correct temperature measurements)
lat.med <- as.vector(c(52.5, 57.5))
lon.med <- as.vector(c(-7.5, -2.5, 2.5))
tmp.medieval.new <- as.data.frame(tmp.medieval.new)
tmp.medieval.new <- data.frame(lon=rep(lon.med, len=6), tmp.medieval.new)
tmp.medieval.new <- data.frame(lat=rep(lat.med, len=6), tmp.medieval.new)
One important question to consider with regards to changing agricultural practice over time is the “medieval warming period.” The article for which the downloaded climate data was originally reconstructed argues for a general warming period from about AD900-1100. The abstract, however, mentions there is some regional variability. I want to know whether this impacted England and whether some regions had visibly different climatic conditions than others. Warmer temperatures could have made some land more accessible/desirable for wheat. In order to compare conditions in the Middle Saxon period to the Late Saxon period, I found the average temperature values for 50 year periods from 755 AD to 1104. I was concerned that setting a time period simply equivalent to the Middle and Late Saxon periods would gloss over any shifts in the climatic conditions. I thought a 50 year period would reduce my data to a much more reasonable size (later to be loaded into QGIS) without masking too much.
Ordinarily, I would try to construct a for loop for this, but I already spent a fair amount of time data wrangling to extract the data I needed and was able to save time with copy/paste.
tmp.medieval.new[,3:52] -> AD755
rowMeans(AD755) -> AD755
tmp.medieval.new[,53:102] -> AD804
rowMeans(AD804) -> AD804
tmp.medieval.new[,103:152] -> AD855
rowMeans(AD855) -> AD855
tmp.medieval.new[,153:202] -> AD904
rowMeans(AD904) -> AD904
tmp.medieval.new[,203:252] -> AD955
rowMeans(AD955) -> AD955
tmp.medieval.new[,253:302] -> AD1004
rowMeans(AD1004) -> AD1004
tmp.medieval.new[,303:352] -> AD1055
rowMeans(AD1055) -> AD1055
tmp.medieval.new[,353:402] -> AD1104
rowMeans(AD1104) -> AD1104
tmp.medieval.agg <- tmp.medieval.new[,1:2]
tmp.medieval.agg$AD755 <- AD755
tmp.medieval.agg$AD804 <- AD804
tmp.medieval.agg$AD855 <- AD855
tmp.medieval.agg$AD904 <- AD904
tmp.medieval.agg$AD955 <- AD955
tmp.medieval.agg$AD1004 <- AD1004
tmp.medieval.agg$AD1055 <- AD1055
tmp.medieval.agg$AD1104 <- AD1104
One of the themes of my project will be bringing data in and out of various formats to be used in the most suitable program. This code is commented out, but it shows how I exported the data into a csv for work in QGIS.
#csvfile <- "medenglandtemp.csv"
#write.table(na.omit(tmp.medieval.agg), csvfile, row.names = FALSE, sep = ",")
I played around with the climate data in QGIS and noticed a shift in all English grid squares to warmer conditions for the AD 905-954 and AD 955-1004 periods before returning to cooler conditions. I decided a line plot with each grid square representing a single line would help me to compare the changes in values over time.
tmp.med.plot <- t(tmp.medieval.agg)
tmp.med.plot <- tmp.med.plot[3:10,1:5]
y <- c(755, 804, 855, 904, 955, 1004, 1055, 1104)
tmp.med.plot <- as.data.frame(tmp.med.plot)
tmp.med.plot$y <- y
ggplot(tmp.med.plot, aes(tmp.med.plot$y)) + geom_line(aes(y=V1), col="red") + geom_line(aes(y=V2), col="yellow") + geom_line(aes(y=V3), col="green") + geom_line(aes(y=V4), col="blue")