setwd("~/Desktop/inPrep/pew/github/erp/Rmd")

Temperature extraction

This R markdown document has been developed to extract HADISST (Global sea ice and Sea Surface Temperature analyses) information for specific shapefiles.

This specific aim of this document is to extract SST information for the area covering the distribution of North East Atlantic Mackerel.

Two areas are examined across two time periods. The first area is the whole stock distribution (subareas 1-8 and 14 and Division 9.a) and the second is an area covering the key spawning areas (West of Ireland, Bay of Biscay, and Northern North Sea). Temperature time series for these areas are extracted as 1) an annual average and 2) a March-May average as this is the key spawning time for Mackerel.

Functions for extracting temperature data

The code below documents the functions developed for the code pipeline.

# Used to generate annual average temps by getting the mean of 12 columns in a sequence
byapply <- function(x, by, fun, ...)
{
  # Create index list
  if (length(by) == 1)
  {
    nc <- ncol(x)
    split.index <- rep(1:ceiling(nc / by), each = by, length.out = nc)
  } else # 'by' is a vector of groups
  {
    nc <- length(by)
    split.index <- by
  }
  index.list <- split(seq(from = 1, to = nc), split.index)
  
  # Pass index list to fun using sapply() and return object
  sapply(index.list, function(i)
  {
    do.call(fun, list(x[, i], ...))
  })
}

#function to change number of x axis ticks
number_ticks <- function(n) {function(limits) pretty(limits, n)} 

#function to plot a moving average
ma <- function(x, n = 5){filter(x, rep(1 / n, n), sides = 2)}

Read in raster file for temperature data

# Read in nc temperature file #####
#HADISST nc file from https://www.metoffice.gov.uk/hadobs/hadisst/data/download.html

hann<-brick("../data/inputs/ecosystem/HadISST_sst.nc") 

# Raster to point conversion (create spatial and numeric dataframes)

hpts<-rasterToPoints(hann, fun=NULL, spatial=TRUE)
hptsvals<-rasterToPoints(hann, fun=NULL, spatial=FALSE)

Outputs for the whole stock area

This section provides an image of the overlap between the whole stock area and temperature reference points from HADISST, along with plots of the average changes in temperature for the area.

# Read in  shapefile

shapefile_path<-file.path("../data/inputs/ecosystem", "Mackerel_area.shp")

range1<-readOGR(dsn=shapefile_path,layer="Mackerel_area")  # dsn = directory layer= shapefile without extension

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/laurie/Desktop/inPrep/pew/github/data/inputs/ecosystem/Mackerel_area.shp", layer: "Mackerel_area"
## with 1 features
## It has 10 fields
## Integer64 fields read as strings:  OBJECTID_1 OBJECTID

range1 <- spTransform(range1, CRS(proj4string(hann))) # match coordinate reference systems

# Check CRS of raster
print(proj4string(hann))

## [1] "+proj=longlat +datum=WGS84 +no_defs"

# Check CRS of shapefile
print(proj4string(range1))

## [1] "+proj=longlat +datum=WGS84 +no_defs"

# Extract annual temperature values from HADISST overlapping area

hrangepts<-!is.na(over(hpts,range1))
hptsv<-data.frame(hptsvals)
hrangevals<-hptsv[hrangepts==TRUE,]
# Replace -1000 with NA in the specified column
hrangevals<- replace(hrangevals, hrangevals == -1000, NA)
hrangevals<-subset(hrangevals,!is.na(x))

plot(range1)  #this is plotting the shapefile area
points(hrangevals$x,hrangevals$y)   #this is plotting the temperature points extracted in the area  overlapping in the plotspace

Figure 1. Overlay between mackerel stock area shapefile and temperature extraction points from the HADISST dataset.

#Create dataframe for temperature time series at each overlapping point

hrangevals_mar_may <- hrangevals[, grepl("\\b.03.|\\b.04.|\\b.05.\\b", names(hrangevals))]

# March:May area average SST

thrange_mar_may<-t(hrangevals_mar_may[,1:ncol(hrangevals_mar_may)]) # Need to transpose data to give year columns

# Calculate the row means for every three rows
num_rows <- nrow(thrange_mar_may)
group_size <- 3
num_groups <- num_rows %/% group_size

tssst_mar_may <- numeric(num_groups)  # Initialize a vector to store the results

for (i in 1:num_groups) {
  start_row <- (i - 1) * group_size + 1
  end_row <- start_row + group_size - 1
  tssst_mar_may[i] <- mean(thrange_mar_may[start_row:end_row, ], na.rm = TRUE)
}

years_mar_may<-1870:2020 # change if data span changes
tempdf_mar_may<-as.data.frame(cbind(years_mar_may,tssst_mar_may)) #create single df

