Tracking the Thaw: Modeling the Antarctic Ice Cap with Python

The Antarctic ice cap is one of the most critical barometers for our planet’s climate health. But when we look at 45 years of satellite data, how do we differentiate between yearly "noise" and a definitive trend? Today, we’re using Linear and Quadratic Regression to analyze ice cap area data from 1979 to 2024 and project what the horizon looks like for 2030.

​Starting in 1979, the September ice cap area sat at roughly 7.2 Million sq km. Fast forward to the 2020s, and we are frequently seeing measurements dip below 4.5 Million sq km.

​To make sense of this, we’ve written a script to calculate the "Goodness of Fit" (R^2) and project future values.

import numpy as np
import matplotlib.pyplot as plt
plt.cla()
def r_squared(actual,predicted):
    ss_res=np.sum((actual-predicted)**2)
    ss_tot=np.sum((actual-np.mean(actual))**2)
    return 1-(ss_res/ss_tot)

def curve_fit(x,y,deg=1):
    return np.polyfit(x,y, deg)
    
#antarctic ice cap in septembe4 
ice_cap_area=[
    #Survey Yr,Ice area(M sq km) 
   (1979, 7.2),
   (1980, 7.9), 
   (1981, 7.3), 
   (1982, 7.5), 
   (1983, 7.5),
   (1984, 7.2), 
   (1985, 6.9), 
   (1986, 7.5), 
   (1987, 7.5), 
   (1988, 7.5)
,
  (1989, 7.0),
   (1990,6.2), 
   (1991,6.6), 
   (1992, 7.6), 
  (1993, 6.5), 
  (1994,7.2 ),
  (1995, 6.1), 
  (1996, 7.9), 
  (1997,6.7),
  (1998, 6.6), 
  (1999, 6.2), 
  (2000, 6.3),
  (2001, 6.8), 
  (2002, 5.9),
   (2003, 6.2),
(2004, 6.1), 
(2005 ,5.6), 
  (2006, 5.9),
  (2007,4.3 ),
  (2008, 4.7), 
   (2009, 5.4),
(2010, 4.9), 
  (2011, 4.6),
(2012, 3.6), #record 
(2013, 5.2),
(2014, 5.3),
(2015, 4.6),
(2016, 4.4),
(2017, 4.8),
(2018, 4.7),
(2019, 4.3),
(2020, 4.2),
(2021, 4.9),
(2022, 4.9),
(2023, 4.4),
(2024, 4.4),
  ]
predict_year=2030
#take 1979 as 0
time_1979=np.array([x-1979 for (x,y) in ice_cap_area])
#print(time_1979)

cap_area=np.array([y for (x,y) in ice_cap_area])
#print(cap_area)

#try linear fit
print("Try a linear fit between data points")
deg_lin=1
lin_fit_coeffs=curve_fit(time_1979,cap_area,deg_lin)
#print(lin_fit_coeffs)
slope=lin_fit_coeffs[0]
intercept=lin_fit_coeffs[1]
print(f"Eqn: ice_cap_area={slope:.4f}*time+{intercept:.4f}")
pred_cap_area_lin=slope*time_1979+intercept
r2_fit_lin=r_squared(cap_area,pred_cap_area_lin)
print(f"R2:{r2_fit_lin:.5f}")
ice_area_2030_lin=slope*(predict_year-1979)+intercept
print(f"Predicted Ice cap area in {predict_year}:{ice_area_2030_lin:.3f} Msqkm")



print()
#nasa asks for quadratic regression also
print("Try a quadratic fit between data")
deg_quad=2
quad_fit_coeffs=curve_fit(time_1979,cap_area,deg_quad)
#print(quad_fit_coeffs)
coeff_x2=quad_fit_coeffs[0]
coeff_x=quad_fit_coeffs[1]
coeff_k=quad_fit_coeffs[2]
print(f"ice_cap_area={coeff_x2:.4f}t2{coeff_x2:.4f}t+{coeff_k:.4f}")
pred_cap_area_quad=coeff_x2*time_1979**2+coeff_x*time_1979+coeff_k
r2_fit_quad=r_squared(cap_area,pred_cap_area_quad)
print(f"R2:{r2_fit_quad:.5f}")
ice_area_2030_quad=coeff_x2*(predict_year-1979)**2+coeff_x*(predict_year-1979)+coeff_k

print(f"Predicted Ice cap area in {predict_year}:{ice_area_2030_quad:.3f} Msqmt")



# Create a smooth range of x-values for the curves
x_range = np.linspace(0, time_1979[-1], 100)

# Calculate y-values for both models
y_lin = slope * x_range + intercept
y_quad = coeff_x2 * x_range**2 + coeff_x * x_range + coeff_k

# Plotting
plt.scatter(time_1979, cap_area, color='gray', alpha=0.5, label='Actual Data')
plt.plot(x_range, y_lin, color='blue', label=f'Linear (R²={r2_fit_lin:.3f})')
plt.plot(x_range, y_quad, color='red', linestyle='--', label=f'Quadratic (R²={r2_fit_quad:.3f})')

plt.xlabel('Years since 1979')
plt.ylabel('Ice Cap Area (M sq km)')
plt.title('Antarctic Ice Cap: Linear vs Quadratic Trend')
plt.legend()
plt.grid(True)
plt.show()