Daily data analysis

Reshape

To analyze the climate data, we sometimes want to reshape the dimensions of an array. Usually, climate data have only one dimension for the time coordinate. However, we can apply a simple analysis method by making two dimensions for the time coordinate: year and day. Once we define the Xarray DataArray or DataSet, it is more challenging to reshape the coordinate. In this example, we first reshape the data array based on NumPy and then define the Xarray to make a plot.

Download Sample Files

• daily_reshape-v2.py

• GHCN_sample.csv

Packages

This script uses the basic packages. The datetime package helps us to make the time coordinate easily.

import numpy as np
import pandas as pd
import xarray as xr
import os
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime

Read CSV file

This script gets climate data from the CSV file. Please check here for more details.

# ----- Parameter setting ------
# == Date ==
dateparse = lambda x: datetime.strptime(x, '%Y-%m-%d')
#dateparse = lambda x: datetime.strptime(x, '%Y-%m-%d %H:%M:%S')
#dateparse = lambda x: datetime.strptime(x, '%m/%d/%y %H:%M')

# == CSV file name and location"
fcsv = 'GHCN_sample.csv'

# == Your definition of variable names
## Day,Precipitation, Multi-Day Precipitation, Snow Depth, Snow Fall, Min Temperature, Max Temperature, Ref Evapotranspiration,
vnames = ['date','precip','mdpr','sdep','sfall','tmin','tmax','refe']  

# == location of header. It starts from 0
nhead = 13

# == Figure name ==
fnFIG = os.path.splitext(os.path.basename(__file__))[0]

# ------------------------------

print(' ----- Read CSV file ------')
# additional parameters: encoding = 'ISO-8859-1'
fin = pd.read_csv(fcsv, 
	 delimiter = ',', header = nhead,  na_values = ['M','T','S'],
         names = vnames, parse_dates = ['date'], date_parser = dateparse)
print(fin)
print(' ----------------------------')

Reshaping

This part is the main description for formating the data array. First, we define two coordinates for the time dimension: year and day. The array “date” comes from the CSV file (Line 40). The next line defines “yyyy”, which corresponds to the NumPy array from the first to the last years in the “date” array with a 1-year interval. The NumPy array “days” indicates a day of the year from 1 to 365 days. We ignore all leap days (i.e., Feb 29) here. We also define the NumPy array “mmdd” for mm-dd format (i.e., 05-01) to make the readable time axis. The important point here is that all coordinates are NumPy arrays.

# ---- Make coordinates ----
date = fin.date
yyyy = np.arange(date[0].year, date[date.size-1].year+1,dtype='int')
days = np.arange(1,366,dtype='int')
mmdd = pd.date_range(start='2001-01-01', periods=365, freq='D').strftime('%m-%d').to_numpy(dtype='str')
#print(yyyy)
#print(days)
#print(mmdd)

def reshape_data(date,series,yyyy,days,mmdd):
  ndat = np.empty((yyyy.size,days.size))
  ndat[:] = np.NaN
  for it in range(0, date.size):
    #print(' it = ',it,' date = ',date[it], ' values = ',series.iloc[it])
    iy, = np.where(yyyy == date[it].year)
    id, = np.where(days == date[it].day_of_year)
    if date[it].is_leap_year and date[it].day_of_year > 60:
       #print('  adjusting day-1 because of a leap year')
       id = id - 1
    ndat[iy,id] = series.iloc[it]

  odat = xr.DataArray(ndat, dims=('year','day'), coords={'year': yyyy, 'day': days, 'mmdd': ('day', mmdd)})
  odat.set_index(days="mmdd")
  return(odat)


tmax=reshape_data(date,fin.tmax,yyyy,days,mmdd)
tmin=reshape_data(date,fin.tmin,yyyy,days,mmdd)

print(tmax)
print(tmin)

Data analysis

#--- calculate climatology and standard deviations for daily basis
tmax_clm=tmax.mean(dim='year', skipna=True)
tmax_sd=tmax.std(dim='year', skipna=True)

Plot a figure

# ------- Draw Figure -------
fig, ax = plt.subplots()
ax.plot(mmdd,tmax_clm, '-', color='darkgray', linewidth=2)
ax.plot(mmdd,tmax_clm+tmax_sd, '-', color='darkgray', linewidth=1)
ax.plot(mmdd,tmax_clm-tmax_sd, '-', color='darkgray', linewidth=1)
ax.plot(mmdd,tmax.sel(year=2020), 'o', markersize = 3, color='C1')
ax.plot(mmdd,tmax.sel(year=2018), 'o', markersize = 3, color='C0')
#ax.plot(mmdd,tmin.sel(year=2020), 'o', markersize = 1)
#ax.plot(mmdd,tmin.sel(year=2018), 'o', markersize = 1)


ax.set_xlabel('Date')
ax.set_ylabel('Temperature [C]')
ax.set_title('Maximum daily air temperature')

# --- Uncomment out for the format of the ticks --
ax.xaxis.set_major_locator(mdates.MonthLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))
ax.xaxis.set_minor_locator(mdates.DayLocator())
# ------------------------

ax.legend(['Clm','+1 std','-1 std', '2020', '2018'])
#ax.legend(['Tmax (clm)','Tmax (2020)','Tmin (2020)','Tmax (2018)', 'Tmin (2018)'])
#ax.axis(xmin='05-01', xmax='08-01')
ax.axis(xmin='01-01', xmax='12-31')
ax.axis(ymin=-10, ymax=45)
ax.grid(True)
fig.autofmt_xdate()


plt.draw()
plt.savefig(fnFIG+".png")