# Create dataframe for temperature time series at each overlapping point

hrangevals_annual <- byapply(hrangevals[3:1812], 12, function(x) rowMeans(x, na.rm = TRUE))              # average every 12 columns to make annual averages

# Annual area average SST

thrange_annual<-t(hrangevals_annual[,1:ncol(hrangevals_annual)]) # Need to transpose data to give year columns
tssst_annual<-rowMeans(thrange_annual,na.rm=TRUE) # average temp for each pixel row  (time series of sea surface temperature : tssst)

temp_df<-cbind(tssst_annual,tssst_mar_may)


years_annual<-1870:2020 # change if data span changes
tempdf_annual<-as.data.frame(cbind(years_annual,tssst_annual)) #create single df
hrangevals<-hptsv[hrangepts==TRUE,]
hrangevals<-subset(hrangevals,!is.na(x))

mid<-mean(tssst_annual) # mean value used for visuals

ma_data<-as.numeric(tempdf_annual$tssst_annual)

# Plot code
temp_plot_annual <- ggplot(tempdf_annual, aes(x=years_annual, y=tssst_annual, fill=tssst_annual)) +
  geom_line(color='black') +
  geom_line(aes(y=ma(ma_data,n=10)), lwd=1.1, color='darkorange', linetype=1) +
  geom_point(size=2, shape=21) +
  scale_fill_gradient2(midpoint=mid, low="blue", mid="white", high="red", space="Lab", guide="none") +  # Remove the legend
  theme_minimal() +
  theme(panel.background = element_rect(fill = "white", color = "black")) +
  scale_x_continuous(breaks = number_ticks(7), name='Year') +
  scale_y_continuous(name="Temperature (\u00B0C)") +
  ggtitle("Mackerel spawning area SST (HADISST, annual)")

temp_plot_annual #view plot

Figure 2. Annual sea surface temperature averaged over the stock area of Atlantic Mackerel.

ggsave("../outputs/figs/Temp_plot_Stock1.png", temp_plot_annual, height = 4, width = 5,dpi = 500)

mid<-mean(tssst_mar_may) # mean value used for visuals

ma_data<-as.numeric(tempdf_mar_may$tssst_mar_may)

# Plot code
temp_plot_may_mar <- ggplot(tempdf_mar_may, aes(x=years_mar_may, y=tssst_mar_may, fill=tssst_mar_may)) +
  geom_line(color='black') +
  geom_line(aes(y=ma(ma_data,n=10)), lwd=1.1, color='darkorange', linetype=1) +
  geom_point(size=2, shape=21) +
  scale_fill_gradient2(midpoint=mid, low="blue", mid="white", high="red", space="Lab", guide="none") +  # Remove the legend
  theme_minimal() +
  theme(panel.background = element_rect(fill = "white", color = "black")) +
  scale_x_continuous(breaks = number_ticks(7), name='Year') +
  scale_y_continuous(name="Temperature (\u00B0C)") +
  ggtitle("Mackerel stock area SST (HADISST, March-May)")

temp_plot_may_mar #view plot

Figure 3. March-May sea surface temperature averaged over the stock area of Atlantic Mackerel.

ggsave("../outputs/figs/Temp_plot_Stock2.png", temp_plot_may_mar, height = 4, width = 5,dpi = 500)

Outputs for the stocks spawning area

This section provides an image of the overlap between the stocks spawning area and temperature reference points from HADISST, along with plots of the average changes in temperature for the area.

# Read in  shapefile

shapefile_path<-file.path("../data/inputs/ecosystem", "Mackerel_s.shp")

range1<-readOGR(dsn=shapefile_path,layer="Mackerel_s")  # dsn = directory layer= shapefile without extension

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/laurie/Desktop/inPrep/pew/github/data/inputs/ecosystem/Mackerel_s.shp", layer: "Mackerel_s"
## with 1 features
## It has 10 fields
## Integer64 fields read as strings:  OBJECTID_1 OBJECTID

range1 <- spTransform(range1, CRS(proj4string(hann))) # match coordinate reference systems

# Check CRS of raster
print(proj4string(hann))

## [1] "+proj=longlat +datum=WGS84 +no_defs"

# Check CRS of shapefile
print(proj4string(range1))

## [1] "+proj=longlat +datum=WGS84 +no_defs"

# Extract annual temperature values from HADISST overlapping area

hrangepts<-!is.na(over(hpts,range1))
hptsv<-data.frame(hptsvals)
hrangevals<-hptsv[hrangepts==TRUE,]
# Replace -1000 with NA in the specified column
hrangevals<- replace(hrangevals, hrangevals == -1000, NA)
hrangevals<-subset(hrangevals,!is.na(x))

plot(range1)  #this is plotting the shapefile area
points(hrangevals$x,hrangevals$y)   #this is plotting the temperature points extracted in the area  overlapping in the plotspace

