Modeling global warming with a sinusoidal fit in Python

Question

I would like to model global warming from temperature records recorded daily from June 1920 to October 2019 in Montélimar on Python. To do this, I would first like to model these seasonal variations by a sinusoidal fit.

Here is the sinusoidal model that I have reproduced: $$T(t) = A \sin(\omega t + \phi) + B$$ where the parameters A (amplitude), 𝜙 (phase) and B (average temperature) are fitted to the data with $$\omega = \frac{2 \pi}{365}$$

However, such a model fitted to the whole data set does not give any increase in average temperature. I therefore try to apply a sinusoidal fit for each decade.

I first plotted the data in the data file like this:

# Import of modules
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit

# Read the montelimar_temperature containing the data
Date, Temperature = np.loadtxt('montelimar_temperature.dat', unpack = True)
Date = Date + 2400000.5 # Conversion of dates from Modified Julian Days to Julian Days
date_new = pd.to_datetime(Date, unit = 'D', origin = "julian") # Convert Julian Days dates to datetime64 format

# Size of the graph
plt.figure(figsize = (25, 9))

# Plot of the curve
plt.plot(date_new, Temperature)

# Title and others
plt.title("Evolution of daily recorded temperatures in Montélimar between 1920 and 2020", fontsize = 30) # Title
plt.xlabel("Year", fontsize = 20) # Title of the absissa
plt.ylabel("Temperature (°C)", fontsize = 20) # Title of the ordinates
plt.grid()
plt.show()

Then I created a time variable so I could do my decadal average. I applied the sine fit to all the decades in the data file, then plotted the entire graph with the fit. You will find my code below. It works flawlessly, but I feel like I'm repeating myself a lot with the code. I feel like I could do this without so much repetition, but every time I try I get errors and my graphs don't plot correctly. That's why I would like to know if there is a way to rewrite this in a cleaner way.

Here is the code I would like to simplify:

# Definition of the function for sinusoidal adjustment
def sinLaw(t, A, phi, B):
    omega = (2 * np.pi) / 365
    return A * np.sin(omega * t + phi) + B

import datetime

current_decade = np.datetime64(date_new[0], 'Y')
time_for_B = np.linspace(1930, 2020, 10)
#print(time_for_B)
count_time = np.array([]) # Variable to store the temperatures of a decade
count_date = np.array([])
B_list = np.array([])
n = 1

for i in range(0, len(date_new)): # Browse the date_new indexes to simply access the date_new value and the corresponding temperature
    if np.datetime64(date_new[i], 'Y') >= current_decade + np.timedelta64(10, 'Y'): # Arrived at a new decade
        current_decade = current_decade + np.timedelta64(10, 'Y') # Change in the current decade
        n = n + 1
        plt.figure(n)
        plt.plot(count_date, count_time)
        
        N = len(count_date)
        time_model = np.linspace(0, N, N)
        # Fit of the linear model
        solution = curve_fit(sinLaw, time_model, count_time)
        # Identification of the parameters
        A, phi, B = solution[0]
        # Display the result
        #print('A = {:4.2f} amplitude'.format(A))
        #print('B = {:4.2f} °C'.format(B))
        #print('phi = {:4.2f} radians'.format(phi))
        # Display the sine fit
        y = sinLaw(time_model, A, phi, B)
        B_list = np.append(B_list, B)
        plt.plot(count_date, y)
        
        errors = 5. * np.ones(y.shape)
        
        # Fit of the linear model
        solution, pcov = curve_fit(sinLaw, time_model, y, sigma = errors, absolute_sigma = True)
        
        # Identification of the model parameters
        A, phi, B = solution
        
        # Calculation of the uncertainty on the fitted parameters
        perr = np.sqrt(np.diag(pcov))
        
        # Display
        print('B = {:5.7f} ± {:5.3f} °C'.format(B, perr[0]))
        
        count_time = np.array([])
        count_date = np.array([])
        
    count_time = np.append(count_time, Temperature[i])
    count_date = np.append(count_date, date_new[i])

n = n + 1
plt.figure(n)
plt.plot(count_date, count_time, '.')

N = len(count_date)
time_model = np.linspace(0, N, N)
# Fit of the linear model
solution = curve_fit(sinLaw, time_model, count_time)
# Identification of the parameters
A, phi, B = solution[0]
# Display the result
#print('A = {:4.2f} amplitude'.format(A))
#print('B = {:4.2f} °C'.format(B))
#print('phi = {:4.2} radians'.format(phi))
# Display the sine fit
y = sinLaw(time_model, A, phi, B)
B_list = np.append(B_list, B)
plt.plot(count_date, y)
#print('B =', B_list)
plt.grid()
plt.figure(n + 1)
plt.plot(time_for_B, B_list)
plt.grid()

# Definition of the table of measurement errors
errors = 0.117 * np.ones(B_list.shape)

solution, pcov = curve_fit(sinLaw, time_model, count_time)

# Identification of the model parameters
A, phi, B = solution

