NetCDF_SeasonalMean.py

"""
Created on Thu Oct 17 16:51:54 2013
@author:Sagy Cohen, Uni of Alabama
Calculates long-term mean raster layer (ESRI ASCII) for the WBMsed model output
layers (NetCDF) 
"""
import scipy    
from scipy.io import netcdf
import numpy as np

sim_first=1990
sim_last=2019
interval=(sim_last-sim_first)+1

nutrient = 'Discharge'#'QsConc' #'QxT_WaterTemp'
months = {'Jan': 0, 'Feb': 1,'Mar': 2,'Apr': 3,'May': 4,'Jun': 5,'Jul': 6,'Aug': 7,'Sep': 8,'Oct': 9,'Nov': 10,'Dec': 11}
seasons = {'DJF':[11,0,1], 'MAM':[2,3,4], 'JJA':[5,6,7], 'SON':[8,9,10]}

NetCDF_fn1='/Volumes/Lacie5/Data3/July6_2021_4p3p6_ReRun_SimSpin/SimSpin/WaterDensity/4p3p6_Dist_SimSpinup2+Dist/06min/Monthly/Global_WaterDensity_4p3p6_Dist_SimSpinup2+Dist_06min_mTS'

Output_fn0 = NetCDF_fn1

# set the output ASCII file header
header = ['ncols  3600\n']
header.append('nrows  1500\n')
# header=['ncols  720\n']# % 30min
# header.append('nrows  360\n')# % 30min
header.append('xllcorner -180\n')
header.append('yllcorner -60\n')
header.append('cellsize  0.1\n')
# header.append('cellsize  0.5\n')# %30min
header.append('NODATA_value  -9999\n')


for season in seasons:
    print (season)
    count = 0
    for sim_yr in range(sim_first, sim_last + 1):
        first = (sim_first + count)
        if first > sim_last:
            break
        last = first + interval
        if last > sim_last:
            last = sim_last
            # open and write the header to the output file
        Output_fn = Output_fn0 + '_'+ season + str(first) + '-' + str(last) + '.asc'
        of = open(Output_fn, 'w')
        of.writelines(header)

    for month in seasons[season]:
        print (month)
        for yr in range(first,last+1):
           # print(yr)
            NetCDF_fn=NetCDF_fn1+str(yr)+'.nc'  # set input netcdf name
            f = netcdf.netcdf_file(NetCDF_fn,'r')
            tmp = f.variables[f.subject.decode()][month]
            if yr==first:
                NetCDFQ = [tmp]
            else:
                NetCDFQ.append(tmp)
            f.close

    means=scipy.mean(NetCDFQ,axis=0)
    out=means[:,:]
    #[-9999 if x=='inf' else x for x in out]
    out=out[::-1] #flip
    tmps=(np.isinf(out))*1
    out[tmps==1]=-9999.0
    tmps=(np.isnan(out))*1
    out[tmps==1]=-9999.0
   # pcolor(out, cmap='RdBu', vmin=0, vmax=1000)
    np.savetxt(of,out,fmt='%e')
    #of.writelines(format(out))
    of.close()
    count=count+interval+1
    #xx=f.variables['SedimentFlux'][0,39,-100]