Here we present a digitalization effort of data from experiments cunducted by the Marine Ecology Laboratory at the Bedford Insstitute of Oceanography. As part of project activities we have acquired many such data sets and are currently working on digitizing them. The code given below is presented in Jupyter notebooks. Once this effort is completed a paper will be published along with the digitized data sets. If you are interested in such data sets please write to us!

Labrador sea 1978 photosynthesis irradiance data¶

Here we present an analysis of photosynthesis irradiance experiments digitized from a data report by B. Irwin and T. Platt. The pdf report is located in the same folder as this document (Data_Rep_83_Labrador_Sea_Oct_1977.pdf). A sample page from the report is given below:

All the photosynthesis irradiance experiments were digitized and saved in an xlsx format, also located in the same folder as this document (Labrador_sea_1978.xlsx).

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import scipy.io
import math
from scipy.optimize import curve_fit

# Loading the data set 
data = scipy.io.loadmat("PI.mat")

In [2]:

# Selecting the experiment by number
N = 35

# Plotting a single P-I curve
plt.figure(figsize=(10,5))
plt.scatter(data['PI'][N,14:24],data['PI'][N,24:34])
plt.xlabel(r'$I \left[\frac{W}{m^2}\right]$', fontsize = 15)
plt.ylabel(r'$p^B\ \left[\frac{ mgC}{mgChl\ h} \right]$', fontsize = 15)
plt.title('Photosynthesis irradiance measurements');

In [4]:

# Defining a photosynthesis irradiance function with initial slope, assimilation number and respiration as parameters
def PBIR(x,aB,PmB,RB): 
    return PmB*(1-np.exp(-aB*x/PmB)) - RB

\begin{equation} p^{B}(I) = P_{m}^{B}\left(1-\exp \left( -\alpha^B I\Big/P_m^B \right)\right) - R^B \end{equation}

In [6]:

# Iradiance [ W m^(-2) ]
I = np.linspace(0,500,500)

# Plotting an example photosynthesis irradiance function
plt.figure(figsize=(10,5))
plt.plot(I,PBIR(I,0.1,10,1), lw = 4, label = r'$p^B(I)$')
plt.title('An example photosynthesis irradiance function $p^B(I)$')
plt.xlabel(r'$I \left[\frac{W}{m^2}\right]$', fontsize = 15)
plt.ylabel(r'$p^B(I)\ \left[\frac{ mgC}{mgChl\ h} \right]$', fontsize = 15)
plt.axhline(0, c = 'lightgrey', lw = 1)
plt.axhline(-1, c = 'red', lw = 3, label = r'$R^B$')
plt.legend(loc=0, fontsize='small');
plt.ylim(-2,10);

In [10]:

# For simplicity declaring measured irradiance as x and measured production as y
x = data['PI'][N,14:24]
y = data['PI'][N,24:34]

# Fitting the photosynthesis irradiance function to the data
poi1, pci1 = curve_fit(PBIR,x,y,p0 = [0.25,5,1])

# Iradiance [ W m^(-2) ]
I = np.linspace(0,np.round(max(x)),100)

# Plotting the photosynthesis irradiance function
plt.figure(figsize=(10,5))
plt.scatter(x,y, c = 'black')
plt.plot(I,PBIR(I,poi1[0],poi1[1],poi1[2]), lw = 4, label = '5')
plt.title('Photosynthesis irradiance function $p^B(I)$')
plt.xlabel(r'$I \left[\frac{W}{m^2}\right]$', fontsize = 15)
plt.ylabel(r'$p^B(I)\ \left[\frac{ mgC}{mgChl\ h} \right]$', fontsize = 15)
plt.ylim(-1.1*poi1[2],1.1*poi1[1]);

In [14]:

# Iradiance [ W m^(-2) ]
I = np.linspace(0.1,1000,10000)

# Plotting the photosynthesis irradiance function
plt.figure(figsize=(10,5))
plt.scatter(x,y, c = 'black')
plt.plot(I,PBIR(I,poi1[0],poi1[1],poi1[2]), lw = 4, label = '5')
plt.title('Photosynthesis irradiance function $p^B(I)$')
plt.xlabel(r'$I \left[\frac{W}{m^2}\right]$', fontsize = 15)
plt.ylabel(r'$p^B(I)\ \left[\frac{ mgC}{mgChl\ h} \right]$', fontsize = 15)
plt.xscale('log')
plt.xlim(0.1,1000);
plt.ylim(-1.1*poi1[2],1.1*poi1[1]);

In [15]:

# Selecting the experiment by number
PmBN = np.zeros(len(data['PI']))
aBN = np.zeros(len(data['PI']))
RBN = np.zeros(len(data['PI']))

IN = np.zeros(len(data['PI'])*10)
PBN = np.zeros(len(data['PI'])*10)
eN = np.zeros(len(data['PI'])*10)
mes = np.zeros(len(data['PI'])*10)
mod = np.zeros(len(data['PI'])*10)