perr = np.sqrt(np.diag(pcov))

# Display
print('B = {:5.7f} ± {:5.3f} °C'.format(B, perr[0]))

# Graphical representation of the data with the error bars
plt.errorbar(time_for_B, B_list, yerr = errors, marker = '+', linestyle = '')

# Graph option
plt.xlabel('Date [year]')
plt.ylabel('B [°C]')
plt.show()

Here is the result I get with this code and which corresponds to what I expect:

I can't put all the data of my file here because it contains 36296 lines and I can't attach a file but you will find below a sample of the first 100 lines of my file which has the name "montelimar_temperature.dat":

22553	22.6
22554	22.6
22555	24.8
22556	17.5
22557	16.5
22558	18.3
22559	18.1
22560	16.8
22561	17.0
22562	17.5
22563	18.0
22564	17.8
22565	17.7
22566	17.6
22567	16.9
22568	18.0
22569	18.4
22570	17.4
22571	19.8
22572	21.0
22573	22.5
22574	20.6
22575	19.6
22576	21.7
22577	22.6
22578	21.7
22579	21.4
22580	21.6
22581	19.6
22582	20.2
22583	18.0
22584	22.4
22585	22.7
22586	19.4
22587	13.6
22588	16.8
22589	16.0
22590	15.2
22591	15.0
22592	14.1
22593	17.7
22594	15.7
22595	16.4
22596	14.4
22597	19.8
22598	17.2
22599	19.4
22600	15.3
22601	14.3
22602	18.6
22603	18.8
22604	20.0
22605	16.6
22606	17.6
22607	16.0
22608	16.3
22609	16.2
22610	17.3
22611	15.4
22612	14.4
22613	13.0
22614	18.4
22615	16.0
22616	13.4
22617	11.9
22618	12.8
22619	11.6
22620	12.2
22621	11.0
22622	13.6
22623	14.0
22624	11.2
22625	10.3
22626	10.2
22627	6.8
22628	10.2
22629	10.0
22630	8.2
22631	8.0
22632	8.8
22633	8.0
22634	9.8
22635	8.0
22636	5.6
22637	7.8
22638	7.2
22639	5.6
22640	5.9
22641	6.9
22642	6.0
22643	6.6
22644	8.7
22645	11.5
22646	11.0
22647	8.9
22648	9.4
22649	8.9
22650	5.1
22651	1.5
22652	5.3

The columns have no names. The first column corresponds to dates in Modified Julian Days and the second column corresponds to temperatures in degrees Celsius.

Answer 1

The original format of the data file is neither a .dat nor a .csv (strictly), but a text file with a descriptive preamble followed by comma-separated data. This file, generated by

1. visiting https://www.ecad.eu/dailydata/customquery.php,
2. selecting non-blended, FRANCE, MONTELIMAR, "Mean temperature";
3. Next, download.
4. Unzip the file "ECA_non-blended_custom.zip"
5. Use "TG_SOUID106765.txt"

is first of all richer than the .dat that you have depicted as it contains this:

EUROPEAN CLIMATE ASSESSMENT & DATASET (ECA&D), file created on 14-03-2022
THESE DATA CAN BE USED FREELY PROVIDED THAT THE FOLLOWING SOURCE IS ACKNOWLEDGED:

Klein Tank, A.M.G. and Coauthors, 2002. Daily dataset of 20th-century surface
air temperature and precipitation series for the European Climate Assessment.
Int. J. of Climatol., 22, 1441-1453.
Data and metadata available at http://www.ecad.eu

FILE FORMAT (MISSING VALUE CODE IS -9999):

01-06 STAID: Station identifier
08-13 SOUID: Source identifier
15-22 DATE : Date YYYYMMDD
24-28 TG   : mean temperature in 0.1 °C
30-34 Q_TG : Quality code for TG (0='valid'; 1='suspect'; 9='missing')

This is the series (SOUID: 106765) of FRANCE, MONTELIMAR (STAID: 786)
See file sources.txt for more info.

and second of all includes five columns instead of just two. Pandas can process this trivially, via

def load(filename: str) -> tuple[
    pd.DataFrame,
    str,  # citation
]:
    df = pd.read_csv(
        filepath_or_buffer=filename,
        skiprows=21,
        skipinitialspace=True,
        parse_dates=['DATE'],
        infer_datetime_format=True,
    )
    valid_quality = 0
    df = df[df.Q_TG == valid_quality]

    with open(filename) as f:
        for _ in range(3):
            next(f)
        citation = ' '.join(line.rstrip() for line in itertools.islice(f, 3))

    return df, citation

This is far more preferable compared to the Julian-calendar processing that you have shown. As I mentioned in the comments, if you don't acknowledge that source in a citation you're violating the terms of use of the data set.

Add type hints.

Use n += 1 in-place addition rather than n = n + 1.

Are you sure that this:

np.linspace(0, N, N)

does what you think? From 0 to N, spaced one apart, there would be N+1 and not N points. So it's likely that you actually just want np.arange(N).

