4
\$\begingroup\$

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: output

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.

\$\endgroup\$
4
  • 1
    \$\begingroup\$ Making some guesses, downloading Montelimar mean temperature with no blending, the return format is NOT the format that you've shown, which means that you're likely pre-processing and have not described this. The returned format is a five-column CSV with full headers and a significant preamble. So... what are you actually doing? \$\endgroup\$
    – Reinderien
    Commented Mar 14, 2022 at 15:21
  • 1
    \$\begingroup\$ p.s. the downloaded data from your indicated website do not use Julian time either. So this question is incomplete. \$\endgroup\$
    – Reinderien
    Commented Mar 14, 2022 at 15:48
  • 1
    \$\begingroup\$ p.p.s. You're violating the terms of use of the dataset, which require a citation to the Int. J. of Climatol. \$\endgroup\$
    – Reinderien
    Commented Mar 14, 2022 at 16:03
  • \$\begingroup\$ Thanks for the feedback @Reinderien. I just added a sample of the data file. I did provide my entire code. The original format of the data file is indeed a .dat format and not a .csv. If the link does not work anymore I will remove it. I apologize for this inconvenience \$\endgroup\$
    – Donny
    Commented Mar 14, 2022 at 18:17

1 Answer 1

2
\$\begingroup\$

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()
\$\endgroup\$
2
  • \$\begingroup\$ Thanks for your advice and suggestions. np.linspace(0, N, N) does what I want. It allows me to have N elements precisely. Regarding the code example you provided, the load function gives the following error: TypeError: 'type' object is not subscriptable because the first column is a date object. \$\endgroup\$
    – Donny
    Commented Mar 17, 2022 at 14:15
  • \$\begingroup\$ Update your version of Python. Your diagnosis of the problem is not correct. Your version of Python does not support type hints. \$\endgroup\$
    – Reinderien
    Commented Mar 17, 2022 at 14:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.