-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
- Loading branch information
There are no files selected for viewing
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,57 @@ | ||
from netCDF4 import Dataset | ||
import matplotlib.pyplot as plt | ||
import numpy as np | ||
import vor_fast | ||
import vor_fast_setup | ||
import xarray as xr | ||
import sys | ||
|
||
#input variables | ||
plot = False | ||
save2netcdf = True | ||
resolution = 'full' | ||
|
||
# Read in NetCDF file with geopotential height values | ||
print('opening') | ||
ncin = xr.open_dataset('zg_daily_CMAM_CMAM30-SD_r1i1p1_19790101-20101231_10hPa.nc') | ||
gph = ncin.zg.values | ||
lons = ncin.lon.values | ||
lats = ncin.lat.values | ||
days = ncin.time.values#[:100] | ||
ncin.close() | ||
|
||
# Set up cartesian mapping xypoints and restrict to NH | ||
gph_nh, lats_nh, xypoints = vor_fast_setup.setup(gph,lats,lons,'NH') | ||
|
||
|
||
|
||
# Calculate diagnostics for each day | ||
print('Calculating for resolution: '+resolution) | ||
for iday, day in enumerate(days): | ||
if iday % 1000 == 0: | ||
print('Calculating moments for day '+str(iday), day) | ||
moments = vor_fast.calc_moments(gph_nh[iday,:,:],lats_nh,lons,xypoints, | ||
hemisphere='NH',field_type='GPH', | ||
edge=3.02e4,resolution=resolution) | ||
if iday == 0: | ||
ds = xr.Dataset(moments, coords={'time': [day]}) | ||
else: | ||
temp = xr.Dataset(moments, coords={'time': [day]}) | ||
ds = xr.concat([ds,temp], dim='time') | ||
|
||
|
||
if save2netcdf: | ||
print('saving') | ||
ds.to_netcdf('moment_calculation_w_obj-area_CMAM.nc') | ||
|
||
# Plot timeseries | ||
if plot: | ||
print('plotting') | ||
fig, axes = plt.subplots(nrows=2) | ||
ds['aspect_ratio'].plot(ax = axes[0]) | ||
ds['aspect_ratio'].where(ds.aspect_ratio >= 2.4).plot(ax = axes[0], color= 'red') | ||
ds['centroid_latitude'].plot(ax = axes[1]) | ||
ds['centroid_latitude'].where(ds.centroid_latitude < 66).plot(ax = axes[1]) | ||
plt.tight_layout() | ||
plt.savefig('ar_centroid_full.pdf', bbox_inches = 'tight') | ||
|
Large diffs are not rendered by default.
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,37 @@ | ||
import xarray as xr | ||
import glob | ||
from aostools.climate import ComputeRefractiveIndex | ||
|
||
|
||
cesta = '/mnt/4data/CMAM/0A.daily/' | ||
ta_list = sorted(glob.glob(cesta+'ta/ta_6hrPlev_CMAM_CMAM30-SD_r1i1p1_*-*18.nc')) | ||
ua_list = sorted(glob.glob(cesta+'ua/ua_6hrPlev_CMAM_CMAM30-SD_r1i1p1_*-*18.nc')) | ||
Earth_radius = 6.378e+6 | ||
|
||
for i,(ta_file, ua_file) in enumerate(zip(ta_list[:], ua_list[:])): | ||
print(ta_file, ua_file) | ||
suffix = ta_file[len(cesta)+5:] | ||
print(suffix) | ||
ds_ta = xr.open_dataset(ta_file) | ||
ds_ua = xr.open_dataset(ua_file) | ||
|
||
if i==0: | ||
lat = ds_ta.lat.values | ||
pres = ds_ta.plev.values/100. | ||
|
||
uz = ds_ua.ua.mean(['lon'])# | ||
Tz = ds_ta.ta.mean(['lon']) | ||
|
||
n_k_ls = [] | ||
for k in range(1,4): | ||
n_k = ComputeRefractiveIndex(lat,pres,uz,Tz,k) | ||
n_k_ls.append(n_k) | ||
|
||
n_k_xa = xr.concat(n_k_ls, dim='wavenumber') | ||
n_k_xa['wavenumber'] = range(1,4) | ||
|
||
n_k_xa.name = 'refr_index' | ||
outfile = '{0}refr_index/refr_index{1}'.format(cesta, suffix) | ||
print(outfile) | ||
n_k_xa.to_netcdf(outfile) | ||
|