0

I am very new to python and am working on the Jupyter notebook. I Have the following packages imported:

import warnings

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

import cartopy.crs as ccrs
import cartopy.mpl.ticker as cticker
from cartopy.util import add_cyclic_point

I have data with a monthly frequency and units per second and I'm trying to get them to be per year. Any help would be appreciated. I tried to resample the data before but have decided I cant rely on the inner functions of python.

#Constants
DayPerYear=360 # Calendars vary by modelling group
path='C:/mypath/'
PercipFile='myfile'
ds=xr.open_dataset(path+PercipFile)
ds 
Time_raw=ds.variables['time'][:]
print(Time_raw)

tt = np.size("time")
Time_m = np.zeros((tt))
SAT_m = np.zeros((tt))
ty = int(tt / 12)

# Time in years
Time_y = np.arange(1911, 1911 + ty, 1) + 0.5

SAT_y = np.zeros((ty))
yc = 0  # year count

Months = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
Months_leap = [31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]

for j in range(0, tt, 12):
    # Get year length, need to change if leap years
    if (j + 12 < tt):
        lyear = Time_raw[j + 12] - Time_raw[j]  # length of year
    else:
        lyear = Time_raw[-1] - Time_raw[j - 1]  # Works BC leap is in 
February

    # Set pre count variables
    YearSum = 0

    # Enter year loop
    for k in range(0, 12):
        if (lyear == 360):
            lmonth = 30
        elif (lyear == 365):
            lmonth = Months[k]
        elif (lyear == 366):
            lmonth = Months_leap[k]
        else:
            print('Year check')
            lmonth = 30  # set a default value
    

        YearSum = YearSum + SAT_m[j + k] * 1month

    # Exit month loop divide by year length
    SAT_y[yc] = YearSum / year

    # Increase year counter
    yc = yc + 1

Im getting IndexError: index 1 is out of bounds for axis 0 with size 1 on line 38 if that helps. thank you!

Violet Silva
  • 1
  • 1
  • Welcome to SO. Can you tell us what error is thrown and by what line? It looks like most of the code comes after the resampling and is irrelevant. Can you remove that code? Otherwise people will just assume the final line is throwing the error – Robert Wilson May 09 '23 at 06:59
  • Hello! I updated the code to what I've managed to get now as well as adding the error! Thanks for commenting, I completely forgot to add it! – Violet Silva May 09 '23 at 17:05
  • I don’t follow the code. Why are you defining the time steps based on the size of the character string ‘time’? – Robert Wilson May 09 '23 at 18:27

1 Answers1

0

I managed to solve it with some help from a friend. This is our solution:

#Import Modules
import warnings

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

#Define the location and name of the netCDF file
path='C:/MyPath/'
PercipFile='MyFile.nc'
ds=xr.open_dataset(path+PercipFile)
ds 

DayPerYear=360 # Calendars vary by modelling group
StartYear=1911
Time_raw=ds.variables['time'][:]

# COMPUTE MONTHLY GLOBAL AVERAGE PR
  
tt=np.size(Time_raw)
Time_m=np.zeros((tt))
PR_m=np.zeros((tt))

    
#YearZero=np.floor(Time_raw[0]/DayPerYear)+1850
    
for ii in range(0,tt):
Time_m[ii]=ii*(1./12)+(1./24)+StartYear
PR_m[ii]=(np.sum(ds.variables['pr'][ii,:,:]))
        
#plt.figure(1)
plt.plot(Time_m,PR_m,'k')

# COMPUTE YEARLY GLOBAL AVERAGE PR
ty=tt/12
# Time in years
Time_y=np.arange(StartYear,StartYear+ty,1)+0.5

PR_y=np.zeros((int(ty)))   
yc=0 # year count


Months=[31,28,31,30,31,30,31,31,30,31,30,31]
Months_leap=[31,29,31,30,31,30,31,31,30,31,30,31]

for j in range(0,tt,12):
    
    # Get year length,need to change if leap years
    if(j+12<tt):
        lyear=Time_raw[j+12]-Time_raw[j] # length of year
    else:
        lyear=Time_raw[-1]-Time_raw[j-1] # Works BC leap is in Febuary
    
    
    # Set pre count variables
    YearSum=0
    
    #enter year loop
    for k in range(0,12):
        if(lyear==360):
            month1=30
        elif(lyear==365):
            month1=Months[k]
        elif(lyear==366): 
            month1=Months_leap[k]
        else:
            #print('Year check')
            month1=30
            #stop
            
        YearSum=YearSum+PR_m[j+k]*month1

    # exit month loop divide by year length
    PR_y[yc]=YearSum/(lyear.astype('timedelta64[s]').item() * (365.25 * 24 * 3600))
    #print("PR_y[{}]: {}".format(yc, PR_y[yc]))   
    # increase year counter
    yc=yc+1
    #print("yc: {}, PR_y: {}".format(yc, PR_y))
    #exit year loop
plt.figure(2)
plt.plot(Time_y,PR_y)
plt.ylabel('PR')
plt.xlabel('Time')
#print(PR_y)

This gave me the graphs I was looking for both the monthly and yearly averages.

Violet Silva
  • 1
  • 1