基于Matplotlib在Python中绘制回归线并计算斜率的技术求助

Great to see you working with satellite data to study lake levels! I'll help you add regression lines to your plot and calculate the annual water level increase using linear regression. Here's how to modify your existing code to achieve this:

First, we need to convert datetime values to numeric ordinal values (since linear regression requires numerical inputs). Then we'll use numpy.polyfit to compute the slope and intercept for both your Total Beams and Hydrolabs datasets, generate the regression lines, and plot them alongside your original data. Finally, we'll calculate the annual increase from the daily slope values.

Modified Code with Regression Analysis #

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates

# Load and prepare data
heights = pd.read_excel('Qinghai_dates_heights.xlsx')

# Extract dates and convert to numeric ordinal values for regression
dates = pd.to_datetime(heights.loc[:, 'Date'])
dates_num = dates.apply(lambda x: x.toordinal())
total_beams = heights.loc[:, 'Total']

reference_dates = pd.to_datetime(heights.loc[:, 'Date.1'])
reference_dates_num = reference_dates.apply(lambda x: x.toordinal())
reference_heights = heights.loc[:, 'Hydrolabs']

# Perform linear regression for Total Beams dataset
slope_total, intercept_total = np.polyfit(dates_num, total_beams, 1)
regression_total = slope_total * dates_num + intercept_total

# Perform linear regression for Hydrolabs reference dataset
slope_ref, intercept_ref = np.polyfit(reference_dates_num, reference_heights, 1)
regression_ref = slope_ref * reference_dates_num + intercept_ref

# Calculate annual water level increase (slope is per day, multiply by 365 days)
annual_increase_total = slope_total * 365
annual_increase_ref = slope_ref * 365

# Print results to console
print(f"Total Beams Annual Water Level Increase: {annual_increase_total:.4f} meters")
print(f"Hydrolabs Reference Annual Water Level Increase: {annual_increase_ref:.4f} meters")

# Plotting setup
locator = mdates.MonthLocator()  # Show tick marks for every month
fmt = mdates.DateFormatter('%b')  # Format dates as Jan/Feb/Mar etc.

plt.figure(figsize=(10,6))
# Plot raw data points (no connecting lines for clarity)
plt.plot(dates, total_beams, label='Total Beams', marker="o", linestyle='None')
plt.plot(reference_dates, reference_heights, label='Reference height (Hydrolabs)', marker="o", linestyle='None')
# Plot regression lines
plt.plot(dates, regression_total, label=f'Total Beams Regression (slope: {slope_total:.6f} m/day)', linestyle='--', color='orange')
plt.plot(reference_dates, regression_ref, label=f'Hydrolabs Regression (slope: {slope_ref:.6f} m/day)', linestyle='--', color='green')

# Format x-axis
X = plt.gca().xaxis
X.set_major_locator(locator)
X.set_major_formatter(fmt)

# Add labels and title
plt.xlabel('Date [months]')
plt.ylabel('Elevation [m]')
plt.title("Water-Height Qinghai from November 2018 - November 2019 ")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Key Explanations: #

• Date Conversion: We use toordinal() to turn datetime objects into integers (days since January 1, 1 AD), which is compatible with linear regression algorithms.
• Linear Regression: np.polyfit(x, y, 1) computes the slope (first output) and intercept (second output) of the best-fit line for your data.
• Regression Lines: We generate line values using the basic linear equation y = slope * x + intercept.
• Annual Increase: The slope represents the daily change in water level. Multiplying by 365 gives the estimated annual change in meters.

Quick Notes: #

• If your Excel dates aren't automatically parsed as datetime objects, pd.to_datetime() ensures they're converted correctly.
• I adjusted the plot to show raw data points without connecting lines to make the regression lines clearer, but you can revert this by removing linestyle='None' if needed.

内容的提问来源于stack exchange，提问作者Wytse Petrie