Linear Regression across time for Lat Lon Map¶

Author: Abdullah A. Fahad (a.fahad@nasa.gov)

import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt

!pip install cartopy
import cartopy.crs as ccrs

!wget https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2.highres/sst.mon.mean.nc

# Read the data using xarray
ds = xr.open_dataset('/content/sst.mon.mean.nc').sel(time=slice('1982','2022'))

#ds = xr.decode_cf(ds, use_cftime=True) # Ignore this line. It is only needed if time coordinates are given in floating data and you need to convert it to cftime format
ds

<xarray.Dataset>
Dimensions:  (time: 492, lat: 720, lon: 1440)
Coordinates:
  * time     (time) datetime64[ns] 1982-01-01 1982-02-01 ... 2022-12-01
  * lat      (lat) float32 -89.88 -89.62 -89.38 -89.12 ... 89.38 89.62 89.88
  * lon      (lon) float32 0.125 0.375 0.625 0.875 ... 359.1 359.4 359.6 359.9
Data variables:
    sst      (time, lat, lon) float32 ...
Attributes:
    Conventions:    CF-1.5
    title:          NOAA/NCEI 1/4 Degree Daily Optimum Interpolation Sea Surf...
    institution:    NOAA/National Centers for Environmental Information
    source:         NOAA/NCEI https://www.ncei.noaa.gov/data/sea-surface-temp...
    References:     https://www.psl.noaa.gov/data/gridded/data.noaa.oisst.v2....
    dataset_title:  NOAA Daily Optimum Interpolation Sea Surface Temperature
    version:        Version 2.1
    comment:        Reynolds, et al.(2007) Daily High-Resolution-Blended Anal...

array(['1982-01-01T00:00:00.000000000', '1982-02-01T00:00:00.000000000',
       '1982-03-01T00:00:00.000000000', ..., '2022-10-01T00:00:00.000000000',
       '2022-11-01T00:00:00.000000000', '2022-12-01T00:00:00.000000000'],
      dtype='datetime64[ns]')

array([-89.875, -89.625, -89.375, ...,  89.375,  89.625,  89.875],
      dtype=float32)

array([1.25000e-01, 3.75000e-01, 6.25000e-01, ..., 3.59375e+02, 3.59625e+02,
       3.59875e+02], dtype=float32)

[510105600 values with dtype=float32]

PandasIndex(DatetimeIndex(['1982-01-01', '1982-02-01', '1982-03-01', '1982-04-01',
               '1982-05-01', '1982-06-01', '1982-07-01', '1982-08-01',
               '1982-09-01', '1982-10-01',
               ...
               '2022-03-01', '2022-04-01', '2022-05-01', '2022-06-01',
               '2022-07-01', '2022-08-01', '2022-09-01', '2022-10-01',
               '2022-11-01', '2022-12-01'],
              dtype='datetime64[ns]', name='time', length=492, freq=None))

PandasIndex(Float64Index([-89.875, -89.625, -89.375, -89.125, -88.875, -88.625, -88.375,
              -88.125, -87.875, -87.625,
              ...
               87.625,  87.875,  88.125,  88.375,  88.625,  88.875,  89.125,
               89.375,  89.625,  89.875],
             dtype='float64', name='lat', length=720))

PandasIndex(Float64Index([  0.125,   0.375,   0.625,   0.875,   1.125,   1.375,   1.625,
                1.875,   2.125,   2.375,
              ...
              357.625, 357.875, 358.125, 358.375, 358.625, 358.875, 359.125,
              359.375, 359.625, 359.875],
             dtype='float64', name='lon', length=1440))

# converting the data to numpy array
lon=ds.lon.data
lat=ds.lat.data
time=ds.time.data
nlon=len(lon); nlat=len(lat); nt=len(time)
sst=ds.sst.data
sst.shape

(492, 720, 1440)

nmonth=12
nyear=nt//nmonth
year=np.arange(1982,2022+1)

#reshaping data to year, mon, lat, lon
sst=np.reshape(sst,(nyear,nmonth,nlat,nlon))
sst.shape

(41, 12, 720, 1440)

#making annual mean data by taking average ver the month dimension
sst_annual=np.nanmean(sst,1)
sst_annual.shape

<ipython-input-25-7dbb5e229cc7>:2: RuntimeWarning: Mean of empty slice
  sst_annual=np.nanmean(sst,1)

(41, 720, 1440)

#plotting Climatology of Annual mean SST

plt.figure(figsize=(10, 6))
plt.contourf(lon, lat, np.nanmean(sst_annual,0),cmap='YlOrRd')
plt.colorbar(label='SST')
plt.title('Sea Surface Temperature')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.show()

<ipython-input-27-a48716cadc65>:4: RuntimeWarning: Mean of empty slice
  plt.contourf(lon, lat, np.nanmean(sst_annual,0),cmap='YlOrRd')

## function for linear regression
from scipy import stats
from numpy import *
def l_trend(var,lon,lat,time,sig=False):
    nlon=len(lon)
    nlat=len(lat)
    nt=len(time)
    vart=zeros(nlat*nlon)
    varp=zeros(nlat*nlon)

    if len(var.shape)== 3:
        var=reshape(var,(nt,nlat*nlon))
        print('l_trend: assuming variable as 3D [time,lat,lon]')
        for i in range(nlat*nlon):
            v=var[:,i]
            vart[i], intercept, r_value, varp[i], std_err=stats.linregress(time,v)

        vart=reshape(vart,(nlat,nlon))
        varp=reshape(varp,(nlat,nlon))
        #return (vart,varp)

    else:
        raise ValueError('Variable shape is not 2D or 3D. plese instert variable in this format var[time,lat,lon] or var[time,lon*lat]')

    if sig==False:
        return (vart, varp)
    else:
        for i in range(nlat):
            for j in range (nlon):
                if varp[i,j]>sig:
                  vart[i,j]=nan
        return (vart, varp)

## calculating SST linear trend and their pvalues

SST_trend, SST_pvalue=l_trend(sst_annual, lon, lat, time=year)

l_trend: assuming variable as 3D [time,lat,lon]

#plotting SST annual mean 41 years of trend, stippling pvalue where it is significant at 95%

plt.figure(figsize=(10, 6))
plt.contourf(lon,lat,SST_pvalue,[0,0.05],colors='none',hatches='...')
plt.contourf(lon, lat, SST_trend*41,cmap='coolwarm',levels=np.arange(-2,2.1,.25),extend='both',alpha=.85)
plt.colorbar(label='degree C')
plt.title('Sea Surface Temperature Trend per 41 years')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.show()