for n in range(len(data['PI'])):

    # For simplicity declaring measured irradiance as x and measured production as y
    x = data['PI'][n,14:24]
    y = data['PI'][n,24:34]

    # Fitting the photosynthesis irradiance function to the data
    poi1, pci1 = curve_fit(PBIR,x,y,p0 = [0.25,5,1])
    
    aBN[n] = poi1[0]
    PmBN[n] = poi1[1]
    RBN[n] = poi1[2]
    
    for m in range(10):
        IN[n*10 + m] = aBN[n]*x[m]/PmBN[n]
        PBN[n*10 + m] = (y[m]+RBN[n])/PmBN[n]
        eN[n*10 + m] = (y[m] - PBIR(x[m],aBN[n],PmBN[m],RBN[n]))
        mes[n*10 + m] = y[m]
        mod[n*10 + m] = PBIR(x[m],poi1[0],poi1[1],poi1[2])
    

C:\Users\Kovac\AppData\Local\Temp\ipykernel_16344\380507425.py:3: RuntimeWarning: divide by zero encountered in scalar divide
  return PmB*(1-np.exp(-aB*x/PmB)) - RB

In [19]:

# Iradiance [ W m^(-2) ]
I = np.linspace(0,np.round(max(IN)),100)

# Plotting the photosynthesis irradiance function
plt.figure(figsize=(10,5))
plt.plot(I,PBIR(I,1,1,0), lw = 4)
plt.scatter(IN,PBN, c = 'black', s = 12)
plt.title('Photosynthesis irradiance function $p^B(I)$')
plt.xlabel(r'$I/I_k$', fontsize = 15)
plt.ylabel(r'$P^B/P_m^B$', fontsize = 15)
plt.xlim(0,np.round(max(IN)));
plt.ylim(0,1.2);

In [23]:

# Iradiance [ W m^(-2) ]
I = np.linspace(0.01,1000,100000)

# Plotting the photosynthesis irradiance function
plt.figure(figsize=(10,5))
plt.plot(I,PBIR(I,1,1,0), lw = 4)
plt.scatter(IN,PBN, c = 'black', s = 12)
plt.title('Photosynthesis irradiance function $p^B(I)$')
plt.xscale('log')
plt.xlabel(r'$I/I_k$', fontsize = 15)
plt.ylabel(r'$P^B/P_m^B$', fontsize = 15)
plt.xlim(0.01,100);
plt.ylim(0,1.3);

In [24]:

# Plotting the histograms of the photosynthesis parameters
fig, axs = plt.subplots(2, 2, figsize=(10, 7)) 
plt.subplots_adjust(hspace=0.4)

axs[0,0].hist(aBN, bins=15, edgecolor='black', color = 'orange')
axs[0,0].axvline(np.mean(aBN), lw = 4, c = 'black', label = 'Selected initial slope')
axs[0,0].set_xlabel(r'$\alpha^B \left[\frac{mgC}{mgChl\ W m^{-2} h}  \right]$')
axs[0,0].set_ylabel('Frequency')
axs[0,0].set_title('Histogram of the initial slope');

axs[0,1].hist(PmBN, bins=15, edgecolor='black', color = 'purple')
axs[0,1].axvline(np.mean(PmBN), lw = 4, c = 'black', label = 'Selected assimilation number')
axs[0,1].set_xlabel(r'$P_m^B \left[\frac{ mgC}{mgChl\ h}  \right]$')
axs[0,1].set_ylabel('Frequency')
axs[0,1].set_title('Histogram of the assimilation number');

axs[1,0].hist(PmBN/aBN, bins=15, edgecolor='black')
axs[1,0].axvline(np.mean(PmBN/aBN), lw = 4, c = 'black', label = 'Selected initial slope')
axs[1,0].set_xlabel(r'$I_k \left[\frac{W}{m^{-2}} \right]$')
axs[1,0].set_ylabel('Frequency')
axs[1,0].set_title('Histogram of the photoadaptation parameter');

axs[1,1].hist(RBN[RBN>0], bins=15, edgecolor='black', color = 'green')
axs[1,1].axvline(np.mean(RBN), lw = 4, c = 'black', label = 'Selected initial slope')
axs[1,1].set_xlabel(r'$P_m^B \left[\frac{ mgC}{mgChl\ h}  \right]$')
axs[1,1].set_ylabel('Frequency')
axs[1,1].set_title('Histogram of the respiration parameter');

In [25]:

# Plotting the errors histogram 
plt.hist(eN, bins = 25, edgecolor = 'black')
plt.xlabel('Relative error')
plt.ylabel('Frequency');

In [30]:

# Model versus measurements on a normaln scale
plt.figure(figsize=(8,6))
plt.scatter(mes,mod, s = 10)
plt.plot(np.linspace(0.01,3,100),np.linspace(0.01,3,100), c = 'black')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Measured ' + r'$P^B \left[\frac{ mgC}{mgChl\ h} \right]$')
plt.ylabel('Model ' + r'$P^B \left[\frac{ mgC}{mgChl\ h} \right]$');

In [31]:

# Plotting the relative error as a function of irradiance
plt.figure(figsize=(10,5))
plt.scatter(IN,eN, c = 'black', s = 12)
plt.title('Relative error')
plt.xscale('log')
plt.xlabel(r'$I/I_k$', fontsize = 15)
plt.ylabel('Relative error', fontsize = 15)
plt.xlim(0.01,100);

Workshops

Digitalization of archived data

Labrador sea 1978 photosynthesis irradiance data¶

Measuring primary production

Modelling primary production