Figure 4. Overlay between mackerel spawning area shapefile and temperature extraction points from the HADISST dataset.

#Create dataframe for temperature time series at each overlapping point

hrangevals_mar_may <- hrangevals[, grepl("\\b.03.|\\b.04.|\\b.05.\\b", names(hrangevals))]

# March:May area average SST

thrange_mar_may<-t(hrangevals_mar_may[,1:ncol(hrangevals_mar_may)]) # Need to transpose data to give year columns

# Calculate the row means for every three rows
num_rows <- nrow(thrange_mar_may)
group_size <- 3
num_groups <- num_rows %/% group_size

tssst_mar_may <- numeric(num_groups)  # Initialize a vector to store the results

for (i in 1:num_groups) {
  start_row <- (i - 1) * group_size + 1
  end_row <- start_row + group_size - 1
  tssst_mar_may[i] <- mean(thrange_mar_may[start_row:end_row, ], na.rm = TRUE)
}

years_mar_may<-1870:2020 # change if data span changes
tempdf_mar_may<-as.data.frame(cbind(years_mar_may,tssst_mar_may)) #create single df

# Create dataframe for temperature time series at each overlapping point

hrangevals_annual <- byapply(hrangevals[3:1812], 12, function(x) rowMeans(x, na.rm = TRUE))              # average every 12 columns to make annual averages

# Annual area average SST

thrange_annual<-t(hrangevals_annual[,1:ncol(hrangevals_annual)]) # Need to transpose data to give year columns
tssst_annual<-rowMeans(thrange_annual,na.rm=TRUE) # average temp for each pixel row  (time series of sea surface temperature : tssst)

temp_df<-cbind(tssst_annual,tssst_mar_may)


years_annual<-1870:2020 # change if data span changes
tempdf_annual<-as.data.frame(cbind(years_annual,tssst_annual)) #create single df
hrangevals<-hptsv[hrangepts==TRUE,]
hrangevals<-subset(hrangevals,!is.na(x))

mid<-mean(tssst_annual) # mean value used for visuals

ma_data<-as.numeric(tempdf_annual$tssst_annual)

# Plot code
temp_plot_annual <- ggplot(tempdf_annual, aes(x=years_annual, y=tssst_annual, fill=tssst_annual)) +
  geom_line(color='black') +
  geom_line(aes(y=ma(ma_data,n=10)), lwd=1.1, color='darkorange', linetype=1) +
  geom_point(size=2, shape=21) +
  scale_fill_gradient2(midpoint=mid, low="blue", mid="white", high="red", space="Lab", guide="none") +  # Remove the legend
  theme_minimal() +
  theme(panel.background = element_rect(fill = "white", color = "black")) +
  scale_x_continuous(breaks = number_ticks(7), name='Year') +
  scale_y_continuous(name="Temperature (\u00B0C)") +
  ggtitle("Mackerel spawning area SST (HADISST, annual)")

temp_plot_annual #view plot

Figure 5. Annual sea surface temperature averaged over the spawning area of Atlantic Mackerel.

ggsave("../outputs/figs/Temp_plot_spawning1.png", temp_plot_annual, height = 4, width = 5,dpi = 500)

mid<-mean(tssst_mar_may) # mean value used for visuals

ma_data<-as.numeric(tempdf_mar_may$tssst_mar_may)

# Plot code
temp_plot_may_mar<- 
  ggplot(tempdf_mar_may, aes(x=years_mar_may, y=tssst_mar_may, fill=tssst_mar_may)) +
  geom_line(color='black') +
  geom_line(aes(y=ma(ma_data,n=10)), lwd=1.1, color='darkorange', linetype=1) +
  geom_point(size=2, shape=21) +
  scale_fill_gradient2(midpoint=mid, low="blue", mid="white", high="red", space="Lab", guide="none") +  # Remove the legend
  theme_minimal() +
  theme(panel.background = element_rect(fill = "white", color = "black")) +
  scale_x_continuous(breaks = number_ticks(7), name='Year') +
  scale_y_continuous(name="Temperature (\u00B0C)") +
  ggtitle("Mackerel spawning area SST (HADISST, March-May)")

temp_plot_may_mar #view plot

Figure 6. March-May sea surface temperature averaged over the spawning area of Atlantic Mackerel.

ggsave("../outputs/figs/Temp_plot_Spawning2.png", temp_plot_may_mar, height = 4, width = 5,dpi = 500)

Pew Ecological Reference Points: HADISST temperature extraction

Dr Jacob Bentley

Last compiled on November 21, 2023

Temperature extraction

Functions for extracting temperature data

Read in raster file for temperature data

Outputs for the whole stock area

Outputs for the stocks spawning area