5. * np.ones(y.shape) should use np.full_like instead of multiplication. Same for 0.117 * np.ones(B_list.shape).

Use f-strings for your progress prints.

It is an anti-pattern to initialize a zero-sized Numpy array and iteratively grow it, as in count_time and count_date. These should be pre-allocated at the correct size.

Your code needs more functions. I've demonstrated how to add some but there need to be more.

Suggested

Example code covering some of the above:

import itertools

from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd


def load(filename: str) -> tuple[
    pd.DataFrame,
    str,  # citation
]:
    df = pd.read_csv(
        filepath_or_buffer=filename,
        skiprows=21,
        skipinitialspace=True,
        parse_dates=['DATE'],
        infer_datetime_format=True,
    )
    valid_quality = 0
    df = df[df.Q_TG == valid_quality]

    with open(filename) as f:
        for _ in range(3):
            next(f)
        citation = ' '.join(line.rstrip() for line in itertools.islice(f, 3))

    return df, citation


def plot_simple(df: pd.DataFrame, citation: str) -> plt.Figure:
    fig: plt.Figure
    ax: plt.Axes
    fig, ax = plt.subplots(figsize=(25, 9))
    ax.plot(df.DATE, df.TG)

    fig.suptitle("Evolution of daily recorded temperatures in Montélimar between 1920 and 2020")
    ax.set_title(citation)
    ax.set_xlabel("Year")
    ax.set_ylabel("Temperature (°C)")
    ax.grid()
    return fig


def plot_process(df: pd.DataFrame) -> None:
    date_new = df.DATE
    temperature = df.TG

    # Definition of the function for sinusoidal adjustment
    def sin_law(t: np.ndarray, A: float, phi: float, B: float) -> np.ndarray:
        omega = (2 * np.pi) / 365
        return A * np.sin(omega * t + phi) + B

    current_decade = np.datetime64(date_new.iloc[0], 'Y')
    time_for_B = np.linspace(1930, 2020, 10)

    count_time = np.array([])  # Variable to store the temperatures of a decade
    count_date = np.array([])
    B_list = np.array([])
    n = 1

    # Browse the date_new indexes to simply access the date_new value and the corresponding temperature
    for i in range(0, len(date_new)):
        if np.datetime64(date_new.iloc[i], 'Y') >= current_decade + np.timedelta64(10, 'Y'):  # Arrived at a new decade
            current_decade = current_decade + np.timedelta64(10, 'Y')  # Change in the current decade
            n += 1
            plt.figure(n)
            plt.plot(count_date, count_time)

            N = len(count_date)
            time_model = np.arange(N)

            solution = curve_fit(sin_law, time_model, count_time)

            # Identification of the parameters
            A, phi, B = solution[0]
            y = sin_law(time_model, A, phi, B)
            B_list = np.append(B_list, B)
            plt.plot(count_date, y)

            errors = np.full_like(y, 5)

            # Fit of the linear model
            solution, pcov = curve_fit(sin_law, time_model, y, sigma=errors, absolute_sigma=True)

            # Identification of the model parameters
            A, phi, B = solution

            # Calculation of the uncertainty on the fitted parameters
            perr = np.sqrt(np.diag(pcov))

            # Display
            print(f'B = {B:5.7f} ± {perr[0]:5.3f} °C')

            count_time = np.array([])
            count_date = np.array([])

        count_time = np.append(count_time, temperature.iloc[i])
        count_date = np.append(count_date, date_new.iloc[i])

    n += 1
    plt.figure(n)
    plt.plot(count_date, count_time, '.')

    N = len(count_date)
    time_model = np.arange(N)
    # Fit of the linear model
    solution = curve_fit(sin_law, time_model, count_time)
    # Identification of the parameters
    A, phi, B = solution[0]

    # Display the sine fit
    y = sin_law(time_model, A, phi, B)
    B_list = np.append(B_list, B)
    plt.plot(count_date, y)

    plt.grid()
    plt.figure(n + 1)
    plt.plot(time_for_B, B_list[:len(time_for_B)])
    plt.grid()

    # Definition of the table of measurement errors
    errors = np.full_like(B_list, 0.117)

    solution, pcov = curve_fit(sin_law, time_model, count_time)

    # Identification of the model parameters
    A, phi, B = solution

    perr = np.sqrt(np.diag(pcov))

    # Display
    print('B = {:5.7f} ± {:5.3f} °C'.format(B, perr[0]))

    # Graphical representation of the data with the error bars
    plt.errorbar(
        time_for_B,
        B_list[:len(time_for_B)],
        yerr=errors[:len(time_for_B)],
        marker='+', linestyle='',
    )

    # Graph option
    plt.xlabel('Date [year]')
    plt.ylabel('B [°C]')


def main() -> None:
    df, citation = load('TG_SOUID106765.txt')
    plot_simple(df, citation)
    plot_process(df)
    plt.show()


if __name__ == '__main__':
    main()