esmlab icon indicating copy to clipboard operation
esmlab copied to clipboard

Function to compute true monthly mean from daily mean.

Open dabail10 opened this issue 5 years ago • 3 comments

Thanks for submitting an issue!

Here's a quick checklist in what to include:

  • [x] Include a detailed description of the bug or suggestion

  • [ ] conda list of the conda environment you are using

  • [x] Minimal, self-contained copy-pastable example that generates the issue if possible. Please be concise with code posted. See guidelines below on how to provide a good bug report:

    Bug reports that follow these guidelines are easier to diagnose, and so are often handled much more quickly.

Description

I was trying to fix some bad monthly and daily means from the CESM-CICE history files. @jessluo came up with a script to do a group my year and month. We thought this would be a nice addition to esmlab.

What I Did

The path of the data is:

/gpfs/fs1/p/cesm/pcwg/timeseries-cmip6/b.e21.B1850.f09_g17.CMIP6-piControl.001/ice/proc/tseries/day_1

case = "b.e21.B1850.f09_g17.CMIP6-piControl.001"
year1 = "0061"
year2 = "0099"
year1m = "0001"
year2m = "0099"

filename_month = "../month_1/"+case+".cice.h.aicen."+year1m+"01-"+year2m+"12.nc"
filename_day = case+".cice.h1.aicen_d."+year1+"0101-"+year2+"1231.nc"
filename1 = case+".cice.h1.aicen_d_fix."+year1+"0101-"+year2+"1231.nc"
filename2 = case+".cice.h.aicen_fix."+year1m+"01-"+year2m+"12.nc"

ds = xr.open_dataset(filename1)
ds_day = xr.open_dataset(filename_day)
ds_month = xr.open_dataset(filename_month)

aicen_d = ds['aicen_d']

times = ds_day.time_bounds.values
print(times)

leftbounds_yr = [x[0].timetuple()[0] for x in times]
leftbounds_mo = [x[0].timetuple()[1] for x in times]

yrmo = [cftime.datetime(y, m, 15) for y,m in zip(leftbounds_yr,leftbounds_mo)]
yrmo = np.array(yrmo)

print(yrmo.shape)

aicen_d_val = ds.aicen_d.values

aicen_d = xr.DataArray(aicen_d_val, coords={'time':yrmo, 'NCAT':ds.NCAT, 'TLON':ds.TLON, 'TLAT':ds.TLAT},
                       dims=('time','nc','nj','ni'))

aicen2 = aicen_d.groupby('time').mean(dim='time')

aicen = ds_month['aicen']

aicen.values[720:][:][:][:] = np.where(aicen2.values[:][:][:][:]>1.,1.0e30,aicen2.values[:][:][:][:])

print(aicen.values[0][4][370][180])

aicen.to_netcdf(filename2)
Paste the command(s) you ran and the output.
If there was a crash, please include the traceback here.

dabail10 avatar Feb 28 '19 22:02 dabail10