From c6c5791cf246c423654d53024b8b73e50a5b5226 Mon Sep 17 00:00:00 2001 From: Kaitlin Naughten Date: Wed, 15 Mar 2017 16:57:16 +1100 Subject: [PATCH] Common code for monthly and seasonal averages; put common functions in their own file instead of copy-paste --- adv_amery_tsplots_indiv.py | 193 --------------------- adv_ross_tsplots.py | 112 +------------ aice_hi_seasonal.py | 170 ++----------------- aice_hs_seasonal.py | 123 +------------- cice_evap_era_monthly.py | 232 -------------------------- compare_nic_seasonal.py | 196 ++-------------------- interp_era2roms.py | 54 ++++++ interp_lon_roms.py | 106 ++++++++++++ interp_lon_sose.py | 37 +++++ monthly_avg_cice.py | 117 +++++++++++++ monthly_avg_roms.py | 112 +++++++++++++ nsidc_aice_monthly.py | 134 ++------------- nsidc_aice_seasonal.py | 166 ++----------------- romscice_atm_monthly.py | 52 +----- romscice_atm_subdaily.py | 63 +------ romscice_evap.py | 53 +----- romscice_ini_ecco.py | 8 +- romscice_nbc.py | 10 +- seasonal_avg_cice.py | 163 ++++++++++++++++++ seasonal_avg_roms.py | 180 ++++++++++++++++++++ sose_roms_seasonal.py | 315 ++--------------------------------- sose_seasonal.py | 38 +---- ssflux_monthly.py | 102 +----------- ssflux_tamura_monthly.py | 165 +++--------------- ssflux_tamura_seasonal.py | 198 ++++------------------ temp_salt_seasonal.py | 331 ++----------------------------------- temp_salt_slice.py | 109 +----------- timeseries_dpt.py | 2 +- zonal_plot.py | 106 +----------- 29 files changed, 939 insertions(+), 2708 deletions(-) delete mode 100644 adv_amery_tsplots_indiv.py delete mode 100644 cice_evap_era_monthly.py create mode 100644 interp_era2roms.py create mode 100644 interp_lon_roms.py create mode 100644 interp_lon_sose.py create mode 100644 monthly_avg_cice.py create mode 100644 monthly_avg_roms.py create mode 100644 seasonal_avg_cice.py create mode 100644 seasonal_avg_roms.py diff --git a/adv_amery_tsplots_indiv.py b/adv_amery_tsplots_indiv.py deleted file mode 100644 index e3955f8..0000000 --- a/adv_amery_tsplots_indiv.py +++ /dev/null @@ -1,193 +0,0 @@ -from netCDF4 import Dataset -from numpy import * -from matplotlib.pyplot import * -from calc_z import * - -# For each advection experiment, plot zonal slices of temperature and salinity -# through 71E (Amery Ice Shelf) at the end of the simulation. -def adv_amery_tsplots_indiv (): - - num_simulations = 6 - # Paths to simulation directories - paths = ['/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/c4_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/c4_highdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/a4_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/a4_highdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/u3_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/u3limiters_lowdif/'] - # End of figure names for each simulation - labels = ['_c4_lowdif.png', '_c4_highdif.png', '_a4_lowdif.png', '_a4_highdif.png', '_u3.png', '_u3_lim.png'] - # Name of ocean output file to read - ocn_file = 'ocean_avg_0001.nc' - # Timestep to plot (average over last day in 1992) - tstep = 366 - # Longitude to plot - lon0 = 71 - # Deepest depth to plot - depth_min = -500 - # Bounds on colour scale for each variable - temp_bounds = [-2, 3] - salt_bounds = [33.8, 34.8] - # Bounds on latitudes to plot - lat_min = -72 - lat_max = -50 - - # Grid parameters - theta_s = 4.0 - theta_b = 0.9 - hc = 40 - N = 31 - - # Build titles for each variable based on longitude - if lon0 < 0: - temp_title = r'Temperature ($^{\circ}$C) at ' + str(int(round(-lon0))) + r'$^{\circ}$W' - salt_title = r'Salinity (psu) at ' + str(int(round(-lon0))) + r'$^{\circ}$W' - # Edit longitude to be between0 and 360, following ROMS convention - lon0 += 360 - else: - temp_title = r'Temperature ($^{\circ}$C) at ' + str(int(round(lon0))) + r'$^{\circ}$E' - salt_title = r'Salinity (psu) at ' + str(int(round(lon0))) + r'$^{\circ}$E' - - # Loop over simulations - for sim in range(num_simulations): - # Loop over variables - for var_name in ['temp', 'salt']: - # Read variable, sea surface height, and grid variables - id = Dataset(paths[sim] + ocn_file, 'r') - data_3d = id.variables[var_name][tstep-1,:,:-15,:] - zeta = id.variables['zeta'][tstep-1,:-15,:] - if sim == 0 and var_name == 'temp': - # Grid variables are the same for all simulations so we - # only need to read them once - h = id.variables['h'][:-15,:] - zice = id.variables['zice'][:-15,:] - lon_2d = id.variables['lon_rho'][:-15,:] - lat_2d = id.variables['lat_rho'][:-15,:] - id.close() - # Get a 3D array of z-coordinates - z_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N, zeta) - # Interpolate the variable, z, and latitude to lon0 - data, z, lat = interp_lon(data_3d, z_3d, lat_2d, lon_2d, lon0) - # Set up colour levels for plotting - if var_name == 'temp': - lev = linspace(temp_bounds[0], temp_bounds[1], num=40) - elif var_name == 'salt': - lev = linspace(salt_bounds[0], salt_bounds[1], num=40) - # Plot - fig = figure(figsize=(12,6)) - contourf(lat, z, data, lev, cmap='jet', extend='both') - colorbar() - if var_name == 'temp': - title(temp_title) - elif var_name == 'salt': - title(salt_title) - xlabel('Latitude') - ylabel('Depth (m)') - xlim([lat_min, lat_max]) - ylim([depth_min, 0]) - # Save plot - fig.savefig(var_name + labels[sim]) - - -# Linearly interpolate data, z, and latitude to the specified longitude. -# Input: -# data_3d = array of data, dimension depth x lat x lon -# z_3d = array of depth values (negative, in metres), dimension depth x lat x lon -# lat_2d = array of latitudevalues, dimension lat x lon -# lon_2d = array of longitude values, dimension lat x lon (between -180 and 180) -# lon0 = longitude to interpolate to (between -180 and 180) -# Output: -# data = array of data interpolated to lon0, dimension depth x lat -# z = array of depth values interpolated to lon0, dimension depth x lat -# lat = array of latitude values interpolated to lon0, dimension depth x lat -def interp_lon (data_3d, z_3d, lat_2d, lon_2d, lon0): - - # Save dimensions - num_depth = size(data_3d, 0) - num_lat = size(data_3d, 1) - num_lon = size(data_3d, 2) - # Set up output arrays - data = ma.empty([num_depth, num_lat]) - z = ma.empty([num_depth, num_lat]) - lat = ma.empty([num_depth, num_lat]) - - # Loop over latitudes; can't find a cleaner way to do this - for j in range(num_lat): - # Extract the longitude values of this slice - lon_tmp = lon_2d[j,:] - # Get indices and coefficients for interpolation - ie, iw, coeffe, coeffw = interp_lon_helper(lon_tmp, lon0) - data[:,j] = coeffe*data_3d[:,j,ie] + coeffw*data_3d[:,j,iw] - z[:,j] = coeffe*z_3d[:,j,ie] + coeffw*z_3d[:,j,iw] - lat[:,j] = coeffe*lat_2d[j,ie] + coeffw*lat_2d[j,iw] - - return data, z, lat - - -# Calculate indices and coefficients for linear interpolation of longitude. -# This takes care of all the mod 360 nonsense. -# Input: -# lon = 1D array of longitude values (straight out of ROMS i.e. between slightly < 0 and slightly > 360) -# lon0 = longitude to interpolate to (between 0 and 360) -# Output: -# ie, iw, coeffe, coeffw = integers (ie and iw) and coefficients (coeffe and coeffw) such that -# coeffe*lon[ie] + coeffw*lon[iw] = lon0, which will also hold for any -# variable on this longitude grid. ie is the index of the nearest point -# to the east of lon0; iw the nearest point to the west. -def interp_lon_helper (lon, lon0): - - if lon0 < amin(lon) or lon0 > amax(lon): - # Special case: lon0 on periodic boundary - # Be careful with mod 360 here - - # Find the periodic boundary - dlon = lon[1:] - lon[0:-1] - bdry = argmax(abs(dlon)) - if dlon[bdry] < -300: - # Jumps from almost 360 to just over 0 - iw = bdry - ie = bdry + 1 - else: - # Periodic boundary lines up with the array boundary - iw = size(lon) - 1 - ie = 0 - # Calculate difference between lon0 and lon[iw], mod 360 if necessary - dlon_num = lon0 - lon[iw] - if dlon_num < -300: - dlon_num += 360 - # Calculate difference between lon[ie] and lon[iw], mod 360 - dlon_den = lon[ie] - lon[iw] + 360 - - else: - # General case - - # Add or subtract 360 from longitude values which wrap around - # so that longitude increases monotonically from west to east - i = arange(1, size(lon)+1) - index1 = nonzero((i > 1200)*(lon < 100)) - lon[index1] = lon[index1] + 360 - index2 = nonzero((i < 200)*(lon > 300)) - lon[index2] = lon[index2] - 360 - - # Take mod 360 of lon0 if necessary - if all(lon < lon0): - lon0 -= 360 - if all(lon > lon0): - lon0 += 360 - - # Find the first index eastward of lon0 - ie = nonzero(lon > lon0)[0][0] - # The index before it will be the last index westward of lon0 - iw = ie - 1 - - dlon_num = lon0 - lon[iw] - dlon_den = lon[ie] - lon[iw] - - if dlon_num > 5 or dlon_den > 5: - print 'interp_lon_helper: Problem at periodic boundary' - return - coeff1 = dlon_num/dlon_den - coeff2 = 1 - coeff1 - - return ie, iw, coeff1, coeff2 - - -# Command-line interface -if __name__ == "__main__": - - adv_amery_tsplots_indiv() diff --git a/adv_ross_tsplots.py b/adv_ross_tsplots.py index 9da35a2..a50788c 100644 --- a/adv_ross_tsplots.py +++ b/adv_ross_tsplots.py @@ -2,6 +2,7 @@ from numpy import * from matplotlib.pyplot import * from calc_z import * +from interp_lon_roms import * # Create a 2x2 plot showing zonal slices of temperature and salinity through # 180E (Ross Sea) at the end of the U3_LIM and C4_LD simulations. @@ -54,7 +55,7 @@ def adv_ross_tsplots (): # Calculate 3D z-coordinates z_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N, zeta) # Interpolate temp/salt, z-coordinates, and latitude to 180E - data, z, lat = interp_lon(data_3d, z_3d, lat_2d, lon_2d, lon0) + data, z, lat = interp_lon_roms(data_3d, z_3d, lat_2d, lon_2d, lon0) ax = fig.add_subplot(2, 2, 2*var+sim+1) # Shade data (pcolor not contourf so we don't misrepresent the # model grid) @@ -92,113 +93,8 @@ def adv_ross_tsplots (): # Main title suptitle(r'180$^{\circ}$E, 31 December (Ross Sea)', fontsize=30) - #fig.show() - fig.savefig('adv_ross_tsplots.png') - - - - -# Linearly interpolate data, z, and latitude to the specified longitude. -# Input: -# data_3d = array of data, dimension depth x lat x lon -# z_3d = array of depth values (negative, in metres), dimension depth x lat x lon -# lat_2d = array of latitudevalues, dimension lat x lon -# lon_2d = array of longitude values, dimension lat x lon (between -180 and 180) -# lon0 = longitude to interpolate to (between -180 and 180) -# Output: -# data = array of data interpolated to lon0, dimension depth x lat -# z = array of depth values interpolated to lon0, dimension depth x lat -# lat = array of latitude values interpolated to lon0, dimension depth x lat -def interp_lon (data_3d, z_3d, lat_2d, lon_2d, lon0): - - # Save dimensions - num_depth = size(data_3d, 0) - num_lat = size(data_3d, 1) - num_lon = size(data_3d, 2) - # Set up output arrays - data = ma.empty([num_depth, num_lat]) - z = ma.empty([num_depth, num_lat]) - lat = ma.empty([num_depth, num_lat]) - - # Loop over latitudes; can't find a cleaner way to do this - for j in range(num_lat): - # Extract the longitude values of this slice - lon_tmp = lon_2d[j,:] - # Get indices and coefficients for interpolation - ie, iw, coeffe, coeffw = interp_lon_helper(lon_tmp, lon0) - data[:,j] = coeffe*data_3d[:,j,ie] + coeffw*data_3d[:,j,iw] - z[:,j] = coeffe*z_3d[:,j,ie] + coeffw*z_3d[:,j,iw] - lat[:,j] = coeffe*lat_2d[j,ie] + coeffw*lat_2d[j,iw] - - return data, z, lat - - -# Calculate indices and coefficients for linear interpolation of longitude. -# This takes care of all the mod 360 nonsense. -# Input: -# lon = 1D array of longitude values (straight out of ROMS i.e. between slightly < 0 and slightly > 360) -# lon0 = longitude to interpolate to (between 0 and 360) -# Output: -# ie, iw, coeffe, coeffw = integers (ie and iw) and coefficients (coeffe and coeffw) such that -# coeffe*lon[ie] + coeffw*lon[iw] = lon0, which will also hold for any -# variable on this longitude grid. ie is the index of the nearest point -# to the east of lon0; iw the nearest point to the west. -def interp_lon_helper (lon, lon0): - - if lon0 < amin(lon) or lon0 > amax(lon): - # Special case: lon0 on periodic boundary - # Be careful with mod 360 here - - # Find the periodic boundary - dlon = lon[1:] - lon[0:-1] - bdry = argmax(abs(dlon)) - if dlon[bdry] < -300: - # Jumps from almost 360 to just over 0 - iw = bdry - ie = bdry + 1 - else: - # Periodic boundary lines up with the array boundary - iw = size(lon) - 1 - ie = 0 - # Calculate difference between lon0 and lon[iw], mod 360 if necessary - dlon_num = lon0 - lon[iw] - if dlon_num < -300: - dlon_num += 360 - # Calculate difference between lon[ie] and lon[iw], mod 360 - dlon_den = lon[ie] - lon[iw] + 360 - - else: - # General case - - # Add or subtract 360 from longitude values which wrap around - # so that longitude increases monotonically from west to east - i = arange(1, size(lon)+1) - index1 = nonzero((i > 1200)*(lon < 100)) - lon[index1] = lon[index1] + 360 - index2 = nonzero((i < 200)*(lon > 300)) - lon[index2] = lon[index2] - 360 - - # Take mod 360 of lon0 if necessary - if all(lon < lon0): - lon0 -= 360 - if all(lon > lon0): - lon0 += 360 - - # Find the first index eastward of lon0 - ie = nonzero(lon > lon0)[0][0] - # The index before it will be the last index westward of lon0 - iw = ie - 1 - - dlon_num = lon0 - lon[iw] - dlon_den = lon[ie] - lon[iw] - - if dlon_num > 5 or dlon_den > 5: - print 'interp_lon_helper: Problem at periodic boundary' - return - coeff1 = dlon_num/dlon_den - coeff2 = 1 - coeff1 - - return ie, iw, coeff1, coeff2 + fig.show() + #fig.savefig('adv_ross_tsplots.png') # Command-line interface diff --git a/aice_hi_seasonal.py b/aice_hi_seasonal.py index 0d5220e..a3d8565 100644 --- a/aice_hi_seasonal.py +++ b/aice_hi_seasonal.py @@ -1,6 +1,7 @@ from numpy import * -from netCDF4 import Dataset, num2date +from netCDF4 import Dataset from matplotlib.pyplot import * +from seasonal_avg_cice import * # Creates a 4x2 plot of seasonally averaged sea ice concentration (top row) and # thickness (bottom row) over the last year of simulation. @@ -13,24 +14,18 @@ # fig_name = if save=True, path to the desired filename for figure def aice_hi_seasonal (cice_file, save=False, fig_name=None): - # Starting and ending months (1-based) for each season - start_month = [12, 3, 6, 9] - end_month = [2, 5, 8, 11] - # Starting and ending days of the month (1-based) for each season - # Assume no leap years, we'll fix this later if we need - start_day = [1, 1, 1, 1] - end_day = [28, 31, 31, 30] - # Number of days in each season (again, ignore leap years for now) - ndays_season = [90, 92, 92, 91] # Season names for titles season_names = ['DJF', 'MAM', 'JJA', 'SON'] # Degrees to radians conversion deg2rad = pi/180.0 - # Read CICE grid and time values + # Read the CICE grid id = Dataset(cice_file, 'r') lon_tmp = id.variables['TLON'][:-15,:] lat_tmp = id.variables['TLAT'][:-15,:] + num_lon = id.variables['TLON'].shape[1] + num_lat = id.variables['TLAT'].shape[0] + id.close() # Wrap the periodic boundary by 1 cell lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1]) lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1]) @@ -38,154 +33,13 @@ def aice_hi_seasonal (cice_file, save=False, fig_name=None): lon[:,-1] = lon_tmp[:,0] lat[:,:-1] = lat_tmp lat[:,-1] = lat_tmp[:,0] - time_id = id.variables['time'] - # Get the year, month, and day (all 1-based) for each output step - # These are 5-day averages marked with the next day's date. - time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - - # Loop backwards through time indices to find the last one we care about - # (which contains 30 November in its averaging period) - end_t = -1 # Missing value flag - for t in range(size(time)-1, -1, -1): - if time[t].month == start_month[0] and time[t].day in range(start_day[0], start_day[0]+5): - end_t = t - break - # Make sure we actually found it - if end_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' - return - - # Continue looping backwards to find the first time index we care about - # (which contains 1 December the previous year in its averaging period) - start_t = -1 # Missing value flag - for t in range(end_t-60, -1, -1): - if time[t].month == start_month[0] and time[t].day in range(start_day[0]+1, start_day[0]+6): - start_t = t - break - # Make sure we actually found it - if start_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' - return - - # Check for leap years - leap_year = False - if mod(time[end_t].year, 4) == 0: - # Years divisible by 4 are leap years - leap_year = True - if mod(time[end_t].year, 100) == 0: - # Unless they're also divisible by 100, in which case they aren't - # leap years - leap_year = False - if mod(time[end_t].year, 400) == 0: - # Unless they're also divisible by 400, in which case they are - # leap years after all - leap_year = True - if leap_year: - # Update last day in February - end_day[0] += 1 - ndays_season[0] += 1 - - # Initialise seasonal averages of CICE output - aice_tmp = ma.empty([4, size(lon_tmp,0), size(lon_tmp,1)]) - aice_tmp[:,:,:] = 0.0 - hi_tmp = ma.empty([4, size(lon_tmp,0), size(lon_tmp,1)]) - hi_tmp[:,:,:] = 0.0 - # Process one season at a time - for season in range(4): - season_days = 0 # Number of days in season; this will be incremented - next_season = mod(season+1, 4) - - # Find starting timestep - start_t_season = -1 - for t in range(start_t, end_t+1): - if time[t].month == start_month[season] and time[t].day in range(start_day[season]+1, start_day[season]+6): - start_t_season = t - break - # Make sure we actually found it - if start_t_season == -1: - print 'Error: could not find starting timestep for season ' + season_names[season] - return - - # Find ending timestep - end_t_season = -1 - for t in range(start_t_season+1, end_t+1): - if time[t].month == start_month[next_season] and time[t].day in range(start_day[next_season], start_day[next_season]+5): - end_t_season = t - break - # Make sure we actually found it - if end_t_season == -1: - print 'Error: could not find ending timestep for season ' + season_names[season] - return - - # Figure out how many of the 5 days averaged in start_t_season are - # actually within this season - if time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 5: - # Starting day is in position 1 of 5; we care about all of them - start_days = 5 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 4: - # Starting day is in position 2 of 5; we care about the last 4 - start_days = 4 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]+ 3: - # Starting day is in position 3 of 5; we care about the last 3 - start_days = 3 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 2: - # Starting day is in position 4 of 5; we care about the last 2 - start_days = 2 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 1: - # Starting day is in position 5 of 5; we care about the last 1 - start_days = 1 - else: - print 'Error for season ' + season_names[season] + ': starting index is month ' + str(time[start_t_season].month) + ', day ' + str(time[start_t_season].day) - return - - # Start accumulating data weighted by days - aice_tmp[season,:,:] += id.variables['aice'][start_t_season,:-15,:]*start_days - hi_tmp[season,:,:] += id.variables['hi'][start_t_season,:-15,:]*start_days - season_days += start_days - - # Between start_t_season and end_t_season, we want all the days - for t in range(start_t_season+1, end_t_season): - aice_tmp[season,:,:] += id.variables['aice'][t,:-15,:]*5 - hi_tmp[season,:,:] += id.variables['hi'][t,:-15,:]*5 - season_days += 5 - - # Figure out how many of the 5 days averaged in end_t_season are - # actually within this season - if time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 4: - # Ending day is in position 1 of 5; we care about the first 1 - end_days = 1 - elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 3: - # Ending day is in position 2 of 5; we care about the first 2 - end_days = 2 - elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 2: - # Ending day is in position 3 of 5; we care about the first 3 - end_days = 3 - elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 1: - # Ending day is in position 4 of 5; we care about the first 4 - end_days = 4 - elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season]: - # Ending day is in position 5 of 5; we care about all 5 - end_days = 5 - else: - print 'Error for season ' + season_names[season] + ': ending index is month ' + str(time[end_t_season].month) + ', day ' + str(time[end_t_season].day) - return - - aice_tmp[season,:,:] += id.variables['aice'][end_t_season,:-15,:]*end_days - hi_tmp[season,:,:] += id.variables['hi'][end_t_season,:-15,:]*end_days - season_days += end_days - - # Check that we got the correct number of days - if season_days != ndays_season[season]: - print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) - return - - # Finished accumulating data, now convert from sum to average - aice_tmp[season,:,:] /= season_days - hi_tmp[season,:,:] /= season_days - - # Finished reading all CICE data - id.close() + # Read seasonally averaged fields + aice_tmp = seasonal_avg_cice(cice_file, 'aice', [num_lat, num_lon]) + hi_tmp = seasonal_avg_cice(cice_file, 'hi', [num_lat, num_lon]) + # Chop off northern boundary + aice_tmp = aice_tmp[:,:-15,:] + hi_tmp = hi_tmp[:,:-15,:] # Wrap the periodic boundary aice = ma.empty([size(aice_tmp,0), size(aice_tmp,1), size(aice_tmp,2)+1]) aice[:,:,:-1] = aice_tmp diff --git a/aice_hs_seasonal.py b/aice_hs_seasonal.py index c5811ae..ea9c486 100644 --- a/aice_hs_seasonal.py +++ b/aice_hs_seasonal.py @@ -1,134 +1,30 @@ from numpy import * -from netCDF4 import Dataset, num2date +from netCDF4 import Dataset from matplotlib.pyplot import * +from seasonal_avg_cice import * def aice_hs_seasonal (cice_file, save=False, fig_name=None): - start_month = [12, 3, 6, 9] - end_month = [2, 5, 8, 11] - start_day = [1, 1, 1, 1] - end_day = [28, 31, 31, 30] - ndays_season = [90, 92, 92, 91] season_names = ['DJF', 'MAM', 'JJA', 'SON'] deg2rad = pi/180.0 id = Dataset(cice_file, 'r') lon_tmp = id.variables['TLON'][:-15,:] lat_tmp = id.variables['TLAT'][:-15,:] + num_lon = id.variables['TLON'].shape[1] + num_lat = id.variables['TLAT'].shape[0] + id.close() lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1]) lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1]) lon[:,:-1] = lon_tmp lon[:,-1] = lon_tmp[:,0] lat[:,:-1] = lat_tmp lat[:,-1] = lat_tmp[:,0] - time_id = id.variables['time'] - time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - - end_t = -1 - for t in range(size(time)-1, -1, -1): - if time[t].month == start_month[0] and time[t].day in range(start_day[0], start_day[0]+5): - end_t = t - break - if end_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' - return - - start_t = -1 - for t in range(end_t-60, -1, -1): - if time[t].month == start_month[0] and time[t].day in range(start_day[0]+1, start_day[0]+6): - start_t = t - break - if start_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' - return - - leap_year = False - if mod(time[end_t].year, 4) == 0: - leap_year = True - if mod(time[end_t].year, 100) == 0: - leap_year = False - if mod(time[end_t].year, 400) == 0: - leap_year = True - if leap_year: - end_day[0] += 1 - ndays_season[0] += 1 - - aice_tmp = ma.empty([4, size(lon_tmp,0), size(lon_tmp,1)]) - aice_tmp[:,:,:] = 0.0 - hs_tmp = ma.empty([4, size(lon_tmp,0), size(lon_tmp,1)]) - hs_tmp[:,:,:] = 0.0 - for season in range(4): - season_days = 0 - next_season = mod(season+1, 4) - - start_t_season = -1 - for t in range(start_t, end_t+1): - if time[t].month == start_month[season] and time[t].day in range(start_day[season]+1, start_day[season]+6): - start_t_season = t - break - if start_t_season == -1: - print 'Error: could not find starting timestep for season ' + season_names[season] - return - - end_t_season = -1 - for t in range(start_t_season+1, end_t+1): - if time[t].month == start_month[next_season] and time[t].day in range(start_day[next_season], start_day[next_season]+5): - end_t_season = t - break - if end_t_season == -1: - print 'Error: could not find ending timestep for season ' + season_names[season] - return - - if time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 5: - start_days = 5 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 4: - start_days = 4 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]+ 3: - start_days = 3 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 2: - start_days = 2 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 1: - start_days = 1 - else: - print 'Error for season ' + season_names[season] + ': starting index is month ' + str(time[start_t_season].month) + ', day ' + str(time[start_t_season].day) - return - - aice_tmp[season,:,:] += id.variables['aice'][start_t_season,:-15,:]*start_days - hs_tmp[season,:,:] += id.variables['hs'][start_t_season,:-15,:]*start_days - season_days += start_days - - for t in range(start_t_season+1, end_t_season): - aice_tmp[season,:,:] += id.variables['aice'][t,:-15,:]*5 - hs_tmp[season,:,:] += id.variables['hs'][t,:-15,:]*5 - season_days += 5 - - if time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 4: - end_days = 1 - elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 3: - end_days = 2 - elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 2: - end_days = 3 - elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 1: - end_days = 4 - elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season]: - end_days = 5 - else: - print 'Error for season ' + season_names[season] + ': ending index is month ' + str(time[end_t_season].month) + ', day ' + str(time[end_t_season].day) - return - - aice_tmp[season,:,:] += id.variables['aice'][end_t_season,:-15,:]*end_days - hs_tmp[season,:,:] += id.variables['hs'][end_t_season,:-15,:]*end_days - season_days += end_days - - if season_days != ndays_season[season]: - print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) - return - - aice_tmp[season,:,:] /= season_days - hs_tmp[season,:,:] /= season_days - - id.close() + aice_tmp = seasonal_avg_cice(cice_file, 'aice', [num_lat, num_lon]) + hs_tmp = seasonal_avg_cice(cice_file, 'hs', [num_lat, num_lon]) + aice_tmp = aice_tmp[:,:-15,:] + hs_tmp = hs_tmp[:,:-15,:] aice = ma.empty([size(aice_tmp,0), size(aice_tmp,1), size(aice_tmp,2)+1]) aice[:,:,:-1] = aice_tmp aice[:,:,-1] = aice_tmp[:,:,0] @@ -139,7 +35,6 @@ def aice_hs_seasonal (cice_file, save=False, fig_name=None): x = -(lat+90)*cos(lon*deg2rad+pi/2) y = (lat+90)*sin(lon*deg2rad+pi/2) - # Find boundaries for each side of plot based on extent of NSIDC grid bdry1 = -35 bdry2 = 35 bdry3 = -33 diff --git a/cice_evap_era_monthly.py b/cice_evap_era_monthly.py deleted file mode 100644 index ca17bdd..0000000 --- a/cice_evap_era_monthly.py +++ /dev/null @@ -1,232 +0,0 @@ -from numpy import * -from netCDF4 import Dataset, num2date -from matplotlib.pyplot import * - -def cice_evap_era_monthly (cice_file, month, colour_bounds=None, save=False, fig_name=None): - - era_head = '/short/y99/kaa561/CMIP5_forcing/atmos/climatology/ERA_Interim_monthly/ER_' - era_tail = '_monthly_orig.nc' - # Starting and ending days in each month - # Assume no leap years, we'll fix this later if we need - start_day = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] - end_day = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] - # Name of each month, for the title - month_name = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] - # Degrees to radians conversion - deg2rad = pi/180.0 - # Conversion factor: m/12h to cm/day - era_conv = 100*2 - - # Read the CICE grid and time values - id = Dataset(cice_file, 'r') - cice_lon_tmp = id.variables['TLON'][:-15,:] - cice_lat_tmp = id.variables['TLAT'][:-15,:] - # Wrap the periodic boundary by 1 cell - cice_lon = ma.empty([size(cice_lon_tmp,0), size(cice_lon_tmp,1)+1]) - cice_lat = ma.empty([size(cice_lat_tmp,0), size(cice_lat_tmp,1)+1]) - cice_lon[:,:-1] = cice_lon_tmp - cice_lon[:,-1] = cice_lon_tmp[:,0] - cice_lat[:,:-1] = cice_lat_tmp - cice_lat[:,-1] = cice_lat_tmp[:,0] - time_id = id.variables['time'] - # Get the year, month, and day (all 1-indexed) for each output step - # These are 5-day averages marked with the next day's date. - cice_time = num2date(time_id[:], units=time_id.units, calendar='standard') - - next_month = mod(month+1,12) - - end_t = -1 - for t in range(size(cice_time)-1, -1, -1): - if cice_time[t].month-1 == next_month and cice_time[t].day in range(start_day[next_month], start_day[next_month]+5): - end_t = t - break - if end_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete ' + month_name[month] - return - - start_t = -1 - for t in range(end_t, -1, -1): - if cice_time[t].month-1 == month and cice_time[t].day in range(start_day[month]+1, start_day[month]+6): - start_t = t - break - if start_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete ' + month_name[month] - return - - leap_year = False - cice_year = cice_time[end_t].year - if month == 11: - cice_year = cice_time[start_t].year - if mod(cice_year, 4) == 0: - leap_year = True - if mod(cice_year, 100) == 0: - leap_year = False - if mod(cice_year, 400) == 0: - leap_year = True - if leap_year: - end_day[1] = 29 - - cice_evap_tmp = ma.empty(shape(cice_lon_tmp)) - cice_evap_tmp[:,:] = 0.0 - num_days = 0 - - # Figure out how many of the 5 days averaged in start_t are actually within - # this month - if cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+5: - # Starting day is in position 1 of 5; we care about all of them - start_days = 5 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+4: - # Starting day is in position 2 of 5; we care about the last 4 - start_days = 4 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+3: - # Starting day is in position 3 of 5; we care about the last 3 - start_days = 3 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+2: - # Starting day is in position 4 of 5; we care about the last 2 - start_days = 2 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+1: - # Starting day is in position 5 of 5; we care about the last 1 - start_days = 1 - else: - print 'Error: starting index is month ' + str(cice_time[start_t].month) + ', day ' + str(cice_time[start_t].day) - return - - cice_evap_tmp += id.variables['evap_ai'][start_t,:-15,:]*start_days - num_days += start_days - - for t in range(start_t+1, end_t): - cice_evap_tmp += id.variables['evap_ai'][t,:-15,:]*5 - num_days += 5 - - # Figure out how many of the 5 days averaged in end_t are actually within - # this month - if cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+4: - # Ending day is in position 1 of 5; we care about the first 1 - end_days = 1 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+3: - # Ending day is in position 2 of 5; we care about the first 2 - end_days = 2 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+2: - # Ending day is in position 3 of 5; we care about the first 3 - end_days = 3 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+1: - # Ending day is in position 4 of 5; we care about the first 4 - end_days = 4 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]: - # Ending day is in position 5 of 5; we care about all 5 - end_days = 5 - else: - print 'Error: ending index is month ' + str(cice_time[end_t].month) + ', day ' + str(cice_time[end_t].day) - return - - cice_evap_tmp += id.variables['evap_ai'][end_t,:-15,:]*end_days - num_days += end_days - - if num_days != end_day[month]: - print 'Error: found ' + str(num_days) + ' days instead of ' + str(end_day[month]) - return - - id.close() - cice_evap_tmp = -1*cice_evap_tmp/num_days - - cice_evap = ma.empty([size(cice_evap_tmp,0), size(cice_evap_tmp,1)+1]) - cice_evap[:,:-1] = cice_evap_tmp - cice_evap[:,-1] = cice_evap_tmp[:,0] - - id = Dataset(era_head + str(cice_year) + era_tail, 'r') - era_lon_1d = id.variables['longitude'][:] - era_lat_1d = id.variables['latitude'][:] - era_lon, era_lat = meshgrid(era_lon_1d, era_lat_1d) - era_evap = id.variables['e'][month,:,:]*era_conv - id.close() - era_evap = -1*era_evap - - cice_x = -(cice_lat+90)*cos(cice_lon*deg2rad+pi/2) - cice_y = (cice_lat+90)*sin(cice_lon*deg2rad+pi/2) - era_x = -(era_lat+90)*cos(era_lon*deg2rad+pi/2) - era_y = (era_lat+90)*sin(era_lon*deg2rad+pi/2) - - if colour_bounds is not None: - lev = linspace(colour_bounds[0], colour_bounds[1], num=50) - else: - min_bound = min(amin(cice_evap), amin(era_evap)) - max_bound = max(amax(cice_evap), amax(era_evap)) - lev = linspace(min_bound, max_bound, num=50) - bdry = -50+90 - - fig = figure(figsize=(20,9)) - fig.add_subplot(1,2,1, aspect='equal') - contourf(era_x, era_y, era_evap, lev, extend='both') - title('ERA-Interim', fontsize=24) - xlim([-bdry, bdry]) - ylim([-bdry, bdry]) - axis('off') - fig.add_subplot(1,2,2, aspect='equal') - img = contourf(cice_x, cice_y, cice_evap, lev, extend='both') - title('CICE', fontsize=24) - xlim([-bdry, bdry]) - ylim([-bdry, bdry]) - axis('off') - cbaxes = fig.add_axes([0.3, 0.04, 0.4, 0.04]) - cbar = colorbar(img, orientation='horizontal', cax=cbaxes) - suptitle(r'Evaporation (cm/day), ' + month_name[month] + ' ' + str(cice_year), fontsize=30) - - if save: - fig.savefig(fig_name) - else: - fig.show() - - -if __name__ == "__main__": - - cice_file = raw_input("Path to CICE file: ") - month = int(raw_input("Month number (1-12): "))-1 - colour_bounds = None - get_bounds = raw_input("Set bounds on colour scale (y/n)? ") - if get_bounds == 'y': - lower_bound = float(raw_input("Lower bound: ")) - upper_bound = float(raw_input("Upper bound: ")) - colour_bounds = [lower_bound, upper_bound] - action = raw_input("Save figure (s) or display on screen (d)? ") - if action == 's': - save = True - fig_name = raw_input("File name for figure: ") - elif action == 'd': - save = False - fig_name = None - cice_evap_era_monthly(cice_file, month, colour_bounds, save, fig_name) - - while True: - # Repeat until the user wants to exit - repeat = raw_input("Make another plot (y/n)? ") - if repeat == 'y': - while True: - # Ask for changes to the parameters until the user is done - changes = raw_input("Enter a parameter to change: (1) file path, (2) month number, (3) colour bounds, (4) save/display; or enter to continue: ") - if len(changes) == 0: - # No more changes to parameters - break - else: - if int(changes) == 1: - # New CICE file - cice_file = raw_input("Path to CICE file: ") - elif int(changes) == 2: - # New month - month = int(raw_input("Month number (1-12): "))-1 - elif int(changes) == 3: - colour_bounds = None - get_bounds = raw_input("Set bounds on colour scale (y/n)? ") - if get_bounds == 'y': - lower_bound = float(raw_input("Lower bound: ")) - upper_bound = float(raw_input("Upper bound: ")) - colour_bounds = [lower_bound, upper_bound] - elif int(changes) == 4: - # Switch from save to display, or vice versa - save = not save - if save: - # Get a new figure name - fig_name = raw_input("File name for figure: ") - # Make the plot - cice_evap_era_monthly(cice_file, month, colour_bounds, save, fig_name) - else: - break diff --git a/compare_nic_seasonal.py b/compare_nic_seasonal.py index ecd2a90..1707f74 100644 --- a/compare_nic_seasonal.py +++ b/compare_nic_seasonal.py @@ -1,7 +1,8 @@ from numpy import * -from netCDF4 import Dataset, num2date +from netCDF4 import Dataset from matplotlib.pyplot import * from rotate_vector_cice import * +from seasonal_avg_cice import * def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, fig_name=None): @@ -40,15 +41,6 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f lat_name = 'TLAT' id.close() - # Starting and ending months (1-based) for each season - start_month = [12, 3, 6, 9] - end_month = [2, 5, 8, 11] - # Starting and ending days of the month (1-based) for each season - # Assume no leap years, we'll fix this later if we need - start_day = [1, 1, 1, 1] - end_day = [28, 31, 31, 30] - # Number of days in each season (again, ignore leap years for now) - ndays_season = [90, 92, 92, 91] # Number of days in each month (this is just for Nic's output) ndays_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] # Season names for titles @@ -56,10 +48,13 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f # Degrees to radians conversion deg2rad = pi/180.0 - # Read CICE grid and time values + # Read the CICE grid id = Dataset(cice_file, 'r') cice_lon_tmp = id.variables[lon_name][:-15,:] cice_lat_tmp = id.variables[lat_name][:-15,:] + num_lon = id.variables['TLON'].shape[1] + num_lat = id.variables['TLAT'].shape[0] + id.close() # Wrap the periodic boundary by 1 cell cice_lon = ma.empty([size(cice_lon_tmp,0), size(cice_lon_tmp,1)+1]) cice_lat = ma.empty([size(cice_lat_tmp,0), size(cice_lat_tmp,1)+1]) @@ -67,175 +62,20 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f cice_lon[:,-1] = cice_lon_tmp[:,0] cice_lat[:,:-1] = cice_lat_tmp cice_lat[:,-1] = cice_lat_tmp[:,0] - time_id = id.variables['time'] - # Get the year, month, and day (all 1-based) for each output step - # These are 5-day averages marked with the next day's date. - cice_time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - - # Loop backwards through time indices to find the last one we care about - # (which contains 30 November in its averaging period) - end_t = -1 # Missing value flag - for t in range(size(cice_time)-1, -1, -1): - if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+5): - end_t = t - break - # Make sure we actually found it - if end_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' - return - - # Continue looping backwards to find the first time index we care about - # (which contains 1 December the previous year in its averaging period) - start_t = -1 # Missing value flag - for t in range(end_t-60, -1, -1): - if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0]+1, start_day[0]+6): - start_t = t - break - # Make sure we actually found it - if start_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' - return - - # Check for leap years - leap_year = False - if mod(cice_time[end_t].year, 4) == 0: - # Years divisible by 4 are leap years - leap_year = True - if mod(cice_time[end_t].year, 100) == 0: - # Unless they're also divisible by 100, in which case they aren't - # leap years - leap_year = False - if mod(cice_time[end_t].year, 400) == 0: - # Unless they're also divisible by 400, in which case they are - # leap years after all - leap_year = True - if leap_year: - # Update last day in February - end_day[0] += 1 - ndays_season[0] += 1 - # Don't update ndays_month because Nic's setup has no leap years - - # Initialise seasonal averages of CICE output - cice_data_tmp = ma.empty([4, size(cice_lon_tmp,0), size(cice_lon_tmp,1)]) - cice_data_tmp[:,:,:] = 0.0 - # Process one season at a time - for season in range(4): - season_days = 0 # Number of days in season; this will be incremented - next_season = mod(season+1, 4) - - # Find starting timestep - start_t_season = -1 - for t in range(start_t, end_t+1): - if cice_time[t].month == start_month[season] and cice_time[t].day in range(start_day[season]+1, start_day[season]+6): - start_t_season = t - break - # Make sure we actually found it - if start_t_season == -1: - print 'Error: could not find starting timestep for season ' + season_names[season] - return - - # Find ending timestep - end_t_season = -1 - for t in range(start_t_season+1, end_t+1): - if cice_time[t].month == start_month[next_season] and cice_time[t].day in range(start_day[next_season], start_day[next_season]+5): - end_t_season = t - break - # Make sure we actually found it - if end_t_season == -1: - print 'Error: could not find ending timestep for season ' + season_names[season] - return - - # Figure out how many of the 5 days averaged in start_t_season are - # actually within this season - if cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 5: - # Starting day is in position 1 of 5; we care about all of them - start_days = 5 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 4: - # Starting day is in position 2 of 5; we care about the last 4 - start_days = 4 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]+ 3: - # Starting day is in position 3 of 5; we care about the last 3 - start_days = 3 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 2: - # Starting day is in position 4 of 5; we care about the last 2 - start_days = 2 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 1: - # Starting day is in position 5 of 5; we care about the last 1 - start_days = 1 - else: - print 'Error for season ' + season_names[season] + ': starting index is month ' + str(cice_time[start_t_season].month) + ', day ' + str(cice_time[start_t_season].day) - return - # Start accumulating data weighted by days - if rotate: - this_cmp = id.variables[var_name][start_t_season,:-15,:] - other_cmp = id.variables[other_name][start_t_season,:-15,:] - if cmp_flag == 'x': - data_tmp, other_tmp = rotate_vector_cice(this_cmp, other_cmp, angle) - elif cmp_flag == 'y': - other_tmp, data_tmp = rotate_vector_cice(other_cmp, this_cmp, angle) - cice_data_tmp[season,:,:] += data_tmp*start_days - else: - cice_data_tmp[season,:,:] += id.variables[var_name][start_t_season,:-15,:]*start_days - season_days += start_days - - # Between start_t_season and end_t_season, we want all the days - for t in range(start_t_season+1, end_t_season): - if rotate: - this_cmp = id.variables[var_name][t,:-15,:] - other_cmp = id.variables[other_name][t,:-15,:] - if cmp_flag == 'x': - data_tmp, other_tmp = rotate_vector_cice(this_cmp, other_cmp, angle) - elif cmp_flag == 'y': - other_tmp, data_tmp = rotate_vector_cice(other_cmp, this_cmp, angle) - cice_data_tmp[season,:,:] += data_tmp*5 - else: - cice_data_tmp[season,:,:] += id.variables[var_name][t,:-15,:]*5 - season_days += 5 - - # Figure out how many of the 5 days averaged in end_t_season are - # actually within this season - if cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 4: - # Ending day is in position 1 of 5; we care about the first 1 - end_days = 1 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 3: - # Ending day is in position 2 of 5; we care about the first 2 - end_days = 2 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 2: - # Ending day is in position 3 of 5; we care about the first 3 - end_days = 3 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 1: - # Ending day is in position 4 of 5; we care about the first 4 - end_days = 4 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season]: - # Ending day is in position 5 of 5; we care about all 5 - end_days = 5 - else: - print 'Error for season ' + season_names[season] + ': ending index is month ' + str(cice_time[end_t_season].month) + ', day ' + str(cice_time[end_t_season].day) - return - - if rotate: - this_cmp = id.variables[var_name][end_t_season,:-15,:] - other_cmp = id.variables[other_name][end_t_season,:-15,:] - if cmp_flag == 'x': - data_tmp, other_tmp = rotate_vector_cice(this_cmp, other_cmp, angle) - elif cmp_flag == 'y': - other_tmp, data_tmp = rotate_vector_cice(other_cmp, this_cmp, angle) - cice_data_tmp[season,:,:] += data_tmp*end_days - else: - cice_data_tmp[season,:,:] += id.variables[var_name][end_t_season,:-15,:]*end_days - season_days += end_days - - # Check that we got the correct number of days - if season_days != ndays_season[season]: - print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) - return - - # Finished accumulating data, now convert from sum to average - cice_data_tmp[season,:,:] /= season_days + # Get seasonal averages of CICE data + if rotate: + this_cmp = seasonal_avg_cice(cice_file, var_name, [num_lat, num_lon]) + other_cmp = seasonal_avg_cice(cice_file, other_name, [num_lat, num_lon]) + if cmp_flag == 'x': + cice_data_tmp, other_tmp = rotate_vector_cice(this_cmp, other_cmp, angle) + elif cmp_flag == 'y': + other_tmp, cice_data_tmp = rotate_vector_cice(other_cmp, this_cmp, angle) + else: + cice_data_tmp = seasonal_avg_cice(cice_file, var_name, [num_lat, num_lon]) - # Finished reading all CICE data - id.close() + # Chop off northern boundary + cice_data_tmp = cice_data_tmp[:,:-15,:] # Wrap the periodic boundary cice_data = ma.empty([size(cice_data_tmp,0), size(cice_data_tmp,1), size(cice_data_tmp,2)+1]) cice_data[:,:,:-1] = cice_data_tmp diff --git a/interp_era2roms.py b/interp_era2roms.py new file mode 100644 index 0000000..1e508b4 --- /dev/null +++ b/interp_era2roms.py @@ -0,0 +1,54 @@ +from numpy import * +from scipy.interpolate import LinearNDInterpolator, RectBivariateSpline + +# Given an array on the ERA-Interim grid, interpolate any missing values, and +# then interpolate to the ROMS grid. +# Input: +# A = array of size nxm containing values on the ERA-Interim grid (dimension +# longitude x latitude) +# lon_era = array of length n containing ERA-Interim longitude values +# lat_era = array of length m containing ERA-Interim latitude values +# lon_roms = array of size pxq containing ROMS longitude values +# lat_roms = array of size pxq containing ROMS latitude values +# Output: +# B = array of size pxq containing values on the ROMS grid (dimension latitude x +# longitude) +def interp_era2roms (A, lon_era, lat_era, lon_roms, lat_roms): + + # Save the sizes of ROMS axes + num_lon = size(lon_roms, 1) + num_lat = size(lon_roms, 0) + + # Missing values are something <<0, but these can change when the offset + # and scale factor attributes are automatically applied. Either way, + # missing values will be the minimum values in A. + flag = amin(A) + + # Interpolate missing values + # I got this bit of code from Stack Exchange + # It seems to work, not exactly sure how + valid_mask = A > flag + coords = array(nonzero(valid_mask)).T + values = A[valid_mask] + fill_function = LinearNDInterpolator(coords, values) + Afill = fill_function(list(ndindex(A.shape))).reshape(A.shape) + # Fill any still-missing values with the mean + Afill[isnan(Afill)] = nanmean(Afill) + + # Now interpolate from ERA-Interim grid to ROMS grid + # First build a function to approximate A with 2D splines + # Note that latitude has to be multiplied by -1 so both axes are strictly + # ascending + interp_function = RectBivariateSpline(lon_era, -lat_era, Afill) + B = zeros(shape(lon_roms)) + # Call this function for each grid point individually - if you try to do + # it all at once it throws a MemoryError + for i in range(num_lon): + for j in range(num_lat): + B[j,i] = interp_function(lon_roms[j,i], -lat_roms[j,i]) + + # Enforce periodic boundary + B[:,0] = B[:,-2] + B[:,-1] = B[:,1] + + return B diff --git a/interp_lon_roms.py b/interp_lon_roms.py new file mode 100644 index 0000000..28f7cad --- /dev/null +++ b/interp_lon_roms.py @@ -0,0 +1,106 @@ +from numpy import * + +# Linearly interpolate data, z, and latitude to the specified longitude. +# Input: +# data_3d = array of data, dimension depth x lat x lon +# z_3d = array of depth values (negative, in metres), dimension depth x lat x lon +# lat_2d = array of latitudevalues, dimension lat x lon +# lon_2d = array of longitude values, dimension lat x lon (between -180 and 180) +# lon0 = longitude to interpolate to (between -180 and 180) +# Output: +# data = array of data interpolated to lon0, dimension depth x lat +# z = array of depth values interpolated to lon0, dimension depth x lat +# lat = array of latitude values interpolated to lon0, dimension depth x lat +def interp_lon_roms (data_3d, z_3d, lat_2d, lon_2d, lon0): + + # Save dimensions + num_depth = size(data_3d, 0) + num_lat = size(data_3d, 1) + num_lon = size(data_3d, 2) + # Set up output arrays + data = ma.empty([num_depth, num_lat]) + z = ma.empty([num_depth, num_lat]) + lat = ma.empty([num_depth, num_lat]) + + # Loop over latitudes; can't find a cleaner way to do this + for j in range(num_lat): + # Extract the longitude values of this slice + lon_tmp = lon_2d[j,:] + # Get indices and coefficients for interpolation + ie, iw, coeffe, coeffw = interp_lon_helper(lon_tmp, lon0) + data[:,j] = coeffe*data_3d[:,j,ie] + coeffw*data_3d[:,j,iw] + z[:,j] = coeffe*z_3d[:,j,ie] + coeffw*z_3d[:,j,iw] + lat[:,j] = coeffe*lat_2d[j,ie] + coeffw*lat_2d[j,iw] + + return data, z, lat + + +# Calculate indices and coefficients for linear interpolation of longitude. +# This takes care of all the mod 360 nonsense. +# Input: +# lon = 1D array of longitude values (straight out of ROMS i.e. between slightly < 0 and slightly > 360) +# lon0 = longitude to interpolate to (between 0 and 360) +# Output: +# ie, iw, coeffe, coeffw = integers (ie and iw) and coefficients +# (coeffe and coeffw) such that +# coeffe*lon[ie] + coeffw*lon[iw] = lon0, +# which will also hold for any variable on this +# longitude grid. ie is the index of the nearest +# point to the east of lon0; iw the nearest point +# to the west. +def interp_lon_helper (lon, lon0): + + if lon0 < amin(lon) or lon0 > amax(lon): + # Special case: lon0 on periodic boundary + # Be careful with mod 360 here + + # Find the periodic boundary + dlon = lon[1:] - lon[0:-1] + bdry = argmax(abs(dlon)) + if dlon[bdry] < -300: + # Jumps from almost 360 to just over 0 + iw = bdry + ie = bdry + 1 + else: + # Periodic boundary lines up with the array boundary + iw = size(lon) - 1 + ie = 0 + # Calculate difference between lon0 and lon[iw], mod 360 if necessary + dlon_num = lon0 - lon[iw] + if dlon_num < -300: + dlon_num += 360 + # Calculate difference between lon[ie] and lon[iw], mod 360 + dlon_den = lon[ie] - lon[iw] + 360 + + else: + # General case + + # Add or subtract 360 from longitude values which wrap around + # so that longitude increases monotonically from west to east + i = arange(1, size(lon)+1) + index1 = nonzero((i > 1200)*(lon < 100)) + lon[index1] = lon[index1] + 360 + index2 = nonzero((i < 200)*(lon > 300)) + lon[index2] = lon[index2] - 360 + + # Take mod 360 of lon0 if necessary + if all(lon < lon0): + lon0 -= 360 + if all(lon > lon0): + lon0 += 360 + + # Find the first index eastward of lon0 + ie = nonzero(lon > lon0)[0][0] + # The index before it will be the last index westward of lon0 + iw = ie - 1 + + dlon_num = lon0 - lon[iw] + dlon_den = lon[ie] - lon[iw] + + if dlon_num > 5 or dlon_den > 5: + print 'interp_lon_helper: Problem at periodic boundary' + return + coeff1 = dlon_num/dlon_den + coeff2 = 1 - coeff1 + + return ie, iw, coeff1, coeff2 diff --git a/interp_lon_sose.py b/interp_lon_sose.py new file mode 100644 index 0000000..e5db36e --- /dev/null +++ b/interp_lon_sose.py @@ -0,0 +1,37 @@ +from numpy import * + +# Linearly interpolate SOSE data to the specified longitude. +# Input: +# data_3d = arary of data, dimension depth x lat x lon +# lon = 1D array of longitude values (between 0 and 360) +# lon0 = longitude to interpolate to (between 0 and 360) +# Output: +# data = array of data interpolated to lon0, dimension depth x lat +def interp_lon_sose (data_3d, lon, lon0): + + if lon0 < lon[0] or lon0 > lon[-1]: + # Special case: lon0 on periodic boundary + # Be careful with mod 360 here + iw = size(lon)-1 + ie = 0 + # Calculate difference between lon0 and lon[iw], mod 360 if necessary + dlon_num = lon0 - lon[iw] + if dlon_num < -300: + dlon_num += 360 + # Calculate difference between lon[ie] and lon[iw], mod 360 + dlon_den = lon[ie] - lon[iw] + 360 + else: + # General case + # Find the first index eastwards of lon0 + ie = nonzero(lon > lon0)[0][0] + # The index before it will be the last index westward of lon0 + iw = ie - 1 + dlon_num = lon0 - lon[iw] + dlon_den = lon[ie] - lon[iw] + # Coefficients for interpolation + coeff1 = dlon_num/dlon_den + coeff2 = 1 - coeff1 + + # Interpolate + data = coeff1*data_3d[:,:,ie] + coeff2*data_3d[:,:,iw] + return data diff --git a/monthly_avg_cice.py b/monthly_avg_cice.py new file mode 100644 index 0000000..e743d09 --- /dev/null +++ b/monthly_avg_cice.py @@ -0,0 +1,117 @@ +from netCDF4 import Dataset, num2date +from numpy import * + +def monthly_avg_cice (file_path, var, shape, month): + + # Starting and ending days in each month + # Assume no leap years, we'll fix this later if we need + start_day = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] + end_day = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] + month_name = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] + + # Set up array to hold monthly data + monthly_data = ma.empty(shape) + monthly_data[:] = 0.0 + + # Read the time values + id = Dataset(file_path, 'r') + time_id = id.variables['time'] + # Get the year, month, and day (all 1-indexed) for each output step + # These are 5-day averages marked with the next day's date. + cice_time = num2date(time_id[:], units=time_id.units, calendar='standard') + + next_month = mod(month+1,12) + + end_t = -1 + for t in range(size(cice_time)-1, -1, -1): + if cice_time[t].month-1 == next_month and cice_time[t].day in range(start_day[next_month], start_day[next_month]+5): + end_t = t + break + if end_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete ' + month_name[month] + return + + start_t = -1 + for t in range(end_t, -1, -1): + if cice_time[t].month-1 == month and cice_time[t].day in range(start_day[month]+1, start_day[month]+6): + start_t = t + break + if start_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete ' + month_name[month] + return + + leap_year = False + cice_year = cice_time[end_t].year + if month == 11: + cice_year = cice_time[start_t].year + if mod(cice_year, 4) == 0: + leap_year = True + if mod(cice_year, 100) == 0: + leap_year = False + if mod(cice_year, 400) == 0: + leap_year = True + if leap_year: + end_day[1] = 29 + + num_days = 0 + + # Figure out how many of the 5 days averaged in start_t are actually within + # this month + if cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+5: + # Starting day is in position 1 of 5; we care about all of them + start_days = 5 + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+4: + # Starting day is in position 2 of 5; we care about the last 4 + start_days = 4 + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+3: + # Starting day is in position 3 of 5; we care about the last 3 + start_days = 3 + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+2: + # Starting day is in position 4 of 5; we care about the last 2 + start_days = 2 + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+1: + # Starting day is in position 5 of 5; we care about the last 1 + start_days = 1 + else: + print 'Error: starting index is month ' + str(cice_time[start_t].month) + ', day ' + str(cice_time[start_t].day) + return + + monthly_data += id.variables[var][start_t,:]*start_days + num_days += start_days + + for t in range(start_t+1, end_t): + monthly_data += id.variables[var][t,:]*5 + num_days +=5 + + # Figure out how many of the 5 days averaged in end_t are actually within + # this month + if cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+4: + # Ending day is in position 1 of 5; we care about the first 1 + end_days = 1 + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+3: + # Ending day is in position 2 of 5; we care about the first 2 + end_days = 2 + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+2: + # Ending day is in position 3 of 5; we care about the first 3 + end_days = 3 + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+1: + # Ending day is in position 4 of 5; we care about the first 4 + end_days = 4 + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]: + # Ending day is in position 5 of 5; we care about all 5 + end_days = 5 + else: + print 'Error: ending index is month ' + str(cice_time[end_t].month) + ', day ' + str(cice_time[end_t].day) + return + + monthly_data += id.variables[var][end_t,:]*end_days + num_days += end_days + + if num_days != end_day[month]: + print 'Error: found ' + str(num_days) + ' days instead of ' + str(end_day[month]) + return + + id.close() + monthly_data /= num_days + + return monthly_data diff --git a/monthly_avg_roms.py b/monthly_avg_roms.py new file mode 100644 index 0000000..7614c32 --- /dev/null +++ b/monthly_avg_roms.py @@ -0,0 +1,112 @@ +from netCDF4 import Dataset, num2date +from numpy import * + +def monthly_avg_roms (file_path, var, shape, month): + + # Starting and ending days in each month + # Assume no leap years, we'll fix this later if we need + start_day = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] + end_day = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] + month_name = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] + + # Set up array to hold monthly data + monthly_data = ma.empty(shape) + monthly_data[:] = 0.0 + + # Read the time values + id = Dataset(file_path, 'r') + time_id = id.variables['ocean_time'] + # Get the year, month, and day (all 1-based) for each output step + # These are 5-day averages marked with the middle day's date + time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) + + next_month = mod(month+1,12) + prev_month = mod(month-1,12) + + end_t = -1 + for t in range(size(time)-1, -1, -1): + if time[t].month-1 == month and time[t].day in range(end_day[month]-2, end_day[month]+1): + end_t = t + break + if time[t].month-1 == next_month and time[t].day in range(start_day[next_month], start_day[next_month]+2): + end_t = t + break + if end_t == -1: + print 'Error: ' + file_path + ' does not contain a complete ' + month_name[month] + return + + start_t = -1 + for t in range(end_t, -1, -1): + if time[t].month-1 == prev_month and time[t].day in range(end_day[prev_month]-1, end_day[prev_month]+1): + start_t = t + break + if time[t].month-1 == month and time[t].day in range(start_day[month], start_day[month]+3): + start_t = t + break + if start_t == -1: + print 'Error: ' + file_path + ' does not contain a complete ' + month_name[month] + return + + leap_year = False + roms_year = time[end_t].year + if month == 11: + roms_year = time[start_t].year + if mod(roms_year, 4) == 0: + leap_year = True + if mod(roms_year, 100) == 0: + leap_year = False + if mod(roms_year, 400) == 0: + leap_year = True + if leap_year: + end_day[1] = 29 + + num_days = 0 + + if time[start_t].month-1 == month and time[start_t].day == start_day[month]+2: + start_days = 5 + elif time[start_t].month-1 == month and time[start_t].day == start_day[month]+1: + start_days = 4 + elif time[start_t].month-1 == month and time[start_t].day == start_day[month]: + start_days = 3 + elif time[start_t].month-1 == prev_month and time[start_t].day == end_day[prev_month]: + start_days = 2 + elif time[start_t].month-1 == prev_month and time[start_t].day == end_day[prev_month]-1: + start_days = 1 + else: + print 'Error: starting index is month ' + str(time[start_t].month) + ', day ' + str(time[start_t].day) + return + + monthly_data += id.variables[var][start_t,:]*start_days + num_days += start_days + + for t in range(start_t+1, end_t): + monthly_data += id.variables[var][t,:]*5 + num_days += 5 + + if time[end_t].month-1 == next_month and time[end_t].day == start_day[next_month]+1: + end_days = 1 + elif time[end_t].month-1 == next_month and time[end_t].day == start_day[next_month]: + end_days = 2 + elif time[end_t].month-1 == month and time[end_t].day == end_day[month]: + end_days = 3 + elif time[end_t].month-1 == month and time[end_t].day == end_day[month]-1: + end_days = 4 + elif time[end_t].month-1 == month and time[end_t].day == end_day[month]-2: + end_days = 5 + else: + print 'Error: ending index is month ' + str(time[end_t].month) + ', day ' + str(time[end_t].day) + return + + monthly_data += id.variables[var][end_t,:]*end_days + num_days += end_days + + if num_days != end_day[month]: + print 'Error: found ' + str(num_days) + ' days instead of ' + str(end_day[month]) + return + + id.close() + monthly_data /= num_days + + return monthly_data + + diff --git a/nsidc_aice_monthly.py b/nsidc_aice_monthly.py index 8e6cbbd..24ee2e2 100644 --- a/nsidc_aice_monthly.py +++ b/nsidc_aice_monthly.py @@ -1,6 +1,7 @@ from numpy import * -from netCDF4 import Dataset, num2date +from netCDF4 import Dataset from matplotlib.pyplot import * +from monthly_avg_cice import * # Make a figure comparing sea ice concentration from NSIDC (1995 data) and # CICE (latest year of spinup under repeated 1995 forcing), for the given @@ -19,19 +20,18 @@ def nsidc_aice_monthly (cice_file, month, save=False, fig_name=None): nsidc_head_0 = nsidc_head + '_f11_' nsidc_head_1 = nsidc_head + '_f13_' nsidc_tail = '_v02r00.nc' - # Starting and ending days in each month - # Assume no leap years, we'll fix this later if we need - start_day = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] - end_day = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] # Name of each month, for the title month_name = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] # Degrees to radians conversion deg2rad = pi/180.0 - # Read CICE grid and time values + # Read CICE grid id = Dataset(cice_file, 'r') cice_lon_tmp = id.variables['TLON'][:-15,:] cice_lat_tmp = id.variables['TLAT'][:-15,:] + num_lon = id.variables['TLON'].shape[1] + num_lat = id.variables['TLAT'].shape[0] + id.close() # Wrap the periodic boundary by 1 cell cice_lon = ma.empty([size(cice_lon_tmp,0), size(cice_lon_tmp,1)+1]) cice_lat = ma.empty([size(cice_lat_tmp,0), size(cice_lat_tmp,1)+1]) @@ -39,125 +39,11 @@ def nsidc_aice_monthly (cice_file, month, save=False, fig_name=None): cice_lon[:,-1] = cice_lon_tmp[:,0] cice_lat[:,:-1] = cice_lat_tmp cice_lat[:,-1] = cice_lat_tmp[:,0] - time_id = id.variables['time'] - # Get the year, month, and day (all 1-based) for each output step - # These are 5-day averages marked with the next day's date. - cice_time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - - # Loop backwards through time indices to find the last one we care about - # (which contains the last day of the month in its averaging period) - end_t = -1 # Missing value flag - for t in range(size(cice_time)-1, -1, -1): - # Note that cice_time[t].month is converted to 0-indexed - if cice_time[t].month-1 == next_month and cice_time[t].day in range(start_day[next_month], start_day[next_month]+5): - end_t = t - break - # Make sure we actually found it - if end_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete ' + month_name[month] - return - - # Continue looping backwards to find the first time index we care about - start_t = -1 # Missing value flag - for t in range(end_t, -1, -1): - if cice_time[t].month-1 == month and cice_time[t].day in range(start_day[month]+1, start_day[month]+6): - start_t = t - break - # Make sure we actually found it - if start_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete ' + month_name[month] - return - - # Check if this is a leap year - leap_year = False - cice_year = cice_time[end_t].year - if month == 11: - cice_year = cice_time[start_t].year - if mod(cice_year, 4) == 0: - # Years divisible by 4 are leap years - leap_year = True - if mod(cice_year, 100) == 0: - # Unless they're also divisible by 100, in which case they aren't - # leap years - leap_year = False - if mod(cice_year, 400) == 0: - # Unless they're also divisible by 400, in which case they are - # leap years after all - leap_year = True - if leap_year: - # Update last day in February - end_day[1] = 29 - - # Calculate monthly average of CICE output - cice_data_tmp = ma.empty(shape(cice_lon_tmp)) - cice_data_tmp[:,:] = 0.0 - num_days = 0 - - # Figure out how many of the 5 days averaged in start_t are actually within - # this month - if cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month] + 5: - # Starting day is in position 1 of 5; we care about all of them - start_days = 5 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month] + 4: - # Starting day is in position 2 of 5; we care about the last 4 - start_days = 4 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+ 3: - # Starting day is in position 3 of 5; we care about the last 3 - start_days = 3 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month] + 2: - # Starting day is in position 4 of 5; we care about the last 2 - start_days = 2 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month] + 1: - # Starting day is in position 5 of 5; we care about the last 1 - start_days = 1 - else: - print 'Error: starting index is month ' + str(cice_time[start_t].month) + ', day ' + str(cice_time[start_t].day) - return - - # Start accumulating data weighted by days - cice_data_tmp[:,:] += id.variables['aice'][start_t,:-15,:]*start_days - num_days += start_days - - # Beween start_t and end_t, we want all the days - for t in range(start_t+1, end_t): - cice_data_tmp[:,:] += id.variables['aice'][t,:-15,:]*5 - num_days += 5 - - # Figure out how many of the 5 days averaged in end_t are actually within - # this month - next_month = mod(month+1, 12) - if cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month] + 4: - # Ending day is in position 1 of 5; we care about the first 1 - end_days = 1 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month] + 3: - # Ending day is in position 2 of 5; we care about the first 2 - end_days = 2 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month] + 2: - # Ending day is in position 3 of 5; we care about the first 3 - end_days = 3 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month] + 1: - # Ending day is in position 4 of 5; we care about the first 4 - end_days = 4 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]: - # Ending day is in position 5 of 5; we care about all 5 - end_days = 5 - else: - print 'Error: ending index is month ' + str(cice_time[end_t].month) + ', day ' + str(cice_time[end_t].day) - return - - cice_data_tmp[:,:] += id.variables['aice'][end_t,:-15,:]*end_days - num_days += end_days - - # Check that we got the correct number of days - if num_days != end_day[month]: - print 'Error: found ' + str(num_days) + ' days instead of ' + str(end_day[month]) - return - - # Finished accumulating data - id.close() - # Convert from sum to average - cice_data_tmp[:,:] /= num_days + # Get monthly average of sea ice area + cice_data_tmp = monthly_avg_cice(cice_file, 'aice', [num_lat, num_lon], month) + # Chop off the northern boundary + cice_data_tmp = cice_data_tmp[:-15,:] # Wrap the periodic boundary cice_data = ma.empty([size(cice_data_tmp,0), size(cice_data_tmp,1)+1]) cice_data[:,:-1] = cice_data_tmp diff --git a/nsidc_aice_seasonal.py b/nsidc_aice_seasonal.py index 59bc4b9..81c7410 100644 --- a/nsidc_aice_seasonal.py +++ b/nsidc_aice_seasonal.py @@ -1,6 +1,7 @@ from numpy import * -from netCDF4 import Dataset, num2date +from netCDF4 import Dataset from matplotlib.pyplot import * +from seasonal_avg_cice import * # Make a figure comparing sea ice concentration from NSIDC (1995 data) and CICE # (latest year of spinup under repeated 1995 forcing) for each season. @@ -17,16 +18,6 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): nsidc_head_0 = nsidc_head + '_f11_' nsidc_head_1 = nsidc_head + '_f13_' nsidc_tail = '_v02r00.nc' - - # Starting and ending months (1-based) for each season - start_month = [12, 3, 6, 9] - end_month = [2, 5, 8, 11] - # Starting and ending days of the month (1-based) for each season - # Assume no leap years, we'll fix this later if we need - start_day = [1, 1, 1, 1] - end_day = [28, 31, 31, 30] - # Number of days in each season (again, ignore leap years for now) - ndays_season = [90, 92, 92, 91] # Number of days in each month (this is just for NSIDC) ndays_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] # Season names for titles @@ -34,10 +25,13 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): # Degrees to radians conversion deg2rad = pi/180.0 - # Read CICE grid and time values + # Read the CICE grid id = Dataset(cice_file, 'r') cice_lon_tmp = id.variables['TLON'][:-15,:] cice_lat_tmp = id.variables['TLAT'][:-15,:] + num_lon = id.variables['TLON'].shape[1] + num_lat = id.variables['TLAT'].shape[0] + id.close() # Wrap the periodic boundary by 1 cell cice_lon = ma.empty([size(cice_lon_tmp,0), size(cice_lon_tmp,1)+1]) cice_lat = ma.empty([size(cice_lat_tmp,0), size(cice_lat_tmp,1)+1]) @@ -45,153 +39,15 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): cice_lon[:,-1] = cice_lon_tmp[:,0] cice_lat[:,:-1] = cice_lat_tmp cice_lat[:,-1] = cice_lat_tmp[:,0] - time_id = id.variables['time'] - # Get the year, month, and day (all 1-based) for each output step - # These are 5-day averages marked with the next day's date. - cice_time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - - # Loop backwards through time indices to find the last one we care about - # (which contains 30 November in its averaging period) - end_t = -1 # Missing value flag - for t in range(size(cice_time)-1, -1, -1): - if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+5): - end_t = t - break - # Make sure we actually found it - if end_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' - return - - # Continue looping backwards to find the first time index we care about - # (which contains 1 December the previous year in its averaging period) - start_t = -1 # Missing value flag - for t in range(end_t-60, -1, -1): - if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0]+1, start_day[0]+6): - start_t = t - break - # Make sure we actually found it - if start_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' - return - - # Check for leap years - leap_year = False - if mod(cice_time[end_t].year, 4) == 0: - # Years divisible by 4 are leap years - leap_year = True - if mod(cice_time[end_t].year, 100) == 0: - # Unless they're also divisible by 100, in which case they aren't - # leap years - leap_year = False - if mod(cice_time[end_t].year, 400) == 0: - # Unless they're also divisible by 400, in which case they are - # leap years after all - leap_year = True - if leap_year: - # Update last day in February - end_day[0] += 1 - ndays_season[0] += 1 - ndays_month[1] += 1 - - # Initialise seasonal averages of CICE output - cice_data_tmp = ma.empty([4, size(cice_lon_tmp,0), size(cice_lon_tmp,1)]) - cice_data_tmp[:,:,:] = 0.0 - # Process one season at a time - for season in range(4): - season_days = 0 # Number of days in season; this will be incremented - next_season = mod(season+1, 4) - - # Find starting timestep - start_t_season = -1 - for t in range(start_t, end_t+1): - if cice_time[t].month == start_month[season] and cice_time[t].day in range(start_day[season]+1, start_day[season]+6): - start_t_season = t - break - # Make sure we actually found it - if start_t_season == -1: - print 'Error: could not find starting timestep for season ' + season_names[season] - return - - # Find ending timestep - end_t_season = -1 - for t in range(start_t_season+1, end_t+1): - if cice_time[t].month == start_month[next_season] and cice_time[t].day in range(start_day[next_season], start_day[next_season]+5): - end_t_season = t - break - # Make sure we actually found it - if end_t_season == -1: - print 'Error: could not find ending timestep for season ' + season_names[season] - return - - # Figure out how many of the 5 days averaged in start_t_season are - # actually within this season - if cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 5: - # Starting day is in position 1 of 5; we care about all of them - start_days = 5 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 4: - # Starting day is in position 2 of 5; we care about the last 4 - start_days = 4 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]+ 3: - # Starting day is in position 3 of 5; we care about the last 3 - start_days = 3 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 2: - # Starting day is in position 4 of 5; we care about the last 2 - start_days = 2 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 1: - # Starting day is in position 5 of 5; we care about the last 1 - start_days = 1 - else: - print 'Error for season ' + season_names[season] + ': starting index is month ' + str(cice_time[start_t_season].month) + ', day ' + str(cice_time[start_t_season].day) - return - - # Start accumulating data weighted by days - cice_data_tmp[season,:,:] += id.variables['aice'][start_t_season,:-15,:]*start_days - season_days += start_days - - # Between start_t_season and end_t_season, we want all the days - for t in range(start_t_season+1, end_t_season): - cice_data_tmp[season,:,:] += id.variables['aice'][t,:-15,:]*5 - season_days += 5 - - # Figure out how many of the 5 days averaged in end_t_season are - # actually within this season - if cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 4: - # Ending day is in position 1 of 5; we care about the first 1 - end_days = 1 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 3: - # Ending day is in position 2 of 5; we care about the first 2 - end_days = 2 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 2: - # Ending day is in position 3 of 5; we care about the first 3 - end_days = 3 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 1: - # Ending day is in position 4 of 5; we care about the first 4 - end_days = 4 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season]: - # Ending day is in position 5 of 5; we care about all 5 - end_days = 5 - else: - print 'Error for season ' + season_names[season] + ': ending index is month ' + str(cice_time[end_t_season].month) + ', day ' + str(cice_time[end_t_season].day) - return - - cice_data_tmp[season,:,:] += id.variables['aice'][end_t_season,:-15,:]*end_days - season_days += end_days - - # Check that we got the correct number of days - if season_days != ndays_season[season]: - print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) - return - - # Finished accumulating data, now convert from sum to average - cice_data_tmp[season,:,:] /= season_days - - # Finished reading all CICE data - id.close() + # Read seasonally averaged fields + cice_data_tmp = seasonal_avg_cice(cice_file, 'aice', [num_lat, num_lon]) + # Chop off northern boundary + cice_data_tmp = cice_data_tmp[:,:-15,:] # Wrap the periodic boundary cice_data = ma.empty([size(cice_data_tmp,0), size(cice_data_tmp,1), size(cice_data_tmp,2)+1]) cice_data[:,:,:-1] = cice_data_tmp - cice_data[:,:,-1] = cice_data_tmp[:,0] + cice_data[:,:,-1] = cice_data_tmp[:,:,0] # Read NSIDC grid from the January file id = Dataset(nsidc_head_0 + '199501' + nsidc_tail, 'r') diff --git a/romscice_atm_monthly.py b/romscice_atm_monthly.py index 38ffe42..f2823a3 100644 --- a/romscice_atm_monthly.py +++ b/romscice_atm_monthly.py @@ -1,6 +1,6 @@ from netCDF4 import Dataset from numpy import * -from scipy.interpolate import LinearNDInterpolator, RectBivariateSpline +from interp_era2roms import * # Convert two ERA-Interim files: # AN_yyyy_monthly_orig.nc: one year of monthly averaged measurements for @@ -173,56 +173,6 @@ def convert_file (year): oppt_fid.close() - -# Given an array on the ERA-Interim grid, interpolate any missing values, and -# then interpolate to the ROMS grid. -# Input: -# A = array of size nxm containing values on the ERA-Interim grid (dimension -# longitude x latitude) -# lon_era = array of length n containing ERA-Interim longitude values -# lat_era = array of length m containing ERA-Interim latitude values -# lon_roms = array of size pxq containing ROMS longitude values -# lat_roms = array of size pxq containing ROMS latitude values -# Output: -# B = array of size pxq containing values on the ROMS grid (dimension latitude x -# longitude) -def interp_era2roms (A, lon_era, lat_era, lon_roms, lat_roms): - - # Save the sizes of ROMS axes - num_lon = size(lon_roms, 1) - num_lat = size(lon_roms, 0) - - # Missing values are something <<0, but these can change when the offset - # and scale factor attributes are automatically applied. Either way, - # missing values will be the minimum values in A. - flag = amin(A) - - # Interpolate missing values - # I got this bit of code from Stack Exchange - # It seems to work, not exactly sure how - valid_mask = A > flag - coords = array(nonzero(valid_mask)).T - values = A[valid_mask] - fill_function = LinearNDInterpolator(coords, values) - Afill = fill_function(list(ndindex(A.shape))).reshape(A.shape) - # Fill any still-missing values with the mean - Afill[isnan(Afill)] = nanmean(Afill) - - # Now interpolate from ERA-Interim grid to ROMS grid - # First build a function to approximate A with 2D splines - # Note that latitude has to be multiplied by -1 so both axes are strictly - # ascending - interp_function = RectBivariateSpline(lon_era, -lat_era, Afill) - B = zeros(shape(lon_roms)) - # Call this function for each grid point individually - if you try to do - # it all at once it throws a MemoryError - for i in range(num_lon): - for j in range(num_lat): - B[j,i] = interp_function(lon_roms[j,i], -lat_roms[j,i]) - - return B - - # Command-line interface if __name__ == "__main__": diff --git a/romscice_atm_subdaily.py b/romscice_atm_subdaily.py index 2e78c51..43f7943 100644 --- a/romscice_atm_subdaily.py +++ b/romscice_atm_subdaily.py @@ -1,6 +1,6 @@ from netCDF4 import Dataset from numpy import * -from scipy.interpolate import LinearNDInterpolator, RectBivariateSpline +from interp_era2roms import * # Convert two ERA-Interim files: # AN_yyyy_subdaily_orig.nc: one year of 6-hour measurements for surface pressure @@ -229,64 +229,3 @@ def convert_file (year, count): log = open(logfile, 'a') log.write('Finished\n') log.close() - - -# Given an array on the ERA-Interim grid, interpolate any missing values, and -# then interpolate to the ROMS grid. -# Input: -# A = array of size nxm containing values on the ERA-Interim grid (dimension -# longitude x latitude) -# lon_era = array of length n containing ERA-Interim longitude values -# lat_era = array of length m containing ERA-Interim latitude values -# lon_roms = array of size pxq containing ROMS longitude values -# lat_roms = array of size pxq containing ROMS latitude values -# Output: -# B = array of size pxq containing values on the ROMS grid (dimension latitude x -# longitude) -def interp_era2roms (A, lon_era, lat_era, lon_roms, lat_roms): - - # Save the sizes of ROMS axes - num_lon = size(lon_roms, 1) - num_lat = size(lon_roms, 0) - - # Missing values are something <<0, but these can change when the offset - # and scale factor attributes are automatically applied. Either way, - # missing values will be the minimum values in A. - flag = amin(A) - - # Interpolate missing values - # I got this bit of code from Stack Exchange - # It seems to work, not exactly sure how - valid_mask = A > flag - coords = array(nonzero(valid_mask)).T - values = A[valid_mask] - fill_function = LinearNDInterpolator(coords, values) - Afill = fill_function(list(ndindex(A.shape))).reshape(A.shape) - # Fill any still-missing values with the mean - Afill[isnan(Afill)] = nanmean(Afill) - - # Now interpolate from ERA-Interim grid to ROMS grid - # First build a function to approximate A with 2D splines - # Note that latitude has to be multiplied by -1 so both axes are strictly - # ascending - interp_function = RectBivariateSpline(lon_era, -lat_era, Afill) - B = zeros(shape(lon_roms)) - # Call this function for each grid point individually - if you try to do - # it all at once it throws a MemoryError - for i in range(num_lon): - for j in range(num_lat): - B[j,i] = interp_function(lon_roms[j,i], -lat_roms[j,i]) - - # Enforce periodic boundary - B[:,0] = B[:,-2] - B[:,-1] = B[:,1] - - return B - - - - - - - - diff --git a/romscice_evap.py b/romscice_evap.py index fef0fac..e4346d4 100644 --- a/romscice_evap.py +++ b/romscice_evap.py @@ -1,6 +1,6 @@ from netCDF4 import Dataset from numpy import * -from scipy.interpolate import LinearNDInterpolator, RectBivariateSpline +from interp_era2roms import * def convert_file (year, count): @@ -81,57 +81,6 @@ def convert_file (year, count): log.close() -# Given an array on the ERA-Interim grid, interpolate any missing values, and -# then interpolate to the ROMS grid. -# Input: -# A = array of size nxm containing values on the ERA-Interim grid (dimension -# longitude x latitude) -# lon_era = array of length n containing ERA-Interim longitude values -# lat_era = array of length m containing ERA-Interim latitude values -# lon_roms = array of size pxq containing ROMS longitude values -# lat_roms = array of size pxq containing ROMS latitude values -# Output: -# B = array of size pxq containing values on the ROMS grid (dimension latitude x -# longitude) -def interp_era2roms (A, lon_era, lat_era, lon_roms, lat_roms): - - # Save the sizes of ROMS axes - num_lon = size(lon_roms, 1) - num_lat = size(lon_roms, 0) - - # Missing values are something <<0, but these can change when the offset - # and scale factor attributes are automatically applied. Either way, - # missing values will be the minimum values in A. - flag = amin(A) - - # Interpolate missing values - # I got this bit of code from Stack Exchange - # It seems to work, not exactly sure how - valid_mask = A > flag - coords = array(nonzero(valid_mask)).T - values = A[valid_mask] - fill_function = LinearNDInterpolator(coords, values) - Afill = fill_function(list(ndindex(A.shape))).reshape(A.shape) - # Fill any still-missing values with the mean - Afill[isnan(Afill)] = nanmean(Afill) - - # Now interpolate from ERA-Interim grid to ROMS grid - # First build a function to approximate A with 2D splines - # Note that latitude has to be multiplied by -1 so both axes are strictly - # ascending - interp_function = RectBivariateSpline(lon_era, -lat_era, Afill) - B = zeros(shape(lon_roms)) - # Call this function for each grid point individually - if you try to do - # it all at once it throws a MemoryError - for i in range(num_lon): - for j in range(num_lat): - B[j,i] = interp_function(lon_roms[j,i], -lat_roms[j,i]) - - # Enforce periodic boundary - B[:,0] = B[:,-2] - B[:,-1] = B[:,1] - - return B diff --git a/romscice_ini_ecco.py b/romscice_ini_ecco.py index 62b2d6e..df32d6c 100644 --- a/romscice_ini_ecco.py +++ b/romscice_ini_ecco.py @@ -3,7 +3,7 @@ # height to zero. Under the ice shelves, extrapolate temperature and salinity # from the ice shelf front. # NB: Users will likely need to edit paths to ECCO2 data! Scroll down below -# the interp_ecco2roms function to find where these filenames are defined. +# the interp_ecco2roms_ini function to find where these filenames are defined. # NB for raijin users: RegularGridInterpolator needs python/2.7.6 but the # default is 2.7.3. Before running this script, switch them as follows: @@ -91,9 +91,9 @@ def run (grid_file, theta_file, salt_file, output_file, Tcline, theta_s, theta_b # Regridding happens here... print 'Interpolating temperature' - temp = interp_ecco2roms(theta, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, -0.5) + temp = interp_ecco2roms_ini(theta, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, -0.5) print 'Interpolating salinity' - salt = interp_ecco2roms(salt, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, 34.5) + salt = interp_ecco2roms_ini(salt, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, 34.5) # Set initial velocities and sea surface height to zero u = zeros((N, num_lat, num_lon-1)) @@ -213,7 +213,7 @@ def run (grid_file, theta_file, salt_file, output_file, Tcline, theta_s, theta_b # Output: # B = array of size pxqxr containing values on the ROMS grid (dimension depth x # latitude x longitude) -def interp_ecco2roms (A, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, fill): +def interp_ecco2roms_ini (A, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, fill): # Radius of the Earth in m r = 6.371e6 diff --git a/romscice_nbc.py b/romscice_nbc.py index f8deb30..d7da630 100644 --- a/romscice_nbc.py +++ b/romscice_nbc.py @@ -8,7 +8,7 @@ # and velocity output, containing boundary conditions at the northern boundary. # Output 12 sets of monthly-averaged data, one for each month of the given year. # NB: Users will likely need to edit paths to ECCO2 data! Scroll down below -# the interp_ecco2roms function to find where these filenames are defined. +# the interp_ecco2roms_nbc function to find where these filenames are defined. # NB: This clamps u and ubar to zero at the northern boundary. # NB for raijin users: RegularGridInterpolator needs python/2.7.6 but the @@ -246,11 +246,11 @@ def convert_file (year): # Regridding happens here... print 'Interpolating temperature' - temp_interp = interp_ecco2roms(theta, lon_ecco, lat_ecco, depth_ecco, lon_rho, lat_rho, z_rho, mean(theta), True) + temp_interp = interp_ecco2roms_nbc(theta, lon_ecco, lat_ecco, depth_ecco, lon_rho, lat_rho, z_rho, mean(theta), True) print 'Interpolating salinity' - salt_interp = interp_ecco2roms(salt, lon_ecco, lat_ecco, depth_ecco, lon_rho, lat_rho, z_rho, mean(salt), True) + salt_interp = interp_ecco2roms_nbc(salt, lon_ecco, lat_ecco, depth_ecco, lon_rho, lat_rho, z_rho, mean(salt), True) print 'Interpolating v' - v_interp = interp_ecco2roms(vvel, lon_ecco, lat_ecco, depth_ecco, lon_v, lat_v, z_v, 0, False) + v_interp = interp_ecco2roms_nbc(vvel, lon_ecco, lat_ecco, depth_ecco, lon_v, lat_v, z_v, 0, False) # Calculate vertical average of v to get vbar # Be sure to treat land mask carefully so we don't divide by 0 @@ -305,7 +305,7 @@ def convert_file (year): # Output: # B = array of size pxr containing values on the ROMS grid (dimension depth x # longitude) -def interp_ecco2roms (A, lon_ecco, lat_ecco, depth_ecco, lon_roms, lat_roms, z_roms, fill_value, wrap): +def interp_ecco2roms_nbc (A, lon_ecco, lat_ecco, depth_ecco, lon_roms, lat_roms, z_roms, fill_value, wrap): # Fill the ECCO2 land mask with a constant value close to the mean of A. # This is less than ideal because it will skew the interpolation slightly diff --git a/seasonal_avg_cice.py b/seasonal_avg_cice.py new file mode 100644 index 0000000..58b1af5 --- /dev/null +++ b/seasonal_avg_cice.py @@ -0,0 +1,163 @@ +from netCDF4 import Dataset, num2date +from numpy import * + +def seasonal_avg_cice (file_path, var, shape): + + # Starting and ending months (1-based) for each season + start_month = [12, 3, 6, 9] + end_month = [2, 5, 8, 11] + # Starting and ending days of the month (1-based) for each season + # Assume no leap years, we'll fix this later if we need + start_day = [1, 1, 1, 1] + end_day = [28, 31, 31, 30] + # Number of days in each season (again, ignore leap years for now) + ndays_season = [90, 92, 92, 91] + season_names = ['DJF', 'MAM', 'JJA', 'SON'] + + # Set up array to hold seasonal data + new_shape = [4] + shape + seasonal_data = ma.empty(new_shape) + seasonal_data[:] = 0.0 + + # Read the time values + id = Dataset(file_path, 'r') + time_id = id.variables['time'] + # Get the year, month, and day (all 1-indexed) for each output step + # These are 5-day averages marked with the next day's date. + cice_time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) + + # Loop backwards through time indices to find the last one we care about + # (which contains 30 November in its averaging period) + end_t = -1 # Missing value flag + for t in range(size(cice_time)-1, -1, -1): + if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+5): + end_t = t + break + # Make sure we actually found it + if end_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' + return + + # Continue looping backwards to find the first time index we care about + # (which contains 1 December the previous year in its averaging period) + start_t = -1 # Missing value flag + for t in range(end_t-60, -1, -1): + if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0]+1, start_day[0]+6): + start_t = t + break + # Make sure we actually found it + if start_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' + return + + # Check for leap years + leap_year = False + cice_year = cice_time[end_t].year + if mod(cice_year, 4) == 0: + # Years divisible by 4 are leap years + leap_year = True + if mod(cice_year, 100) == 0: + # Unless they're also divisible by 100, in which case they aren't + # leap years + leap_year = False + if mod(cice_year, 400) == 0: + # Unless they're also divisible by 400, in which case they are + # leap years after all + leap_year = True + if leap_year: + # Update last day in February + end_day[0] += 1 + ndays_season[0] += 1 + + # Process one season at a time + for season in range(4): + season_days = 0 # Number of days in season; this will be incremented + next_season = mod(season+1, 4) + + # Find starting timestep + start_t_season = -1 + for t in range(start_t, end_t+1): + if cice_time[t].month == start_month[season] and cice_time[t].day in range(start_day[season]+1, start_day[season]+6): + start_t_season = t + break + # Make sure we actually found it + if start_t_season == -1: + print 'Error: could not find starting timestep for season ' + season_names[season] + return + + # Find ending timestep + end_t_season = -1 + for t in range(start_t_season+1, end_t+1): + if cice_time[t].month == start_month[next_season] and cice_time[t].day in range(start_day[next_season], start_day[next_season]+5): + end_t_season = t + break + # Make sure we actually found it + if end_t_season == -1: + print 'Error: could not find ending timestep for season ' + season_names[season] + return + + # Figure out how many of the 5 days averaged in start_t_season are + # actually within this season + if cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 5: + # Starting day is in position 1 of 5; we care about all of them + start_days = 5 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 4: + # Starting day is in position 2 of 5; we care about the last 4 + start_days = 4 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]+ 3: + # Starting day is in position 3 of 5; we care about the last 3 + start_days = 3 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 2: + # Starting day is in position 4 of 5; we care about the last 2 + start_days = 2 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 1: + # Starting day is in position 5 of 5; we care about the last 1 + start_days = 1 + else: + print 'Error for season ' + season_names[season] + ': starting index is month ' + str(cice_time[start_t_season].month) + ', day ' + str(cice_time[start_t_season].day) + return + + # Start accumulating data weighted by days + seasonal_data[season,:] += id.variables[var][start_t_season,:]*start_days + season_days += start_days + + # Between start_t_season and end_t_season, we want all the days + for t in range(start_t_season+1, end_t_season): + seasonal_data[season,:] += id.variables[var][t,:]*5 + season_days += 5 + + # Figure out how many of the 5 days averaged in end_t_season are + # actually within this season + if cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 4: + # Ending day is in position 1 of 5; we care about the first 1 + end_days = 1 + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 3: + # Ending day is in position 2 of 5; we care about the first 2 + end_days = 2 + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 2: + # Ending day is in position 3 of 5; we care about the first 3 + end_days = 3 + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 1: + # Ending day is in position 4 of 5; we care about the first 4 + end_days = 4 + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season]: + # Ending day is in position 5 of 5; we care about all 5 + end_days = 5 + else: + print 'Error for season ' + season_names[season] + ': ending index is month ' + str(cice_time[end_t_season].month) + ', day ' + str(cice_time[end_t_season].day) + return + + seasonal_data[season,:] += id.variables[var][end_t_season,:]*end_days + season_days += end_days + + # Check that we got the correct number of days + if season_days != ndays_season[season]: + print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) + return + + # Finished accumulating data, now convert from sum to average + seasonal_data[season,:,:] /= season_days + + id.close() + + return seasonal_data diff --git a/seasonal_avg_roms.py b/seasonal_avg_roms.py new file mode 100644 index 0000000..c3d5eca --- /dev/null +++ b/seasonal_avg_roms.py @@ -0,0 +1,180 @@ +from netCDF4 import Dataset, num2date +from numpy import * + +def seasonal_avg_roms (file_path, var, shape): + + # Starting and ending months (1-based) for each season + start_month = [12, 3, 6, 9] + end_month = [2, 5, 8, 11] + # Starting and ending days of the month (1-based) for each season + # Assume no leap years, we'll fix this later if we need + start_day = [1, 1, 1, 1] + end_day = [28, 31, 31, 30] + # Number of days in each season (again, ignore leap years for now) + ndays_season = [90, 92, 92, 91] + season_names = ['DJF', 'MAM', 'JJA', 'SON'] + + # Set up array to hold seasonal data + new_shape = [4] + shape + seasonal_data = ma.empty(new_shape) + seasonal_data[:] = 0.0 + + # Read the time values + id = Dataset(file_path, 'r') + time_id = id.variables['ocean_time'] + # Get the year, month, and day (all 1-based) for each output step + # These are 5-day averages marked with the middle day's date + time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) + + # Loop backwards through time indices to find the last one we care about + # (which contains 30 November in its averaging period) + end_t = -1 # Missing value flag + for t in range(size(time)-1, -1, -1): + if time[t].month == end_month[-1] and time[t].day in range(end_day[-1]-2, end_day[-1]+1): + end_t = t + break + if time[t].month == start_month[0] and time[t].day in range(start_day[0], start_day[0]+2): + end_t = t + break + # Make sure we actually found it + if end_t == -1: + print 'Error: ' + file_path + ' does not contain a complete Dec-Nov period' + return + + # Continue looping backwards to find the first time index we care about + # (which contains 1 December the previous year in its averaging period) + start_t = -1 # Missing value flag + for t in range(end_t-60, -1, -1): + if time[t].month == end_month[-1] and time[t].day in range(end_day[-1]-1, end_day[-1]+1): + start_t = t + break + if time[t].month == start_month[0] and time[t].day in range(start_day[0], start_day[0]+3): + start_t = t + break + # Make sure we actually found it + if start_t == -1: + print 'Error: ' + file_path + ' does not contain a complete Dec-Nov period' + return + + # Check if end_t occurs on a leap year + leap_year = False + if mod(time[end_t].year, 4) == 0: + # Years divisible by 4 are leap years + leap_year = True + if mod(time[end_t].year, 100) == 0: + # Unless they're also divisible by 100, in which case they aren't + # leap years + leap_year = False + if mod(time[end_t].year, 400) == 0: + # Unless they're also divisible by 400, in which case they are + # leap years after all + leap_year = True + if leap_year: + # Update last day in February + end_day[0] += 1 + ndays_season[0] += 1 + + # Process one season at a time + for season in range(4): + print 'Calculating seasonal averages for ' + season_names[season] + season_days = 0 # Number of days in season; this will be incremented + next_season = mod(season+1, 4) + + # Find starting timestep + start_t_season = -1 + for t in range(start_t, end_t+1): + if time[t].month == end_month[season-1] and time[t].day in range(end_day[season-1]-1, end_day[season-1]+1): + start_t_season = t + break + if time[t].month == start_month[season] and time[t].day in range(start_day[season], start_day[season]+3): + start_t_season = t + break + # Make sure we actually found it + if start_t_season == -1: + print 'Error: could not find starting timestep for season ' + season_title[season] + return + + # Find ending timestep + end_t_season = -1 + for t in range(start_t_season+1, end_t+1): + if time[t].month == end_month[season] and time[t].day in range(end_day[season]-2, end_day[season]+1): + end_t_season = t + break + if time[t].month == start_month[next_season] and time[t].day in range(start_day[next_season], start_day[next_season]+2): + end_t_season = t + break + # Make sure we actually found it + if end_t_season == -1: + print 'Error: could not find ending timestep for season ' + season_title[season] + return + + # Figure out how many of the 5 days averaged in start_t_season are + # actually within this season + if time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]+2: + # Starting day is in position 1 of 5; we care about all of them + start_days = 5 + elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]+1: + # Starting day is in position 2 of 5; we care about the last 4 + start_days = 4 + elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]: + # Starting day is in position 3 of 5; we care about the last 3 + start_days = 3 + elif time[start_t_season].month == end_month[season-1] and time[start_t_season].day == end_day[season-1]: + # Starting day is in position 4 of 5; we care about the last 2 + start_days = 2 + elif time[start_t_season].month == end_month[season-1] and time[start_t_season].day == end_day[season-1]-1: + # Starting day is in position 5 of 5; we care about the last 1 + start_days = 1 + else: + print 'Error for season ' + season_title[season] + ': starting index is month ' + str(time[start_t_season].month) + ', day ' + str(time[start_t_season].day) + return + + # Start accumulating data weighted by days + seasonal_data[season,:] += id.variables[var][start_t_season,:]*start_days + season_days += start_days + + # Between start_t_season and end_t_season, we want all the days + for t in range(start_t_season+1, end_t_season): + seasonal_data[season,:] += id.variables[var][t,:]*5 + season_days += 5 + + # Figure out how many of the 5 days averaged in end_t_season are + # actually within this season + if time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season]+1: + # Ending day is in position 1 of 5; we care about the first 1 + end_days = 1 + elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season]: + # Ending day is in position 2 of 5; we care about the first 2 + end_days = 2 + elif time[end_t_season].month == end_month[season] and time[end_t_season].day == end_day[season]: + # Ending day is in position 3 of 5; we care about the first 3 + end_days = 3 + elif time[end_t_season].month == end_month[season] and time[end_t_season].day == end_day[season]-1: + # Ending day is in position 4 of 5; we care about the first 4 + end_days = 4 + elif time[end_t_season].month == end_month[season] and time[end_t_season].day == end_day[season]-2: + # Ending day is in position 5 of 5; we care about all 5 + end_days = 5 + else: + print 'Error for season ' + season_title[season] + ': ending index is month ' + str(time[end_t_season].month) + ', day ' + str(time[end_t_season].day) + return + + seasonal_data[season,:] += id.variables[var][end_t_season,:]*end_days + season_days += end_days + + # Check that we got the correct number of days + if season_days != ndays_season[season]: + print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) + return + + # Finished accumulating data, now convert from sum to average + seasonal_data[season,:] /= season_days + + # Finished reading data + id.close() + + return seasonal_data + + + + diff --git a/sose_roms_seasonal.py b/sose_roms_seasonal.py index 488cc9b..bcddd3f 100644 --- a/sose_roms_seasonal.py +++ b/sose_roms_seasonal.py @@ -2,6 +2,9 @@ from numpy import * from matplotlib.pyplot import * from calc_z import * +from seasonal_avg_roms import * +from interp_lon_roms import * +from interp_lon_sose import * # Make a 4x2 plot comparing lat vs. depth slices of seasonally averaged # temperature or salinity at the given longitude, between ROMS (last year of @@ -27,15 +30,6 @@ def sose_roms_seasonal (file_path, var_name, lon0, depth_bdry, save=False, fig_n hc = 40 N = 31 - # Starting and ending months (1-based) for each season - start_month = [12, 3, 6, 9] - end_month = [2, 5, 8, 11] - # Starting and ending days of the month (1-based) for each season - # Assume no leap years, we'll fix this later if we need - start_day = [1, 1, 1, 1] - end_day = [28, 31, 31, 30] - # Number of days in each season (again, ignore leap years for now) - ndays_season = [90, 92, 92, 91] # Season names for titles season_names = ['DJF', 'MAM', 'JJA', 'SON'] @@ -68,167 +62,21 @@ def sose_roms_seasonal (file_path, var_name, lon0, depth_bdry, save=False, fig_n print 'Processing ROMS data' - # Read grid and time values + # Read grid id = Dataset(file_path, 'r') h = id.variables['h'][:-15,:] zice = id.variables['zice'][:-15,:] lon_roms_2d = id.variables['lon_rho'][:-15,:] lat_roms_2d = id.variables['lat_rho'][:-15,:] - time_id = id.variables['ocean_time'] - # Get the year, month, and day (all 1-based) for each output step - # These are 5-day averages marked with the middle day's date - time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - - # Loop backwards through time indices to find the last one we care about - # (which contains 30 November in its averaging period) - end_t = -1 # Missing value flag - for t in range(size(time)-1, -1, -1): - if time[t].month == end_month[-1] and time[t].day in range(end_day[-1]-2, end_day[-1]+1): - end_t = t - break - if time[t].month == start_month[0] and time[t].day in range(start_day[0], start_day[0]+2): - end_t = t - break - # Make sure we actually found it - if end_t == -1: - print 'Error: ' + file_path + ' does not contain a complete Dec-Nov period' - return - - # Continue looping backwards to find the first time index we care about - # (which contains 1 December the previous year in its averaging period) - start_t = -1 # Missing value falg - for t in range(end_t-60, -1, -1): - if time[t].month == end_month[-1] and time[t].day in range(end_day[-1]-1, end_day[-1]+1): - start_t = t - break - if time[t].month == start_month[0] and time[t].day in range(start_day[0], start_day[0]+3): - start_t = t - break - # Make sure we actually found it - if start_t == -1: - print 'Error: ' + file_path + ' does not contain a complete Dec-Nov period' - return - - # Check if end_t occurs on a leap year - leap_year = False - if mod(time[end_t].year, 4) == 0: - # Years divisible by 4 are leap years - leap_year = True - if mod(time[end_t].year, 100) == 0: - # Unless they're also divisible by 100, in which case they aren't - # leap years - leap_year = False - if mod(time[end_t].year, 400) == 0: - # Unless they're also divisible by 400, in which case they are - # leap years after all - leap_year = True - if leap_year: - # Update last day in February - end_day[0] += 1 - ndays_season[0] += 1 - - # Initialise seasonal averages - var_3d_roms = ma.empty([4, N, size(lon_roms_2d,0), size(lon_roms_2d,1)]) - var_3d_roms[:,:,:,:] = 0.0 - # Process one season at a time - for season in range(4): - print 'Calculating seasonal averages for ' + season_names[season] - season_days = 0 # Number of days in season; this will be incremented - next_season = mod(season+1, 4) - - # Find starting timestep - start_t_season = -1 - for t in range(start_t, end_t+1): - if time[t].month == end_month[season-1] and time[t].day in range(end_day[season-1]-1, end_day[season-1]+1): - start_t_season = t - break - if time[t].month == start_month[season] and time[t].day in range(start_day[season], start_day[season]+3): - start_t_season = t - break - # Make sure we actually found it - if start_t_season == -1: - print 'Error: could not find starting timestep for season ' + season_title[season] - return - - # Find ending timestep - end_t_season = -1 - for t in range(start_t_season+1, end_t+1): - if time[t].month == end_month[season] and time[t].day in range(end_day[season]-2, end_day[season]+1): - end_t_season = t - break - if time[t].month == start_month[next_season] and time[t].day in range(start_day[next_season], start_day[next_season]+2): - end_t_season = t - break - # Make sure we actually found it - if end_t_season == -1: - print 'Error: could not find ending timestep for season ' + season_title[season] - return - - # Figure out how many of the 5 days averaged in start_t_season are - # actually within this season - if time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]+2: - # Starting day is in position 1 of 5; we care about all of them - start_days = 5 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]+1: - # Starting day is in position 2 of 5; we care about the last 4 - start_days = 4 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]: - # Starting day is in position 3 of 5; we care about the last 3 - start_days = 3 - elif time[start_t_season].month == end_month[season-1] and time[start_t_season].day == end_day[season-1]: - # Starting day is in position 4 of 5; we care about the last 2 - start_days = 2 - elif time[start_t_season].month == end_month[season-1] and time[start_t_season].day == end_day[season-1]-1: - # Starting day is in position 5 of 5; we care about the last 1 - start_days = 1 - else: - print 'Error for season ' + season_title[season] + ': starting index is month ' + str(time[start_t_season].month) + ', day ' + str(time[start_t_season].day) - return - - # Start accumulating data weighted by days - var_3d_roms[season,:,:,:] += id.variables[var_name][start_t_season,:,:-15,:]*start_days - season_days += start_days - - # Between start_t_season and end_t_season, we want all the days - for t in range(start_t_season+1, end_t_season): - var_3d_roms[season,:,:,:] += id.variables[var_name][t,:,:-15,:]*5 - season_days += 5 - - # Figure out how many of the 5 days averaged in end_t_season are - # actually within this season - if time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season]+1: - # Ending day is in position 1 of 5; we care about the first 1 - end_days = 1 - elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season]: - # Ending day is in position 2 of 5; we care about the first 2 - end_days = 2 - elif time[end_t_season].month == end_month[season] and time[end_t_season].day == end_day[season]: - # Ending day is in position 3 of 5; we care about the first 3 - end_days = 3 - elif time[end_t_season].month == end_month[season] and time[end_t_season].day == end_day[season]-1: - # Ending day is in position 4 of 5; we care about the first 4 - end_days = 4 - elif time[end_t_season].month == end_month[season] and time[end_t_season].day == end_day[season]-2: - # Ending day is in position 5 of 5; we care about all 5 - end_days = 5 - else: - print 'Error for season ' + season_title[season] + ': ending index is month ' + str(time[end_t_season].month) + ', day ' + str(time[end_t_season].day) - return - - var_3d_roms[season,:,:,:] += id.variables[var_name][end_t_season,:,:-15,:]*end_days - season_days += end_days - - # Check that we got the correct number of days - if season_days != ndays_season[season]: - print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) - return - - # Finished accumulating data, now convert from sum to average - var_3d_roms[season,:,:,:] /= season_days - - # Finished reading data + num_lon = id.variables['lon_rho'].shape[1] + num_lat = id.variables['lat_rho'].shape[0] id.close() + # Calculate seasonal averages of ROMS data + var_3d_roms = seasonal_avg_roms(file_path, var_name, [N, num_lat, num_lon]) + # Chop off northern boundary + var_3d_roms = var_3d_roms[:,:,:-15,:] + # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script z_roms_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N) @@ -300,147 +148,6 @@ def sose_roms_seasonal (file_path, var_name, lon0, depth_bdry, save=False, fig_n fig.show() -# Linearly interpolate ROMS data, z, and latitude to the specified longitude. -# Input: -# data_3d = array of data, dimension depth x lat x lon -# z_3d = array of depth values (negative, in metres), dimension depth x lat x lon -# lat_2d = array of latitude values, dimension lat x lon -# lon_2d = array of longitude values, dimension lat x lon (between -180 and 180) -# lon0 = longitude to interpolate to (between -180 and 180) -# Output: -# data = array of data interpolated to lon0, dimension depth x lat -# z = array of depth values interpolated to lon0, dimension depth x lat -# lat = array of latitude values interpolated to lon0, dimension depth x lat -def interp_lon_roms (data_3d, z_3d, lat_2d, lon_2d, lon0): - - # Save dimensions - num_depth = size(data_3d, 0) - num_lat = size(data_3d, 1) - num_lon = size(data_3d, 2) - # Set up output arrays - data = ma.empty([num_depth, num_lat]) - z = ma.empty([num_depth, num_lat]) - lat = ma.empty([num_depth, num_lat]) - - # Loop over latitudes; can't find a cleaner way to do this - for j in range(num_lat): - # Extract the longitude values of this slice - lon_tmp = lon_2d[j,:] - # Get indices and coefficients for interpolation - ie, iw, coeffe, coeffw = interp_lon_roms_helper(lon_tmp, lon0) - data[:,j] = coeffe*data_3d[:,j,ie] + coeffw*data_3d[:,j,iw] - z[:,j] = coeffe*z_3d[:,j,ie] + coeffw*z_3d[:,j,iw] - lat[:,j] = coeffe*lat_2d[j,ie] + coeffw*lat_2d[j,iw] - - return data, z, lat - - -# Calculate indices and coefficients for linear interpolation of ROMS longitude. -# This takes care of all the mod 360 nonsense. -# Input: -# lon = 1D array of longitude values (straight out of ROMS i.e. between slightly < 0 and slightly > 360) -# lon0 = longitude to interpolate to (between 0 and 360) -# Output: -# ie, iw, coeffe, coeffw = integers (ie and iw) and coefficients (coeffe and -# coeffw) such that coeffe*lon[ie] + coeffw*lon[iw] = -# lon0, which will also hold for any variable on this -# longitude grid. ie is the index of the nearest point -# to the east of lon0; iw the nearest to the west. -def interp_lon_roms_helper (lon, lon0): - - if lon0 < amin(lon) or lon0 > amax(lon): - # Special case: lon0 on periodic boundary - # Be careful with mod 360 here - - # Find the periodic boundary - dlon = lon[1:] - lon[0:-1] - bdry = argmax(abs(dlon)) - if dlon[bdry] < -300: - # Jumps from almost 360 to just over 0 - iw = bdry - ie = bdry + 1 - else: - # Periodic boundary lines up with the array boundary - iw = size(lon) - 1 - ie = 0 - # Calculate difference between lon0 and lon[iw], mod 360 if necessary - dlon_num = lon0 - lon[iw] - if dlon_num < -300: - dlon_num += 360 - # Calculate difference between lon[ie] and lon[iw], mod 360 - dlon_den = lon[ie] - lon[iw] + 360 - - else: - # General case - - # Add or subtract 360 from longitude values which wrap around - # so that longitude increases monotonically from west to east - i = arange(1, size(lon)+1) - index1 = nonzero((i > 1200)*(lon < 100)) - lon[index1] = lon[index1] + 360 - index2 = nonzero((i < 200)*(lon > 300)) - lon[index2] = lon[index2] - 360 - - # Take mod 360 of lon0 if necessary - if all(lon < lon0): - lon0 -= 360 - if all(lon > lon0): - lon0 += 360 - - # Find the first index eastward of lon0 - ie = nonzero(lon > lon0)[0][0] - # The index before it will be the last index westward of lon0 - iw = ie - 1 - - dlon_num = lon0 - lon[iw] - dlon_den = lon[ie] - lon[iw] - - if dlon_num > 5 or dlon_den > 5: - print 'interp_lon_helper: Problem at periodic boundary' - return - coeff1 = dlon_num/dlon_den - coeff2 = 1 - coeff1 - - return ie, iw, coeff1, coeff2 - - -# Linearly interpolate SOSE data to the specified longitude. -# Input: -# data_3d = arary of data, dimension depth x lat x lon -# lon = 1D array of longitude values (between 0 and 360) -# lon0 = longitude to interpolate to (between 0 and 360) -# Output: -# data = array of data interpolated to lon0, dimension depth x lat -def interp_lon_sose (data_3d, lon, lon0): - - if lon0 < lon[0] or lon0 > lon[-1]: - # Special case: lon0 on periodic boundary - # Be careful with mod 360 here - iw = size(lon)-1 - ie = 0 - # Calculate difference between lon0 and lon[iw], mod 360 if necessary - dlon_num = lon0 - lon[iw] - if dlon_num < -300: - dlon_num += 360 - # Calculate difference between lon[ie] and lon[iw], mod 360 - dlon_den = lon[ie] - lon[iw] + 360 - else: - # General case - # Find the first index eastwards of lon0 - ie = nonzero(lon > lon0)[0][0] - # The index before it will be the last index westward of lon0 - iw = ie - 1 - dlon_num = lon0 - lon[iw] - dlon_den = lon[ie] - lon[iw] - # Coefficients for interpolation - coeff1 = dlon_num/dlon_den - coeff2 = 1 - coeff1 - - # Interpolate - data = coeff1*data_3d[:,:,ie] + coeff2*data_3d[:,:,iw] - return data - - # Command-line interface if __name__ == "__main__": diff --git a/sose_seasonal.py b/sose_seasonal.py index a71bbad..50f63a9 100644 --- a/sose_seasonal.py +++ b/sose_seasonal.py @@ -2,6 +2,7 @@ from numpy import * from matplotlib.pyplot import * from calc_z import * +from interp_lon_sose import * # Make a 4x2 plot showing lat vs. depth slices of seasonally averaged # temperature (top) and salinity (bottom) at the given longitude, for the @@ -106,43 +107,6 @@ def sose_seasonal (lon0, depth_bdry, save=False, fig_name=None): fig.show() -# Linearly interpolate SOSE data to the specified longitude. -# Input: -# data_3d = arary of data, dimension depth x lat x lon -# lon = 1D array of longitude values (between 0 and 360) -# lon0 = longitude to interpolate to (between 0 and 360) -# Output: -# data = array of data interpolated to lon0, dimension depth x lat -def interp_lon_sose (data_3d, lon, lon0): - - if lon0 < lon[0] or lon0 > lon[-1]: - # Special case: lon0 on periodic boundary - # Be careful with mod 360 here - iw = size(lon)-1 - ie = 0 - # Calculate difference between lon0 and lon[iw], mod 360 if necessary - dlon_num = lon0 - lon[iw] - if dlon_num < -300: - dlon_num += 360 - # Calculate difference between lon[ie] and lon[iw], mod 360 - dlon_den = lon[ie] - lon[iw] + 360 - else: - # General case - # Find the first index eastwards of lon0 - ie = nonzero(lon > lon0)[0][0] - # The index before it will be the last index westward of lon0 - iw = ie - 1 - dlon_num = lon0 - lon[iw] - dlon_den = lon[ie] - lon[iw] - # Coefficients for interpolation - coeff1 = dlon_num/dlon_den - coeff2 = 1 - coeff1 - - # Interpolate - data = coeff1*data_3d[:,:,ie] + coeff2*data_3d[:,:,iw] - return data - - # Command-line interface if __name__ == "__main__": diff --git a/ssflux_monthly.py b/ssflux_monthly.py index 131ba86..23fec3d 100644 --- a/ssflux_monthly.py +++ b/ssflux_monthly.py @@ -1,112 +1,24 @@ from netCDF4 import Dataset, num2date from numpy import * from matplotlib.pyplot import * +from monthly_avg_roms import * def ssflux_monthly (roms_file, month, bound=None, save=False, fig_name=None): # Month names for plot titles month_name = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] - # Number of days per month - start_day = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] - end_day = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] deg2rad = pi/180 + # Read ROMS grid id = Dataset(roms_file, 'r') lon = id.variables['lon_rho'][:-15,:-1] lat = id.variables['lat_rho'][:-15,:-1] - time_id = id.variables['ocean_time'] - # Get the year, month, and day (all 1-based) for each output step - # These are 5-day averages marked with the middle day's date - time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - - next_month = mod(month+1,12) - prev_month = mod(month-1,12) - - end_t = -1 - for t in range(size(time)-1, -1, -1): - if time[t].month-1 == month and time[t].day in range(end_day[month]-2, end_day[month]+1): - end_t = t - break - if time[t].month-1 == next_month and time[t].day in range(start_day[next_month], start_day[next_month]+2): - end_t = t - break - if end_t == -1: - print 'Error: ' + roms_file + ' does not contain a complete ' + month_name[month] - return - - start_t = -1 - for t in range(end_t, -1, -1): - if time[t].month-1 == prev_month and time[t].day in range(end_day[prev_month]-1, end_day[prev_month]+1): - start_t = t - break - if time[t].month-1 == month and time[t].day in range(start_day[month], start_day[month]+3): - start_t = t - break - if start_t == -1: - print 'Error: ' + roms_file + ' does not contain a complete ' + month_name[month] - return - - leap_year = False - roms_year = time[end_t].year - if month == 11: - roms_year = time[start_t].year - if mod(roms_year, 4) == 0: - leap_year = True - if mod(roms_year, 100) == 0: - leap_year = False - if mod(roms_year, 400) == 0: - leap_year = True - if leap_year: - end_day[1] = 29 - - ssflux = ma.empty(shape(lon)) - ssflux[:,:] = 0.0 - num_days = 0 - - if time[start_t].month-1 == month and time[start_t].day == start_day[month]+2: - start_days = 5 - elif time[start_t].month-1 == month and time[start_t].day == start_day[month]+1: - start_days = 4 - elif time[start_t].month-1 == month and time[start_t].day == start_day[month]: - start_days = 3 - elif time[start_t].month-1 == prev_month and time[start_t].day == end_day[prev_month]: - start_days = 2 - elif time[start_t].month-1 == prev_month and time[start_t].day == end_day[prev_month]-1: - start_days = 1 - else: - print 'Error: starting index is month ' + str(time[start_t].month) + ', day ' + str(time[start_t].day) - return - - ssflux += id.variables['ssflux'][start_t,:-15,:-1]*start_days - num_days += start_days - - for t in range(start_t+1, end_t): - ssflux += id.variables['ssflux'][t,:-15,:-1]*5 - num_days += 5 - - if time[end_t].month-1 == next_month and time[end_t].day == start_day[next_month]+1: - end_days = 1 - elif time[end_t].month-1 == next_month and time[end_t].day == start_day[next_month]: - end_days = 2 - elif time[end_t].month-1 == month and time[end_t].day == end_day[month]: - end_days = 3 - elif time[end_t].month-1 == month and time[end_t].day == end_day[month]-1: - end_days = 4 - elif time[end_t].month-1 == month and time[end_t].day == end_day[month]-2: - end_days = 5 - else: - print 'Error: ending index is month ' + str(time[end_t].month) + ', day ' + str(time[end_t].day) - return - - ssflux += id.variables['ssflux'][end_t,:-15,:-1]*end_days - num_days += end_days - - if num_days != end_day[month]: - print 'Error: found ' + str(num_days) + ' days instead of ' + str(end_day[month]) - return - + num_lon = id.variables['lon_rho'].shape[1] + num_lat = id.variables['lat_rho'].shape[0] id.close() - ssflux /= num_days + + ssflux = monthly_avg_roms(roms_file, 'ssflux', [num_lat, num_lon], month) + ssflux = ssflux[:-15,:-1] ssflux *= 1e6 x = -(lat+90)*cos(lon*deg2rad+pi/2) diff --git a/ssflux_tamura_monthly.py b/ssflux_tamura_monthly.py index 129849f..2468049 100644 --- a/ssflux_tamura_monthly.py +++ b/ssflux_tamura_monthly.py @@ -1,6 +1,7 @@ from numpy import * from netCDF4 import Dataset, num2date from matplotlib.pyplot import * +from monthly_avg_cice import * # Make a figure comparing monthly-averaged sea ice to ocean salt fluxes, # from Tamura's dataset to CICE output. @@ -8,18 +9,15 @@ # cice_file = path to CICE output file with 5-day averages; if it covers more # than one instance of the given month, plot the last one # month = month number (0-indexed) from 0 to 11 +# cice_year = year that this month in cice_file corresponds to # save = optional boolean to save the figure to a file, rather than displaying # it on the screen # fig_name = if save=True, path to the desired filename for figure -def ssflux_tamura_monthly (cice_file, month, save=False, fig_name=None): +def ssflux_tamura_monthly (cice_file, month, cice_year, save=False, fig_name=None): # Beginning and end of Tamura file paths tamura_head = '/short/m68/kaa561/tamura_fluxes/Tamura_ssflux_' tamura_tail = '_monthly.nc' - # Starting and ending days in each month - # Assume no leap years, we'll fix this later if we need - start_day = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] - end_day = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] # Name of each month, for the title month_name = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] # Degrees to radians conversion @@ -31,10 +29,13 @@ def ssflux_tamura_monthly (cice_file, month, save=False, fig_name=None): # Conversion factor: m/s to cm/day mps_to_cmpday = 8.64e6 - # Read the CICE grid and time values + # Read the CICE grid id = Dataset(cice_file, 'r') cice_lon_tmp = id.variables['TLON'][:-15,:] cice_lat_tmp = id.variables['TLAT'][:-15,:] + num_lon = id.variables['TLON'].shape[1] + num_lat = id.variables['TLAT'].shape[0] + id.close() # Wrap the periodic boundary by 1 cell cice_lon = ma.empty([size(cice_lon_tmp,0), size(cice_lon_tmp,1)+1]) cice_lat = ma.empty([size(cice_lat_tmp,0), size(cice_lat_tmp,1)+1]) @@ -42,142 +43,15 @@ def ssflux_tamura_monthly (cice_file, month, save=False, fig_name=None): cice_lon[:,-1] = cice_lon_tmp[:,0] cice_lat[:,:-1] = cice_lat_tmp cice_lat[:,-1] = cice_lat_tmp[:,0] - time_id = id.variables['time'] - # Get the year, month, and day (all 1-indexed) for each output step - # These are 5-day averages marked with the next day's date. - cice_time = num2date(time_id[:], units=time_id.units, calendar='standard') - - # Get index of the next month - next_month = mod(month+1, 12) - # Loop backwards through time indices to find the last one we care about - # (which contains the last day of the month in its averaging period) - end_t = -1 # Missing value flag - for t in range(size(cice_time)-1, -1, -1): - # Note that cice_time[t].month is converted to 0-indexed - if cice_time[t].month-1 == next_month and cice_time[t].day in range(start_day[next_month], start_day[next_month]+5): - end_t = t - break - # Make sure we actually found it - if end_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete ' + month_name[month] - return - - # Continue looping backwards to find the first time index we care about - start_t = -1 # Missing value flag - for t in range(end_t, -1, -1): - if cice_time[t].month-1 == month and cice_time[t].day in range(start_day[month]+1, start_day[month]+6): - start_t = t - break - # Make sure we actually found it - if start_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete ' + month_name[month] - return - - # Check if this is a leap year - leap_year = False - cice_year = cice_time[end_t].year - if month == 11: - # Timestep end_t is likely in the next year, so find the year of start_t - cice_year = cice_time[start_t].year - if mod(cice_year, 4) == 0: - # Years divisible by 4 are leap years - leap_year = True - if mod(cice_year, 100) == 0: - # Unless they're also divisible by 100, in which case they aren't - # leap years - leap_year = False - if mod(cice_year, 400) == 0: - # Unless they're also divisible by 400, in which case they are - # leap years after all - leap_year = True - if leap_year: - # Update last day in February - end_day[1] = 29 - - # Calculate monthly average of CICE output - cice_data_tmp = ma.empty(shape(cice_lon_tmp)) - cice_data_tmp[:,:] = 0.0 - num_days = 0 - - # Figure out how many of the 5 days averaged in start_t are actually within - # this month - if cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+5: - # Starting day is in position 1 of 5; we care about all of them - start_days = 5 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+4: - # Starting day is in position 2 of 5; we care about the last 4 - start_days = 4 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+3: - # Starting day is in position 3 of 5; we care about the last 3 - start_days = 3 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+2: - # Starting day is in position 4 of 5; we care about the last 2 - start_days = 2 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+1: - # Starting day is in position 5 of 5; we care about the last 1 - start_days = 1 - else: - print 'Error: starting index is month ' + str(cice_time[start_t].month) + ', day ' + str(cice_time[start_t].day) - return - - # Read the fields we need at start_t - fresh_ai = id.variables['fresh_ai'][start_t,:-15,:] - sss = id.variables['sss'][start_t,:-15,:] - rain_ai = id.variables['rain_ai'][start_t,:-15,:] - fsalt_ai = id.variables['fsalt_ai'][start_t,:-15,:] - # Start accumulating data weighted by days - # Convert to units of psu m/s (equivalent to kg/m^2/s of salt) - # Subtract rain from freshwater flux, since Tamura doesn't count precip - cice_data_tmp += -1/rho_fw*((fresh_ai-rain_ai)*sss*rho_sw/mps_to_cmpday - fsalt_ai*1e3)*start_days - num_days += start_days - - # Between start_t and end_t, we want all the days - for t in range(start_t+1, end_t): - fresh_ai = id.variables['fresh_ai'][t,:-15,:] - sss = id.variables['sss'][t,:-15,:] - rain_ai = id.variables['rain_ai'][t,:-15,:] - fsalt_ai = id.variables['fsalt_ai'][t,:-15,:] - cice_data_tmp += -1/rho_fw*((fresh_ai-rain_ai)*sss*rho_sw/mps_to_cmpday - fsalt_ai*1e3)*5 - num_days += 5 - - # Figure out how many of the 5 days averaged in end_t are actually within - # this month - if cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+4: - # Ending day is in position 1 of 5; we care about the first 1 - end_days = 1 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+3: - # Ending day is in position 2 of 5; we care about the first 2 - end_days = 2 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+2: - # Ending day is in position 3 of 5; we care about the first 3 - end_days = 3 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+1: - # Ending day is in position 4 of 5; we care about the first 4 - end_days = 4 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]: - # Ending day is in position 5 of 5; we care about all 5 - end_days = 5 - else: - print 'Error: ending index is month ' + str(cice_time[end_t].month) + ', day ' + str(cice_time[end_t].day) - return - - fresh_ai = id.variables['fresh_ai'][end_t,:-15,:] - sss = id.variables['sss'][end_t,:-15,:] - rain_ai = id.variables['rain_ai'][end_t,:-15,:] - fsalt_ai = id.variables['fsalt_ai'][end_t,:-15,:] - cice_data_tmp += -1/rho_fw*((fresh_ai-rain_ai)*sss*rho_sw/mps_to_cmpday - fsalt_ai*1e3)*end_days - num_days += end_days - - # Check that we got the correct number of days - if num_days != end_day[month]: - print 'Error: found ' + str(num_days) + ' days instead of ' + str(end_day[month]) - return - - # Finished accumulating data - id.close() - # Convert from sum to average - cice_data_tmp /= num_days + # Get monthly averages of each variable we need + fresh_ai = monthly_avg_cice(cice_file, 'fresh_ai', [num_lat, num_lon], month) + sss = monthly_avg_cice(cice_file, 'sss', [num_lat, num_lon], month) + rain_ai = monthly_avg_cice(cice_file, 'rain_ai', [num_lat, num_lon], month) + fsalt_ai = monthly_avg_cice(cice_file, 'fsalt_ai', [num_lat, num_lon], month) + cice_data_tmp = -1/rho_fw*((fresh_ai-rain_ai)*sss*rho_sw/mps_to_cmpday - fsalt_ai*1e3) + # Chop off northern boundary + cice_data_tmp = cice_data_tmp[:-15,:] # Multiply by 1e6 so colour bar is easier to read cice_data_tmp *= 1e6 @@ -247,6 +121,7 @@ def ssflux_tamura_monthly (cice_file, month, save=False, fig_name=None): cice_file = raw_input("Path to CICE file: ") # Convert month from 1-indexed to 0-indexed month = int(raw_input("Month number (1-12): "))-1 + cice_year = int(raw_input("Year this month corresponds to in the CICE file: ")) action = raw_input("Save figure (s) or display on screen (d)? ") if action == 's': save = True @@ -255,7 +130,7 @@ def ssflux_tamura_monthly (cice_file, month, save=False, fig_name=None): save = False fig_name = None # Make the plot - ssflux_tamura_monthly(cice_file, month, save, fig_name) + ssflux_tamura_monthly(cice_file, month, cice_year, save, fig_name) while True: # Repeat until the user wants to exit @@ -263,7 +138,7 @@ def ssflux_tamura_monthly (cice_file, month, save=False, fig_name=None): if repeat == 'y': while True: # Ask for changes to the parameters until the user is done - changes = raw_input("Enter a parameter to change: (1) file path, (2) month number, (3) save/display; or enter to continue: ") + changes = raw_input("Enter a parameter to change: (1) file path, (2) month number, (3) year, (4) save/display; or enter to continue: ") if len(changes) == 0: # No more changes to parameters break @@ -275,13 +150,15 @@ def ssflux_tamura_monthly (cice_file, month, save=False, fig_name=None): # New month month = int(raw_input("Month number (1-12): "))-1 elif int(changes) == 3: + cice_year = int(raw_input("Year this month corresponds to in the CICE file: ")) + elif int(changes) == 4: # Switch from save to display, or vice versa save = not save if save: # Get a new figure name fig_name = raw_input("File name for figure: ") # Make the plot - ssflux_tamura_monthly(cice_file, month, save, fig_name) + ssflux_tamura_monthly(cice_file, month, cice_year, save, fig_name) else: break diff --git a/ssflux_tamura_seasonal.py b/ssflux_tamura_seasonal.py index 3bbfa08..ccf9833 100644 --- a/ssflux_tamura_seasonal.py +++ b/ssflux_tamura_seasonal.py @@ -1,6 +1,7 @@ from numpy import * -from netCDF4 import Dataset, num2date +from netCDF4 import Dataset from matplotlib.pyplot import * +from seasonal_avg_cice import * # Make a figure comparing seasonally-averaged sea ice to ocean salt fluxes, # from Tamura's dataset to CICE output. @@ -8,23 +9,15 @@ # cice_file = path to CICE output file with 5-day averages, containing at least # one complete Dec-Nov period (if there are multiple such periods, # this script uses the last one) +# cice_year = year that the last Jan-Nov period in CICE corresponds to # save = optional boolean to save the figure to a file, rather than displaying # it on the screen # fig_name = if save=True, path to the desired filename for figure -def ssflux_tamura_seasonal (cice_file, save=False, fig_name=None): +def ssflux_tamura_seasonal (cice_file, cice_year, save=False, fig_name=None): # Beginning and end of Tamura file paths tamura_head = '/short/m68/kaa561/tamura_fluxes/Tamura_ssflux_' tamura_tail = '_monthly.nc' - # Starting and ending months (1-based) for each season - start_month = [12, 3, 6, 9] - end_month = [2, 5, 8, 11] - # Starting and ending days of the month (1-based) for each season - # Assume no leap years, we'll fix this later if we need - start_day = [1, 1, 1, 1] - end_day = [28, 31, 31, 30] - # Number of days in each season (again, ignore leap years for now) - ndays_season = [90, 92, 92, 91] # Number of days in each month (this is just for Tamura) ndays_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] # Season names for titles @@ -38,10 +31,13 @@ def ssflux_tamura_seasonal (cice_file, save=False, fig_name=None): # Conversion factor: m/s to cm/day mps_to_cmpday = 8.64e6 - # Read the CICE grid and time values + # Read the CICE grid id = Dataset(cice_file, 'r') cice_lon_tmp = id.variables['TLON'][:-15,:] cice_lat_tmp = id.variables['TLAT'][:-15,:] + num_lon = id.variables['TLON'].shape[1] + num_lat = id.variables['TLAT'].shape[0] + id.close() # Wrap the periodic boundary by 1 cell cice_lon = ma.empty([size(cice_lon_tmp,0), size(cice_lon_tmp,1)+1]) cice_lat = ma.empty([size(cice_lat_tmp,0), size(cice_lat_tmp,1)+1]) @@ -49,38 +45,34 @@ def ssflux_tamura_seasonal (cice_file, save=False, fig_name=None): cice_lon[:,-1] = cice_lon_tmp[:,0] cice_lat[:,:-1] = cice_lat_tmp cice_lat[:,-1] = cice_lat_tmp[:,0] - time_id = id.variables['time'] - # Get the year, month, and day (all 1-indexed) for each output step - # These are 5-day averages marked with the next day's date. - cice_time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - # Loop backwards through time indices to find the last one we care about - # (which contains 30 November in its averaging period) - end_t = -1 # Missing value flag - for t in range(size(cice_time)-1, -1, -1): - if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+5): - end_t = t - break - # Make sure we actually found it - if end_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' - return + # Read seasonally averaged fields + print 'Reading fresh_ai' + fresh_ai = seasonal_avg_cice(cice_file, 'fresh_ai', [num_lat, num_lon]) + print 'Reading sss' + sss = seasonal_avg_cice(cice_file, 'sss', [num_lat, num_lon]) + print 'Reading rain_ai' + rain_ai = seasonal_avg_cice(cice_file, 'rain_ai', [num_lat, num_lon]) + print 'Reading fsalt_ai' + fsalt_ai = seasonal_avg_cice(cice_file, 'fsalt_ai', [num_lat, num_lon]) + cice_data_tmp = -1/rho_fw*((fresh_ai-rain_ai)*sss*rho_sw/mps_to_cmpday - fsalt_ai*1e3) + # Multiply by 1e6 so colour bar is easier to read + cice_data_tmp *= 1e6 + # Chop off northern boundary + cice_data_tmp = cice_data_tmp[:,:-15,:] + # Wrap periodic boundary + cice_data = ma.empty([size(cice_data_tmp,0), size(cice_data_tmp,1), size(cice_data_tmp,2)+1]) + cice_data[:,:,:-1] = cice_data_tmp + cice_data[:,:,-1] = cice_data_tmp[:,:,0] - # Continue looping backwards to find the first time index we care about - # (which contains 1 December the previous year in its averaging period) - start_t = -1 # Missing value flag - for t in range(end_t-60, -1, -1): - if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0]+1, start_day[0]+6): - start_t = t - break - # Make sure we actually found it - if start_t == -1: - print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' - return + # Read Tamura grid + id = Dataset(tamura_head + str(cice_year) + tamura_tail, 'r') + tamura_lon = id.variables['longitude'][:,:] + tamura_lat = id.variables['latitude'][:,:] + id.close() # Check for leap years leap_year = False - cice_year = cice_time[end_t].year if mod(cice_year, 4) == 0: # Years divisible by 4 are leap years leap_year = True @@ -94,133 +86,8 @@ def ssflux_tamura_seasonal (cice_file, save=False, fig_name=None): leap_year = True if leap_year: # Update last day in February - end_day[0] += 1 - ndays_season[0] += 1 ndays_month[1] += 1 - # Initialise seasonal averages of CICE output - cice_data_tmp = ma.empty([4, size(cice_lon_tmp,0), size(cice_lon_tmp,1)]) - cice_data_tmp[:,:,:] = 0.0 - - # Process one season at a time - for season in range(4): - season_days = 0 # Number of days in season; this will be incremented - next_season = mod(season+1, 4) - - # Find starting timestep - start_t_season = -1 - for t in range(start_t, end_t+1): - if cice_time[t].month == start_month[season] and cice_time[t].day in range(start_day[season]+1, start_day[season]+6): - start_t_season = t - break - # Make sure we actually found it - if start_t_season == -1: - print 'Error: could not find starting timestep for season ' + season_names[season] - return - - # Find ending timestep - end_t_season = -1 - for t in range(start_t_season+1, end_t+1): - if cice_time[t].month == start_month[next_season] and cice_time[t].day in range(start_day[next_season], start_day[next_season]+5): - end_t_season = t - break - # Make sure we actually found it - if end_t_season == -1: - print 'Error: could not find ending timestep for season ' + season_names[season] - return - - # Figure out how many of the 5 days averaged in start_t_season are - # actually within this season - if cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 5: - # Starting day is in position 1 of 5; we care about all of them - start_days = 5 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 4: - # Starting day is in position 2 of 5; we care about the last 4 - start_days = 4 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]+ 3: - # Starting day is in position 3 of 5; we care about the last 3 - start_days = 3 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 2: - # Starting day is in position 4 of 5; we care about the last 2 - start_days = 2 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 1: - # Starting day is in position 5 of 5; we care about the last 1 - start_days = 1 - else: - print 'Error for season ' + season_names[season] + ': starting index is month ' + str(cice_time[start_t_season].month) + ', day ' + str(cice_time[start_t_season].day) - return - - # Read the fields we need at start_t_season - fresh_ai = id.variables['fresh_ai'][start_t_season,:-15,:] - sss = id.variables['sss'][start_t_season,:-15,:] - rain_ai = id.variables['rain_ai'][start_t_season,:-15,:] - fsalt_ai = id.variables['fsalt_ai'][start_t_season,:-15,:] - # Start accumulating data weighted by days - # Convert to units of psu m/s (equivalent to kg/m^2/s of salt) - # Subtract rain from freshwater flux, since Tamura doesn't count precip - cice_data_tmp[season,:,:] += -1/rho_fw*((fresh_ai-rain_ai)*sss*rho_sw/mps_to_cmpday - fsalt_ai*1e3)*start_days - season_days += start_days - - # Between start_t_season and end_t_season, we want all the days - for t in range(start_t_season+1, end_t_season): - fresh_ai = id.variables['fresh_ai'][t,:-15,:] - sss = id.variables['sss'][t,:-15,:] - rain_ai = id.variables['rain_ai'][t,:-15,:] - fsalt_ai = id.variables['fsalt_ai'][t,:-15,:] - cice_data_tmp[season,:,:] += -1/rho_fw*((fresh_ai-rain_ai)*sss*rho_sw/mps_to_cmpday - fsalt_ai*1e3)*5 - season_days += 5 - - # Figure out how many of the 5 days averaged in end_t_season are - # actually within this season - if cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 4: - # Ending day is in position 1 of 5; we care about the first 1 - end_days = 1 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 3: - # Ending day is in position 2 of 5; we care about the first 2 - end_days = 2 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 2: - # Ending day is in position 3 of 5; we care about the first 3 - end_days = 3 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 1: - # Ending day is in position 4 of 5; we care about the first 4 - end_days = 4 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season]: - # Ending day is in position 5 of 5; we care about all 5 - end_days = 5 - else: - print 'Error for season ' + season_names[season] + ': ending index is month ' + str(cice_time[end_t_season].month) + ', day ' + str(cice_time[end_t_season].day) - return - - fresh_ai = id.variables['fresh_ai'][end_t_season,:-15,:] - sss = id.variables['sss'][end_t_season,:-15,:] - rain_ai = id.variables['rain_ai'][end_t_season,:-15,:] - fsalt_ai = id.variables['fsalt_ai'][end_t_season,:-15,:] - cice_data_tmp[season,:,:] += -1/rho_fw*((fresh_ai-rain_ai)*sss*rho_sw/mps_to_cmpday - fsalt_ai*1e3)*end_days - season_days += end_days - - # Check that we got the correct number of days - if season_days != ndays_season[season]: - print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) - return - - # Finished accumulating data, now convert from sum to average - cice_data_tmp[season,:,:] /= season_days - - id.close() - # Multiply by 1e6 so colour bar is easier to read - cice_data_tmp *= 1e6 - - # Wrap periodic boundary - cice_data = ma.empty([size(cice_data_tmp,0), size(cice_data_tmp,1), size(cice_data_tmp,2)+1]) - cice_data[:,:,:-1] = cice_data_tmp - cice_data[:,:,-1] = cice_data_tmp[:,:,0] - - # Read Tamura grid - id = Dataset(tamura_head + str(cice_year) + tamura_tail, 'r') - tamura_lon = id.variables['longitude'][:,:] - tamura_lat = id.variables['latitude'][:,:] - id.close() - # Set up array for seasonal averages of Tamura data tamura_data = ma.empty([4, size(tamura_lon,0), size(tamura_lon,1)]) tamura_data[:,:,:] = 0.0 @@ -311,6 +178,7 @@ def ssflux_tamura_seasonal (cice_file, save=False, fig_name=None): if __name__ == "__main__": cice_file = raw_input("Path to CICE file, containing at least one complete Dec-Nov period: ") + cice_year = int(raw_input("Year that the last complete Jan-Nov period in CICE file corresponds to: ")) action = raw_input("Save figure (s) or display on screen (d)? ") if action == 's': save = True @@ -318,4 +186,4 @@ def ssflux_tamura_seasonal (cice_file, save=False, fig_name=None): elif action == 'd': save = False fig_name = None - ssflux_tamura_seasonal(cice_file, save, fig_name) + ssflux_tamura_seasonal(cice_file, cice_year, save, fig_name) diff --git a/temp_salt_seasonal.py b/temp_salt_seasonal.py index 2805b9f..6c526e4 100644 --- a/temp_salt_seasonal.py +++ b/temp_salt_seasonal.py @@ -1,7 +1,9 @@ -from netCDF4 import Dataset, num2date +from netCDF4 import Dataset from numpy import * from matplotlib.pyplot import * from calc_z import * +from seasonal_avg_roms import * +from interp_lon_roms import * # Make a 4x2 plot showing lat vs. depth slices of seasonally averaged # temperature (top) and salinity (bottom) at the given longitude, over the @@ -17,24 +19,12 @@ # fig_name = optional string containing filename for figure, if save=True def temp_salt_seasonal (file_path, lon0, depth_bdry, save=False, fig_name=None): - # Path to SOSE seasonal climatology file - sose_file = '../SOSE_seasonal_climatology.nc' - # Grid parameters theta_s = 4.0 theta_b = 0.9 hc = 40 N = 31 - # Starting and ending months (1-based) for each season - start_month = [12, 3, 6, 9] - end_month = [2, 5, 8, 11] - # Starting and ending days of the month (1-based) for each season - # Assume no leap years, we'll fix this later if we need - start_day = [1, 1, 1, 1] - end_day = [28, 31, 31, 30] - # Number of days in each season (again, ignore leap years for now) - ndays_season = [90, 92, 92, 91] # Season names for titles season_names = ['DJF', 'MAM', 'JJA', 'SON'] @@ -55,175 +45,25 @@ def temp_salt_seasonal (file_path, lon0, depth_bdry, save=False, fig_name=None): if lon0 < 0: lon0 += 360 - print 'Processing ROMS data' - - # Read grid and time values + # Read grid id = Dataset(file_path, 'r') h = id.variables['h'][:-15,:] zice = id.variables['zice'][:-15,:] lon_roms_2d = id.variables['lon_rho'][:-15,:] lat_roms_2d = id.variables['lat_rho'][:-15,:] - time_id = id.variables['ocean_time'] - # Get the year, month, and day (all 1-based) for each output step - # These are 5-day averages marked with the middle day's date - time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - - # Loop backwards through time indices to find the last one we care about - # (which contains 30 November in its averaging period) - end_t = -1 # Missing value flag - for t in range(size(time)-1, -1, -1): - if time[t].month == end_month[-1] and time[t].day in range(end_day[-1]-2, end_day[-1]+1): - end_t = t - break - if time[t].month == start_month[0] and time[t].day in range(start_day[0], start_day[0]+2): - end_t = t - break - # Make sure we actually found it - if end_t == -1: - print 'Error: ' + file_path + ' does not contain a complete Dec-Nov period' - return - - # Continue looping backwards to find the first time index we care about - # (which contains 1 December the previous year in its averaging period) - start_t = -1 # Missing value flag - for t in range(end_t-60, -1, -1): - if time[t].month == end_month[-1] and time[t].day in range(end_day[-1]-1, end_day[-1]+1): - start_t = t - break - if time[t].month == start_month[0] and time[t].day in range(start_day[0], start_day[0]+3): - start_t = t - break - # Make sure we actually found it - if start_t == -1: - print 'Error: ' + file_path + ' does not contain a complete Dec-Nov period' - return - - # Check if end_t occurs on a leap year - leap_year = False - if mod(time[end_t].year, 4) == 0: - # Years divisible by 4 are leap years - leap_year = True - if mod(time[end_t].year, 100) == 0: - # Unless they're also divisible by 100, in which case they aren't - # leap years - leap_year = False - if mod(time[end_t].year, 400) == 0: - # Unless they're also divisible by 400, in which case they are - # leap years after all - leap_year = True - if leap_year: - # Update last day in February - end_day[0] += 1 - ndays_season[0] += 1 - - # Initialise seasonal averages - temp_3d_roms = ma.empty([4, N, size(lon_roms_2d,0), size(lon_roms_2d,1)]) - temp_3d_roms[:,:,:,:] = 0.0 - salt_3d_roms = ma.empty([4, N, size(lon_roms_2d,0), size(lon_roms_2d,1)]) - salt_3d_roms[:,:,:,:] = 0.0 - # Process one season at a time - for season in range(4): - print 'Calculating seasonal averages for ' + season_names[season] - season_days = 0 # Number of days in season; this will be incremented - next_season = mod(season+1, 4) - - # Find starting timestep - start_t_season = -1 - for t in range(start_t, end_t+1): - if time[t].month == end_month[season-1] and time[t].day in range(end_day[season-1]-1, end_day[season-1]+1): - start_t_season = t - break - if time[t].month == start_month[season] and time[t].day in range(start_day[season], start_day[season]+3): - start_t_season = t - break - # Make sure we actually found it - if start_t_season == -1: - print 'Error: could not find starting timestep for season ' + season_title[season] - return - - # Find ending timestep - end_t_season = -1 - for t in range(start_t_season+1, end_t+1): - if time[t].month == end_month[season] and time[t].day in range(end_day[season]-2, end_day[season]+1): - end_t_season = t - break - if time[t].month == start_month[next_season] and time[t].day in range(start_day[next_season], start_day[next_season]+2): - end_t_season = t - break - # Make sure we actually found it - if end_t_season == -1: - print 'Error: could not find ending timestep for season ' + season_title[season] - return - - # Figure out how many of the 5 days averaged in start_t_season are - # actually within this season - if time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]+2: - # Starting day is in position 1 of 5; we care about all of them - start_days = 5 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]+1: - # Starting day is in position 2 of 5; we care about the last 4 - start_days = 4 - elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]: - # Starting day is in position 3 of 5; we care about the last 3 - start_days = 3 - elif time[start_t_season].month == end_month[season-1] and time[start_t_season].day == end_day[season-1]: - # Starting day is in position 4 of 5; we care about the last 2 - start_days = 2 - elif time[start_t_season].month == end_month[season-1] and time[start_t_season].day == end_day[season-1]-1: - # Starting day is in position 5 of 5; we care about the last 1 - start_days = 1 - else: - print 'Error for season ' + season_title[season] + ': starting index is month ' + str(time[start_t_season].month) + ', day ' + str(time[start_t_season].day) - return - - # Start accumulating data weighted by days - temp_3d_roms[season,:,:,:] += id.variables['temp'][start_t_season,:,:-15,:]*start_days - salt_3d_roms[season,:,:,:] += id.variables['salt'][start_t_season,:,:-15,:]*start_days - season_days += start_days - - # Between start_t_season and end_t_season, we want all the days - for t in range(start_t_season+1, end_t_season): - temp_3d_roms[season,:,:,:] += id.variables['temp'][t,:,:-15,:]*5 - salt_3d_roms[season,:,:,:] += id.variables['salt'][t,:,:-15,:]*5 - season_days += 5 - - # Figure out how many of the 5 days averaged in end_t_season are - # actually within this season - if time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season]+1: - # Ending day is in position 1 of 5; we care about the first 1 - end_days = 1 - elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season]: - # Ending day is in position 2 of 5; we care about the first 2 - end_days = 2 - elif time[end_t_season].month == end_month[season] and time[end_t_season].day == end_day[season]: - # Ending day is in position 3 of 5; we care about the first 3 - end_days = 3 - elif time[end_t_season].month == end_month[season] and time[end_t_season].day == end_day[season]-1: - # Ending day is in position 4 of 5; we care about the first 4 - end_days = 4 - elif time[end_t_season].month == end_month[season] and time[end_t_season].day == end_day[season]-2: - # Ending day is in position 5 of 5; we care about all 5 - end_days = 5 - else: - print 'Error for season ' + season_title[season] + ': ending index is month ' + str(time[end_t_season].month) + ', day ' + str(time[end_t_season].day) - return - - temp_3d_roms[season,:,:,:] += id.variables['temp'][end_t_season,:,:-15,:]*end_days - salt_3d_roms[season,:,:,:] += id.variables['salt'][end_t_season,:,:-15,:]*end_days - season_days += end_days - - # Check that we got the correct number of days - if season_days != ndays_season[season]: - print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) - return - - # Finished accumulating data, now convert from sum to average - temp_3d_roms[season,:,:,:] /= season_days - salt_3d_roms[season,:,:,:] /= season_days - - # Finished reading data + num_lon = id.variables['lon_rho'].shape[1] + num_lat = id.variables['lat_rho'].shape[0] id.close() + # Calculate seasonal averages of temperature and salinity + print 'Reading temperature' + temp_3d_roms = seasonal_avg_roms(file_path, 'temp', [N, num_lat, num_lon]) + print 'Reading salinity' + salt_3d_roms = seasonal_avg_roms(file_path, 'salt', [N, num_lat, num_lon]) + # Chop off northern boundary + temp_3d_roms = temp_3d_roms[:,:,:-15,:] + salt_3d_roms = salt_3d_roms[:,:,:-15,:] + # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script z_roms_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N) @@ -287,147 +127,6 @@ def temp_salt_seasonal (file_path, lon0, depth_bdry, save=False, fig_name=None): fig.show() -# Linearly interpolate ROMS data, z, and latitude to the specified longitude. -# Input: -# data_3d = array of data, dimension depth x lat x lon -# z_3d = array of depth values (negative, in metres), dimension depth x lat x lon -# lat_2d = array of latitude values, dimension lat x lon -# lon_2d = array of longitude values, dimension lat x lon (between -180 and 180) -# lon0 = longitude to interpolate to (between -180 and 180) -# Output: -# data = array of data interpolated to lon0, dimension depth x lat -# z = array of depth values interpolated to lon0, dimension depth x lat -# lat = array of latitude values interpolated to lon0, dimension depth x lat -def interp_lon_roms (data_3d, z_3d, lat_2d, lon_2d, lon0): - - # Save dimensions - num_depth = size(data_3d, 0) - num_lat = size(data_3d, 1) - num_lon = size(data_3d, 2) - # Set up output arrays - data = ma.empty([num_depth, num_lat]) - z = ma.empty([num_depth, num_lat]) - lat = ma.empty([num_depth, num_lat]) - - # Loop over latitudes; can't find a cleaner way to do this - for j in range(num_lat): - # Extract the longitude values of this slice - lon_tmp = lon_2d[j,:] - # Get indices and coefficients for interpolation - ie, iw, coeffe, coeffw = interp_lon_roms_helper(lon_tmp, lon0) - data[:,j] = coeffe*data_3d[:,j,ie] + coeffw*data_3d[:,j,iw] - z[:,j] = coeffe*z_3d[:,j,ie] + coeffw*z_3d[:,j,iw] - lat[:,j] = coeffe*lat_2d[j,ie] + coeffw*lat_2d[j,iw] - - return data, z, lat - - -# Calculate indices and coefficients for linear interpolation of ROMS longitude. -# This takes care of all the mod 360 nonsense. -# Input: -# lon = 1D array of longitude values (straight out of ROMS i.e. between slightly < 0 and slightly > 360) -# lon0 = longitude to interpolate to (between 0 and 360) -# Output: -# ie, iw, coeffe, coeffw = integers (ie and iw) and coefficients (coeffe and -# coeffw) such that coeffe*lon[ie] + coeffw*lon[iw] = -# lon0, which will also hold for any variable on this -# longitude grid. ie is the index of the nearest point -# to the east of lon0; iw the nearest to the west. -def interp_lon_roms_helper (lon, lon0): - - if lon0 < amin(lon) or lon0 > amax(lon): - # Special case: lon0 on periodic boundary - # Be careful with mod 360 here - - # Find the periodic boundary - dlon = lon[1:] - lon[0:-1] - bdry = argmax(abs(dlon)) - if dlon[bdry] < -300: - # Jumps from almost 360 to just over 0 - iw = bdry - ie = bdry + 1 - else: - # Periodic boundary lines up with the array boundary - iw = size(lon) - 1 - ie = 0 - # Calculate difference between lon0 and lon[iw], mod 360 if necessary - dlon_num = lon0 - lon[iw] - if dlon_num < -300: - dlon_num += 360 - # Calculate difference between lon[ie] and lon[iw], mod 360 - dlon_den = lon[ie] - lon[iw] + 360 - - else: - # General case - - # Add or subtract 360 from longitude values which wrap around - # so that longitude increases monotonically from west to east - i = arange(1, size(lon)+1) - index1 = nonzero((i > 1200)*(lon < 100)) - lon[index1] = lon[index1] + 360 - index2 = nonzero((i < 200)*(lon > 300)) - lon[index2] = lon[index2] - 360 - - # Take mod 360 of lon0 if necessary - if all(lon < lon0): - lon0 -= 360 - if all(lon > lon0): - lon0 += 360 - - # Find the first index eastward of lon0 - ie = nonzero(lon > lon0)[0][0] - # The index before it will be the last index westward of lon0 - iw = ie - 1 - - dlon_num = lon0 - lon[iw] - dlon_den = lon[ie] - lon[iw] - - if dlon_num > 5 or dlon_den > 5: - print 'interp_lon_helper: Problem at periodic boundary' - return - coeff1 = dlon_num/dlon_den - coeff2 = 1 - coeff1 - - return ie, iw, coeff1, coeff2 - - -# Linearly interpolate SOSE data to the specified longitude. -# Input: -# data_3d = arary of data, dimension depth x lat x lon -# lon = 1D array of longitude values (between 0 and 360) -# lon0 = longitude to interpolate to (between 0 and 360) -# Output: -# data = array of data interpolated to lon0, dimension depth x lat -def interp_lon_sose (data_3d, lon, lon0): - - if lon0 < lon[0] or lon0 > lon[-1]: - # Special case: lon0 on periodic boundary - # Be careful with mod 360 here - iw = size(lon)-1 - ie = 0 - # Calculate difference between lon0 and lon[iw], mod 360 if necessary - dlon_num = lon0 - lon[iw] - if dlon_num < -300: - dlon_num += 360 - # Calculate difference between lon[ie] and lon[iw], mod 360 - dlon_den = lon[ie] - lon[iw] + 360 - else: - # General case - # Find the first index eastwards of lon0 - ie = nonzero(lon > lon0)[0][0] - # The index before it will be the last index westward of lon0 - iw = ie - 1 - dlon_num = lon0 - lon[iw] - dlon_den = lon[ie] - lon[iw] - # Coefficients for interpolation - coeff1 = dlon_num/dlon_den - coeff2 = 1 - coeff1 - - # Interpolate - data = coeff1*data_3d[:,:,ie] + coeff2*data_3d[:,:,iw] - return data - - # Command-line interface if __name__ == "__main__": diff --git a/temp_salt_slice.py b/temp_salt_slice.py index 213e291..ae2808e 100644 --- a/temp_salt_slice.py +++ b/temp_salt_slice.py @@ -2,6 +2,7 @@ from numpy import * from matplotlib.pyplot import * from calc_z import * +from interp_lon_roms import * def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=None): @@ -42,8 +43,8 @@ def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=Non if lon0 < 0: lon0 += 360 - temp, z, lat = interp_lon(temp_3d, z_3d, lat_2d, lon_2d, lon0) - salt, z, lat = interp_lon(salt_3d, z_3d, lat_2d, lon_2d, lon0) + temp, z, lat = interp_lon_roms(temp_3d, z_3d, lat_2d, lon_2d, lon0) + salt, z, lat = interp_lon_roms(salt_3d, z_3d, lat_2d, lon_2d, lon0) # Choose latitude bounds based on land mask temp_sum = sum(temp, axis=0) @@ -102,110 +103,6 @@ def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=Non lon0 -= 360 -# Linearly interpolate data, z, and latitude to the specified longitude. -# Input: -# data_3d = array of data, dimension depth x lat x lon -# z_3d = array of depth values (negative, in metres), dimension depth x lat x lon -# lat_2d = array of latitudevalues, dimension lat x lon -# lon_2d = array of longitude values, dimension lat x lon (between -180 and 180) -# lon0 = longitude to interpolate to (between -180 and 180) -# Output: -# data = array of data interpolated to lon0, dimension depth x lat -# z = array of depth values interpolated to lon0, dimension depth x lat -# lat = array of latitude values interpolated to lon0, dimension depth x lat -def interp_lon (data_3d, z_3d, lat_2d, lon_2d, lon0): - - # Save dimensions - num_depth = size(data_3d, 0) - num_lat = size(data_3d, 1) - num_lon = size(data_3d, 2) - # Set up output arrays - data = ma.empty([num_depth, num_lat]) - z = ma.empty([num_depth, num_lat]) - lat = ma.empty([num_depth, num_lat]) - - # Loop over latitudes; can't find a cleaner way to do this - for j in range(num_lat): - # Extract the longitude values of this slice - lon_tmp = lon_2d[j,:] - # Get indices and coefficients for interpolation - ie, iw, coeffe, coeffw = interp_lon_helper(lon_tmp, lon0) - data[:,j] = coeffe*data_3d[:,j,ie] + coeffw*data_3d[:,j,iw] - z[:,j] = coeffe*z_3d[:,j,ie] + coeffw*z_3d[:,j,iw] - lat[:,j] = coeffe*lat_2d[j,ie] + coeffw*lat_2d[j,iw] - - return data, z, lat - - -# Calculate indices and coefficients for linear interpolation of longitude. -# This takes care of all the mod 360 nonsense. -# Input: -# lon = 1D array of longitude values (straight out of ROMS i.e. between slightly < 0 and slightly > 360) -# lon0 = longitude to interpolate to (between 0 and 360) -# Output: -# ie, iw, coeffe, coeffw = integers (ie and iw) and coefficients (coeffe and coeffw) such that -# coeffe*lon[ie] + coeffw*lon[iw] = lon0, which will also hold for any -# variable on this longitude grid. ie is the index of the nearest point -# to the east of lon0; iw the nearest point to the west. -def interp_lon_helper (lon, lon0): - - if lon0 < amin(lon) or lon0 > amax(lon): - # Special case: lon0 on periodic boundary - # Be careful with mod 360 here - - # Find the periodic boundary - dlon = lon[1:] - lon[0:-1] - bdry = argmax(abs(dlon)) - if dlon[bdry] < -300: - # Jumps from almost 360 to just over 0 - iw = bdry - ie = bdry + 1 - else: - # Periodic boundary lines up with the array boundary - iw = size(lon) - 1 - ie = 0 - # Calculate difference between lon0 and lon[iw], mod 360 if necessary - dlon_num = lon0 - lon[iw] - if dlon_num < -300: - dlon_num += 360 - # Calculate difference between lon[ie] and lon[iw], mod 360 - dlon_den = lon[ie] - lon[iw] + 360 - - else: - # General case - - # Add or subtract 360 from longitude values which wrap around - # so that longitude increases monotonically from west to east - i = arange(1, size(lon)+1) - index1 = nonzero((i > 1200)*(lon < 100)) - lon[index1] = lon[index1] + 360 - index2 = nonzero((i < 200)*(lon > 300)) - lon[index2] = lon[index2] - 360 - - # Take mod 360 of lon0 if necessary - if all(lon < lon0): - lon0 -= 360 - if all(lon > lon0): - lon0 += 360 - - # Find the first index eastward of lon0 - ie = nonzero(lon > lon0)[0][0] - # The index before it will be the last index westward of lon0 - iw = ie - 1 - - dlon_num = lon0 - lon[iw] - dlon_den = lon[ie] - lon[iw] - - if dlon_num > 5 or dlon_den > 5: - print 'interp_lon_helper: Problem at periodic boundary' - return - coeff1 = dlon_num/dlon_den - coeff2 = 1 - coeff1 - - return ie, iw, coeff1, coeff2 - - - if __name__ == "__main__": file_path = raw_input("Path to ocean averages file: ") diff --git a/timeseries_dpt.py b/timeseries_dpt.py index 60b1baf..12768ab 100644 --- a/timeseries_dpt.py +++ b/timeseries_dpt.py @@ -3,7 +3,7 @@ from matplotlib.pyplot import * from os.path import * from rotate_vector_roms import * -from zonal_plot import interp_lon_helper +from interp_lon_roms import interp_lon_helper # Calculate and plot timeseries of the Drake Passage transport during a # ROMS simulation. diff --git a/zonal_plot.py b/zonal_plot.py index bacb55b..ddfa3dc 100644 --- a/zonal_plot.py +++ b/zonal_plot.py @@ -3,6 +3,7 @@ from matplotlib.pyplot import * from calc_z import * from rotate_vector_roms import * +from interp_lon_roms import * # Create a zonally averaged or zonally sliced plot (i.e. depth vs latitude) # of the given variable. @@ -110,7 +111,7 @@ def zonal_plot (file_path, var_name, tstep, lon_key, lon0, lon_bounds, depth_min # Interpolate or average data if lon_key == 0: # Interpolate to lon0 - data, z, lat = interp_lon(data_3d, z_3d, lat_2d, lon_2d, lon0) + data, z, lat = interp_lon_roms(data_3d, z_3d, lat_2d, lon_2d, lon0) elif lon_key == 1: # Zonally average over all longitudes # dlon is constant on this grid (0.25 degrees) so this is easy @@ -195,41 +196,6 @@ def zonal_plot (file_path, var_name, tstep, lon_key, lon0, lon_bounds, depth_min lon_bounds[1] -= 360 -# Linearly interpolate data, z, and latitude to the specified longitude. -# Input: -# data_3d = array of data, dimension depth x lat x lon -# z_3d = array of depth values (negative, in metres), dimension depth x lat x lon -# lat_2d = array of latitudevalues, dimension lat x lon -# lon_2d = array of longitude values, dimension lat x lon (between -180 and 180) -# lon0 = longitude to interpolate to (between -180 and 180) -# Output: -# data = array of data interpolated to lon0, dimension depth x lat -# z = array of depth values interpolated to lon0, dimension depth x lat -# lat = array of latitude values interpolated to lon0, dimension depth x lat -def interp_lon (data_3d, z_3d, lat_2d, lon_2d, lon0): - - # Save dimensions - num_depth = size(data_3d, 0) - num_lat = size(data_3d, 1) - num_lon = size(data_3d, 2) - # Set up output arrays - data = ma.empty([num_depth, num_lat]) - z = ma.empty([num_depth, num_lat]) - lat = ma.empty([num_depth, num_lat]) - - # Loop over latitudes; can't find a cleaner way to do this - for j in range(num_lat): - # Extract the longitude values of this slice - lon_tmp = lon_2d[j,:] - # Get indices and coefficients for interpolation - ie, iw, coeffe, coeffw = interp_lon_helper(lon_tmp, lon0) - data[:,j] = coeffe*data_3d[:,j,ie] + coeffw*data_3d[:,j,iw] - z[:,j] = coeffe*z_3d[:,j,ie] + coeffw*z_3d[:,j,iw] - lat[:,j] = coeffe*lat_2d[j,ie] + coeffw*lat_2d[j,iw] - - return data, z, lat - - # Zonally average data, z, and latitude between the specified longitude bounds. # Input: # data_3d = array of data, dimension depth x lat x lon @@ -350,74 +316,6 @@ def average_btw_lons (data_3d, z_3d, lat_2d, lon_2d, lon_bounds): return data, z, lat -# Calculate indices and coefficients for linear interpolation of longitude. -# This takes care of all the mod 360 nonsense. -# Input: -# lon = 1D array of longitude values (straight out of ROMS i.e. between slightly < 0 and slightly > 360) -# lon0 = longitude to interpolate to (between 0 and 360) -# Output: -# ie, iw, coeffe, coeffw = integers (ie and iw) and coefficients (coeffe and coeffw) such that -# coeffe*lon[ie] + coeffw*lon[iw] = lon0, which will also hold for any -# variable on this longitude grid. ie is the index of the nearest point -# to the east of lon0; iw the nearest point to the west. -def interp_lon_helper (lon, lon0): - - if lon0 < amin(lon) or lon0 > amax(lon): - # Special case: lon0 on periodic boundary - # Be careful with mod 360 here - - # Find the periodic boundary - dlon = lon[1:] - lon[0:-1] - bdry = argmax(abs(dlon)) - if dlon[bdry] < -300: - # Jumps from almost 360 to just over 0 - iw = bdry - ie = bdry + 1 - else: - # Periodic boundary lines up with the array boundary - iw = size(lon) - 1 - ie = 0 - # Calculate difference between lon0 and lon[iw], mod 360 if necessary - dlon_num = lon0 - lon[iw] - if dlon_num < -300: - dlon_num += 360 - # Calculate difference between lon[ie] and lon[iw], mod 360 - dlon_den = lon[ie] - lon[iw] + 360 - - else: - # General case - - # Add or subtract 360 from longitude values which wrap around - # so that longitude increases monotonically from west to east - i = arange(1, size(lon)+1) - index1 = nonzero((i > 1200)*(lon < 100)) - lon[index1] = lon[index1] + 360 - index2 = nonzero((i < 200)*(lon > 300)) - lon[index2] = lon[index2] - 360 - - # Take mod 360 of lon0 if necessary - if all(lon < lon0): - lon0 -= 360 - if all(lon > lon0): - lon0 += 360 - - # Find the first index eastward of lon0 - ie = nonzero(lon > lon0)[0][0] - # The index before it will be the last index westward of lon0 - iw = ie - 1 - - dlon_num = lon0 - lon[iw] - dlon_den = lon[ie] - lon[iw] - - if dlon_num > 5 or dlon_den > 5: - print 'interp_lon_helper: Problem at periodic boundary' - return - coeff1 = dlon_num/dlon_den - coeff2 = 1 - coeff1 - - return ie, iw, coeff1, coeff2 - - # Command-line interface if __name__ == "__main__":