Skip to content
GitLab
Commit cc438207 authored by Klaus Zimmermann's avatar Klaus Zimmermann
Browse files

Cleanup legacy part. 



Signed-off-by: Klaus Zimmermann's avatarKlaus Zimmermann <klaus.zimmermann@smhi.se>
parent ab984913
Hide whitespace changes
Inline Side-by-side
Precip_Dry_TEST.py deleted 100644 → 0
View file @ ab984913
import glob
import numpy as np
import iris
import iris.coords
import iris.coord_categorisation
# from iris.analysis import Aggregator
# from iris.util import rolling_window
from cf_units import Unit
import matplotlib.pyplot as plt
MISSVAL = 1.0e20
# Define the parameters of the index.
threshold_precip = 2.0
var = 'pr'
indir = '/home/rossby/imports/obs/UERRA/SURFEX-MESCAN/input/day/'
inroot = '_EUR-055_UERRA-HARMONIE_RegRean_v1d1-v1d2_SURFEX-MESCAN_v1_day_*.nc' # noqa
filepattern = '{}{}{}'.format(indir, var, inroot)
filelist = sorted(glob.glob(filepattern))
creation_date = ''
tracking_id = ''
DS2short = iris.cube.CubeList([])
DS2long = iris.cube.CubeList([])
ix = np.array([589, 677, 677, 600, 643, 381, 161, 1041]) # y coord
jx = np.array([741, 606, 690, 690, 900, 577, 93, 161]) # x coord
ij = 0
f = [filelist[0]]
datacube = iris.load_cube(f)
datacube *= 24 * 60 * 60
datacube.units = Unit('mm')
# creation_date += datacube.attributes.pop('creation_date', None) + ' '
# tracking_id += datacube.attributes.pop('tracking_id', None) + ' '
iris.coord_categorisation.add_day_of_year(datacube,
'time', name='doy')
iris.coord_categorisation.add_month_number(datacube,
'time', name='month')
# Threshold data
gtThreshold = datacube.copy()
gtThreshold.data = (gtThreshold.data >= threshold_precip) * 1
# pinpoint start and end of spells (wet and dry)
taxis = 0
GTstartend = np.diff(gtThreshold.data, axis=taxis)
GTstartend = np.append(GTstartend, np.zeros(nx, ny), axis=taxis)
# plt.figure(); plt.imshow(np.sum(DSstartend, axis=taxis),origin='lower'); plt.colorbar(); plt.title('DSstartend') # noqa
(ntime, nx, ny) = GTstartend.shape
# handle totally dry gridcells: dry spell start first day and ends last day
z = np.sum(np.abs(GTstartend), axis=taxis) == 0
GTstartend[0, z] = -1
GTstartend[-1, z] = 1
# separate the starts and ends (now it matters whether it is dry or wet)
GT1start = np.argmax((GTstartend > 0) * 1, axis=taxis)
GT1end = np.argmax((GTstartend < 0) * 1, axis=taxis)
z = (GTstartend[0, :, :] == 0) & (GT1start > GT1end)
GTstartend[0, z] = 1
z = (GTstartend[0, :, :] == 0) & (GT1start > GT1end)
GTstartend[0, z] = -1
DS1 =
DSstart = DSstartend < 0
DSend = DSstartend > 0
# plt.figure(); plt.imshow(np.sum(DSstart, axis=taxis) - np.sum(DSend, axis=taxis), origin='lower'); plt.colorbar(); plt.title('DS diff 1') # noqa
plt.figure(); plt.imshow(z * 1, origin='lower'); plt.colorbar(); plt.title('z') # noqa
DSstart[0, z] = 1
plt.figure(); plt.imshow(np.sum(DSstart, axis=taxis) - np.sum(DSend, axis=taxis), origin='lower'); plt.colorbar(); plt.title('DS diff 2') # noqa
DSend[ntime, z] = 1
plt.figure(); plt.imshow(np.sum(DSstart, axis=taxis) - np.sum(DSend, axis=taxis), origin='lower'); plt.colorbar(); plt.title('DS diff 3') # noqa
# find first start and first end of dry spell
DS1end = np.argmax(DSend, axis=taxis)
fixstart = DS1start > DS1end
DSstart[0, fixstart] = 1
plt.figure(); plt.imshow(np.sum(DSstart, axis=taxis) - np.sum(DSend, axis=taxis), origin='lower'); plt.colorbar(); plt.title('DS diff 4') # noqa
DS1start = ntime - np.argmax(np.flip(DSstart, axis=taxis))
DS1end = ntime - np.argmax(np.flip(DSend, axis=taxis))
fixstart = DS1start > DS1end
DSend[ntime, fixstart] = 1
plt.figure(); plt.imshow(np.sum(DSstart, axis=taxis) - np.sum(DSend, axis=taxis), origin='lower'); plt.colorbar(); plt.title('DS diff 5') # noqa
"""
fixstart =
###
mask = np.ma.getmaskarray(Vstart.data)
# iris.util.mask_cube(Vstart, mask)
Vend = gtThreshold.collapsed('time', VEG_BOUNDS,
condition=consecutive_days,
start_end='end')
mask = mask | np.ma.getmaskarray(Vend.data)
Vend.data = np.ma.masked_array(Vend.data, mask)
# iris.util.mask_cube(Vend, mask)
Vlen = Vend - Vstart + 1
Vlen.data[mask] = 0
(Vgdd, testgdd) = veg_gdd(datacube.copy(), Vstart.data, Vend.data)
Vgdd.data[mask] = 0
Vstart.remove_coord('doy')
Vend.remove_coord('doy')
Vlen.remove_coord('doy')
Vgdd.remove_coord('doy')
if test != 'test':
Vstart.remove_coord('year')
Vend.remove_coord('year')
Vlen.remove_coord('year')
Vgdd.remove_coord('year')
vegStart.append(Vstart.copy())
vegEnd.append(Vend.copy())
vegLen.append(Vlen.copy())
vegGDD.append(Vgdd.copy())
if test == 'test':
# co = netcdftime.utime('days since 1961-01-01 00:00:00')
tx365 = np.arange(366)
plt.figure(figsize=(8, 8))
plt.imshow(Vstart.data,
origin='lower',
cmap=plt.get_cmap('coolwarm_r'))
plt.colorbar()
plt.title('Start of Vegetation period')
plt.figure(figsize=(8, 8))
plt.imshow(Vend.data,
origin='lower',
cmap=plt.get_cmap('coolwarm'))
plt.colorbar()
plt.title('End of Vegetation period')
plt.figure(figsize=(8, 8))
plt.imshow(Vlen.data,
origin='lower',
cmap=plt.get_cmap('coolwarm'))
plt.colorbar()
plt.title('Length of Vegetation period')
plt.figure(figsize=(8, 8))
plt.imshow(Vgdd.data,
origin='lower',
cmap=plt.get_cmap('coolwarm'))
plt.colorbar()
plt.title('GDD during the Vegetation period')
for ij in range(len(ix)):
plt.text(jx[ij], ix[ij], '{}'.format(ij), color='black')
for ij in range(len(ix)):
tdata = datacube.data[:, ix[ij], jx[ij]]
tdata = np.insert(tdata, 0, tdata[0])
tgdd = testgdd.data[:, ix[ij], jx[ij]]
tgdd = np.insert(tgdd, 0, tdata[0])
L = tdata < threshold_temperature
tdata_t1 = tdata.copy()
tdata_t1[L] = 0
plt.figure(figsize=(13, 8))
plt.step(tx365, tdata)
# plt.step(tx365, tdata_t1)
plt.hlines(5, 0, 365)
plt.fill_between(tx365, tgdd + threshold_temperature,
threshold_temperature, facecolor='gray',
step='pre')
ss = Vstart.data[ix[ij], jx[ij]] - 1
ee = Vend.data[ix[ij], jx[ij]]
ll = Vlen.data[ix[ij], jx[ij]]
gg = Vgdd.data[ix[ij], jx[ij]]
ylim = plt.ylim()
plt.vlines(ss, ylim[0], ylim[1], color='r')
plt.vlines(ee, ylim[0], ylim[1], color='b')
lat = datacube.coords('latitude')[0].points[ix[ij], jx[ij]]
lon = datacube.coords('longitude')[0].points[ix[ij], jx[ij]]
plt.title(('{} : Vegetation indices '
'[i,j=({},{}), '
'lat,lon=({},{})] : '
'start={}, end={}, len={}, GDD={}'
'').format(ij, ix[ij], jx[ij],
lat, lon,
ss, ee, ll, gg),
fontsize=10)
plt.savefig('./vegIX_example_{}.png'.format(ij+1))
dayVegStart5_ANN = vegStart.merge_cube()
dayVegStart5_ANN.units = Unit('days')
dayVegStart5_ANN._FillValue = MISSVAL
dayVegStart5_ANN.valid_min = 1.
dayVegStart5_ANN.valid_max = 183.
dayVegStart5_ANN.standard_name = None
dayVegStart5_ANN.long_name = u'Start of Vegetation Period (first 4 days above 5 degC)' # noqa
dayVegStart5_ANN.var_name = u'dayVegStart5_ANN'
dayVegStart5_ANN.attributes = datacube.attributes
dayVegStart5_ANN.attributes['creation_date'] = creation_date
dayVegStart5_ANN.attributes['tracking_id'] = tracking_id
dayVegEnd5_ANN = vegEnd.merge_cube()
dayVegEnd5_ANN.units = Unit('days')
dayVegEnd5_ANN._FillValue = MISSVAL
dayVegEnd5_ANN.valid_min = 183.
dayVegEnd5_ANN.valid_max = 366.
dayVegEnd5_ANN.standard_name = None
dayVegEnd5_ANN.long_name = u'End of Vegetation Period (last 4 days above 5 degC)' # noqa
dayVegEnd5_ANN.var_name = u'dayVegEnd5_ANN'
dayVegEnd5_ANN.attributes = datacube.attributes
dayVegEnd5_ANN.attributes['creation_date'] = creation_date
dayVegEnd5_ANN.attributes['tracking_id'] = tracking_id
dayVegLen5_ANN = vegLen.merge_cube()
dayVegLen5_ANN.units = Unit('days')
dayVegLen5_ANN._FillValue = MISSVAL
dayVegLen5_ANN.valid_min = 0.
dayVegLen5_ANN.valid_max = 366.
dayVegLen5_ANN.standard_name = None
dayVegLen5_ANN.long_name = u'Length of Vegetation Period' # noqa
dayVegLen5_ANN.var_name = u'dayVegLen5_ANN'
dayVegLen5_ANN.attributes = datacube.attributes
dayVegLen5_ANN.attributes['creation_date'] = creation_date
dayVegLen5_ANN.attributes['tracking_id'] = tracking_id
GDD5_ANN = vegGDD.merge_cube()
GDD5_ANN.units = Unit('degC days')
GDD5_ANN._FillValue = MISSVAL
GDD5_ANN.valid_min = 0.
GDD5_ANN.valid_max = 15000.
GDD5_ANN.standard_name = None
GDD5_ANN.long_name = u'Growing Degree Days (over 5 degC) during Vegetation Period' # noqa
GDD5_ANN.var_name = u'GDD5_ANN'
GDD5_ANN.attributes = datacube.attributes
GDD5_ANN.attributes['creation_date'] = creation_date
GDD5_ANN.attributes['tracking_id'] = tracking_id
outdir = '/home/sm_lbarr/PROJ/B4EST/'
outroot = '_EUR-055_UERRA-HARMONIE_RegRean_v1d1-v1d2_SURFEX-MESCAN_v1_ann_1961-2015.nc' # noqa
if test != 'test':
iris.coord_categorisation.add_year(dayVegStart5_ANN, 'time', name='year') # noqa
iris.coord_categorisation.add_year(dayVegEnd5_ANN, 'time', name='year')
iris.coord_categorisation.add_year(dayVegLen5_ANN, 'time', name='year')
iris.coord_categorisation.add_year(GDD5_ANN, 'time', name='year')
iris.save(dayVegStart5_ANN, '{}{}{}'.format(outdir, 'dayVegStart5_ANN', outroot), fill_value=MISSVAL) # noqa
iris.save(dayVegEnd5_ANN, '{}{}{}'.format(outdir, 'dayVegEnd5_ANN', outroot), fill_value=MISSVAL) # noqa
iris.save(dayVegLen5_ANN, '{}{}{}'.format(outdir, 'dayVegLen5_ANN', outroot), fill_value=MISSVAL) # noqa
iris.save(GDD5_ANN, '{}{}{}'.format(outdir, 'GDD5_ANN', outroot), fill_value=MISSVAL) # noqa
if test != 'test':
print(dayVegStart5_ANN)
print(dayVegStart5_ANN.coord('year'))
c = dayVegStart5_ANN.extract(iris.Constraint(
coord_values={'year': lambda cell: 1960 < cell < 1991}))
dayVegStart5_6190 = c.collapsed('year', iris.analysis.MEAN)
c = dayVegEnd5_ANN.extract(iris.Constraint(
coord_values={'year': lambda cell: 1960 < cell < 1991}))
dayVegEnd5_6190 = c.collapsed('year', iris.analysis.MEAN)
c = dayVegLen5_ANN.extract(iris.Constraint(
coord_values={'year': lambda cell: 1960 < cell < 1991}))
dayVegLen5_6190 = c.collapsed('year', iris.analysis.MEAN)
c = GDD5_ANN.extract(iris.Constraint(
coord_values={'year': lambda cell: 1960 < cell < 1991}))
GDD5_6190 = c.collapsed('year', iris.analysis.MEAN)
iris.save(dayVegStart5_6190, '{}{}{}'.format(outdir, 'dayVegStart5_6190', outroot), fill_value=MISSVAL) # noqa
iris.save(dayVegEnd5_6190, '{}{}{}'.format(outdir, 'dayVegEnd5_6190', outroot), fill_value=MISSVAL) # noqa
iris.save(dayVegLen5_6190, '{}{}{}'.format(outdir, 'dayVegLen5_6190', outroot), fill_value=MISSVAL) # noqa
iris.save(GDD5_6190, '{}{}{}'.format(outdir, 'GDD5_6190', outroot), fill_value=MISSVAL) # noqa
return (dayVegStart5_ANN, dayVegEnd5_ANN, dayVegLen5_ANN, GDD5_ANN)
if __name__ == "__main__":
main()
if __name__ == "__test__":
main('test')
"""
Tmean_Veg_TEST.py deleted 100644 → 0
View file @ ab984913
import glob
import numpy as np
import iris
import iris.coords
import iris.coord_categorisation
from iris.analysis import Aggregator
from iris.util import rolling_window
from cf_units import Unit
import matplotlib.pyplot as plt
MISSVAL = 1.0e20
def veg_bounds(indata, axis, condition, start_end):
"""
Args:
* indata (array):
thresholded (indicator) data.
* axis (int):
number of the array dimension mapping the time sequences.
(Can also be negative, e.g. '-1' means last dimension)
* condition (int):
minimum number of consecutive days to start the vegetation period.
* start_end (string):
"start" | "end" indicating which bound
"""
se = start_end.lower()
if se.startswith('s'):
data = indata
s = True
elif se.startswith('e'):
data = np.flip(indata, axis=axis)
s = False
else:
raise RuntimeError("Argument start_end not understood <{}>."
"".format(start_end))
if axis < 0:
# just cope with negative axis numbers
axis += indata.ndim
year_length = indata.shape[-1]
# Make an array with data values "windowed" along the time axis.
hit_windows = rolling_window(data, window=condition, axis=axis)
# Find the windows "full of True-s" (along the added 'window axis').
full_windows = np.all(hit_windows, axis=axis+1)
# Find first occurence fulfilling the condition (along the time axis).
doy = np.argmax(full_windows, axis=axis)
# mask gridcells where doy == 0 as a result of the condition not being met
ind = map(np.squeeze,
np.split(np.mgrid[:doy.shape[0], :doy.shape[1]], 2, 0)) + [doy]
mask = (~full_windows[ind]) | (doy > np.ceil(year_length / 2))
masked_doy = np.ma.masked_array(np.float32(doy), mask)
# convert to day of year in the interval 1..year_length
# for end calculations count from end of year
if s:
masked_doy = masked_doy + 1
else:
masked_doy = year_length - masked_doy
return masked_doy
def veg_gdd(data, start, end):
doy = np.arange(data.shape[0]) + 1
dayix = doy[:, None, None]
L = (dayix < start[None, :, :]) | (dayix > end[None, :, :])
data.data[L] = 0.0
data.data = np.maximum(data.data - 5, 0.0)
gdd = data.collapsed('time', iris.analysis.SUM)
return gdd, data
def test():
main('test')
def main(test=''):
# Make an aggregator from the user function.
VEG_BOUNDS = Aggregator('veg_bounds',
veg_bounds,
units_func=lambda units: 1)
# Define the parameters of the index.
threshold_temperature = 5.0
consecutive_days = 4
var = 'tas'
indir = '/home/rossby/imports/obs/UERRA/SURFEX-MESCAN/input/day/'
inroot = '_EUR-055_UERRA-HARMONIE_RegRean_v1d1-v1d2_SURFEX-MESCAN_v1_day_*.nc' # noqa
filepattern = '{}{}{}'.format(indir, var, inroot)
filelist = sorted(glob.glob(filepattern))
creation_date = ''
tracking_id = ''
vegStart = iris.cube.CubeList([])
vegEnd = iris.cube.CubeList([])
vegLen = iris.cube.CubeList([])
vegGDD = iris.cube.CubeList([])
ix = np.array([589, 677, 677, 600, 643, 381, 161, 1041]) # y coord
jx = np.array([741, 606, 690, 690, 900, 577, 93, 161]) # x coord
ij = 0
if test == 'test':
filelist = [filelist[0]]
else:
filelist = filelist[0:2]
for f in filelist: # one file / year at the time ...
datacube = iris.load_cube(f)
creation_date += datacube.attributes.pop('creation_date', None) + ' '
tracking_id += datacube.attributes.pop('tracking_id', None) + ' '
datacube.convert_units(Unit('degC'))
iris.coord_categorisation.add_day_of_year(datacube, 'time', name='doy')
iris.coord_categorisation.add_year(datacube, 'time', name='year')
gtThreshold = datacube.copy()
gtThreshold.data = gtThreshold.data > threshold_temperature
Vstart = gtThreshold.collapsed('time', VEG_BOUNDS,
condition=consecutive_days,
start_end='start')
mask = np.ma.getmaskarray(Vstart.data)
# iris.util.mask_cube(Vstart, mask)
Vend = gtThreshold.collapsed('time', VEG_BOUNDS,
condition=consecutive_days,
start_end='end')
mask = mask | np.ma.getmaskarray(Vend.data)
Vend.data = np.ma.masked_array(Vend.data, mask)
# iris.util.mask_cube(Vend, mask)
Vlen = Vend - Vstart + 1
Vlen.data[mask] = 0
(Vgdd, testgdd) = veg_gdd(datacube.copy(), Vstart.data, Vend.data)
Vgdd.data[mask] = 0
Vstart.remove_coord('doy')
Vend.remove_coord('doy')
Vlen.remove_coord('doy')
Vgdd.remove_coord('doy')
if test != 'test':
Vstart.remove_coord('year')
Vend.remove_coord('year')
Vlen.remove_coord('year')
Vgdd.remove_coord('year')
vegStart.append(Vstart.copy())
vegEnd.append(Vend.copy())
vegLen.append(Vlen.copy())
vegGDD.append(Vgdd.copy())
if test == 'test':
# co = netcdftime.utime('days since 1961-01-01 00:00:00')
tx365 = np.arange(366)
plt.figure(figsize=(8, 8))
plt.imshow(Vstart.data,
origin='lower',
cmap=plt.get_cmap('coolwarm_r'))
plt.colorbar()
plt.title('Start of Vegetation period 1961')
plt.savefig('./vegIX_Start_1961.png')
plt.figure(figsize=(8, 8))
plt.imshow(Vend.data,
origin='lower',
cmap=plt.get_cmap('coolwarm'))
plt.colorbar()
plt.title('End of Vegetation period 1961')
plt.savefig('./vegIX_End_1961.png')
plt.figure(figsize=(8, 8))
plt.imshow(Vlen.data,
origin='lower',
cmap=plt.get_cmap('coolwarm'))
plt.colorbar()
plt.title('Length of Vegetation period 1961')
plt.savefig('./vegIX_Length_1961.png')
plt.figure(figsize=(8, 8))
plt.imshow(Vgdd.data,
origin='lower',
cmap=plt.get_cmap('coolwarm'))
plt.colorbar()
plt.title('GDD during the Vegetation period 1961')
plt.savefig('./vegIX_GDD_1961.png')
for ij in range(len(ix)):
plt.text(jx[ij], ix[ij], '{}'.format(ij + 1), color='black')
for ij in range(len(ix)):
tdata = datacube.data[:, ix[ij], jx[ij]]
tdata = np.insert(tdata, 0, tdata[0])
tgdd = testgdd.data[:, ix[ij], jx[ij]]
tgdd = np.insert(tgdd, 0, tdata[0])
L = tdata < threshold_temperature
tdata_t1 = tdata.copy()
tdata_t1[L] = 0
plt.figure(figsize=(13, 8))
plt.step(tx365, tdata)
# plt.step(tx365, tdata_t1)
plt.hlines(5, 0, 365)
plt.fill_between(tx365, tgdd + threshold_temperature,
threshold_temperature, facecolor='gray',
step='pre')
ss = Vstart.data[ix[ij], jx[ij]] - 1
ee = Vend.data[ix[ij], jx[ij]]
ll = Vlen.data[ix[ij], jx[ij]]
gg = Vgdd.data[ix[ij], jx[ij]]
ylim = plt.ylim()
plt.vlines(ss, ylim[0], ylim[1], color='r')
plt.vlines(ee, ylim[0], ylim[1], color='b')
lat = datacube.coords('latitude')[0].points[ix[ij], jx[ij]]
lon = datacube.coords('longitude')[0].points[ix[ij], jx[ij]]
plt.title(('{} : Vegetation indices '
'[i,j=({},{}), '
'lat,lon=({},{})] : '
'start={}, end={}, len={}, GDD={}'
'').format(ij + 1, ix[ij], jx[ij],
lat, lon,
ss, ee, ll, gg),
fontsize=10)
plt.savefig('./vegIX_example_{}.png'.format(ij+1))
dayVegStart5_ANN = vegStart.merge_cube()
dayVegStart5_ANN.units = Unit('days')
dayVegStart5_ANN._FillValue = MISSVAL
dayVegStart5_ANN.valid_min = 1.
dayVegStart5_ANN.valid_max = 183.
dayVegStart5_ANN.standard_name = None
dayVegStart5_ANN.long_name = u'Start of Vegetation Period (first 4 days above 5 degC)' # noqa
dayVegStart5_ANN.var_name = u'dayVegStart5_ANN'
dayVegStart5_ANN.attributes = datacube.attributes
dayVegStart5_ANN.attributes['creation_date'] = creation_date
dayVegStart5_ANN.attributes['tracking_id'] = tracking_id
dayVegEnd5_ANN = vegEnd.merge_cube()
dayVegEnd5_ANN.units = Unit('days')
dayVegEnd5_ANN._FillValue = MISSVAL
dayVegEnd5_ANN.valid_min = 183.
dayVegEnd5_ANN.valid_max = 366.
dayVegEnd5_ANN.standard_name = None
dayVegEnd5_ANN.long_name = u'End of Vegetation Period (last 4 days above 5 degC)' # noqa
dayVegEnd5_ANN.var_name = u'dayVegEnd5_ANN'
dayVegEnd5_ANN.attributes = datacube.attributes
dayVegEnd5_ANN.attributes['creation_date'] = creation_date
dayVegEnd5_ANN.attributes['tracking_id'] = tracking_id
dayVegLen5_ANN = vegLen.merge_cube()
dayVegLen5_ANN.units = Unit('days')
dayVegLen5_ANN._FillValue = MISSVAL
dayVegLen5_ANN.valid_min = 0.
dayVegLen5_ANN.valid_max = 366.
dayVegLen5_ANN.standard_name = None
dayVegLen5_ANN.long_name = u'Length of Vegetation Period' # noqa
dayVegLen5_ANN.var_name = u'dayVegLen5_ANN'
dayVegLen5_ANN.attributes = datacube.attributes
dayVegLen5_ANN.attributes['creation_date'] = creation_date
dayVegLen5_ANN.attributes['tracking_id'] = tracking_id
GDD5_ANN = vegGDD.merge_cube()
GDD5_ANN.units = Unit('degC days')
GDD5_ANN._FillValue = MISSVAL
GDD5_ANN.valid_min = 0.
GDD5_ANN.valid_max = 15000.
GDD5_ANN.standard_name = None
GDD5_ANN.long_name = u'Growing Degree Days (over 5 degC) during Vegetation Period' # noqa
GDD5_ANN.var_name = u'GDD5_ANN'
GDD5_ANN.attributes = datacube.attributes
GDD5_ANN.attributes['creation_date'] = creation_date
GDD5_ANN.attributes['tracking_id'] = tracking_id
outdir = '/home/sm_lbarr/PROJ/B4EST/'
outroot = '_EUR-055_UERRA-HARMONIE_RegRean_v1d1-v1d2_SURFEX-MESCAN_v1_ann_1961-2015.nc' # noqa
if test != 'test':
iris.coord_categorisation.add_year(dayVegStart5_ANN, 'time', name='year') # noqa
iris.coord_categorisation.add_year(dayVegEnd5_ANN, 'time', name='year')
iris.coord_categorisation.add_year(dayVegLen5_ANN, 'time', name='year')
iris.coord_categorisation.add_year(GDD5_ANN, 'time', name='year')
iris.save(dayVegStart5_ANN, '{}{}{}'.format(outdir, 'dayVegStart5_ANN', outroot), fill_value=MISSVAL) # noqa
iris.save(dayVegEnd5_ANN, '{}{}{}'.format(outdir, 'dayVegEnd5_ANN', outroot), fill_value=MISSVAL) # noqa
iris.save(dayVegLen5_ANN, '{}{}{}'.format(outdir, 'dayVegLen5_ANN', outroot), fill_value=MISSVAL) # noqa
iris.save(GDD5_ANN, '{}{}{}'.format(outdir, 'GDD5_ANN', outroot), fill_value=MISSVAL) # noqa
if test != 'test':
print(dayVegStart5_ANN)
print(dayVegStart5_ANN.coord('year'))
c = dayVegStart5_ANN.extract(iris.Constraint(
coord_values={'year': lambda cell: 1960 < cell < 1991}))
dayVegStart5_6190 = c.collapsed('year', iris.analysis.MEAN)
c = dayVegEnd5_ANN.extract(iris.Constraint(
coord_values={'year': lambda cell: 1960 < cell < 1991}))
dayVegEnd5_6190 = c.collapsed('year', iris.analysis.MEAN)
c = dayVegLen5_ANN.extract(iris.Constraint(
coord_values={'year': lambda cell: 1960 < cell < 1991}))
dayVegLen5_6190 = c.collapsed('year', iris.analysis.MEAN)
c = GDD5_ANN.extract(iris.Constraint(
coord_values={'year': lambda cell: 1960 < cell < 1991}))
GDD5_6190 = c.collapsed('year', iris.analysis.MEAN)
iris.save(dayVegStart5_6190, '{}{}{}'.format(outdir, 'dayVegStart5_6190', outroot), fill_value=MISSVAL) # noqa
iris.save(dayVegEnd5_6190, '{}{}{}'.format(outdir, 'dayVegEnd5_6190', outroot), fill_value=MISSVAL) # noqa
iris.save(dayVegLen5_6190, '{}{}{}'.format(outdir, 'dayVegLen5_6190', outroot), fill_value=MISSVAL) # noqa
iris.save(GDD5_6190, '{}{}{}'.format(outdir, 'GDD5_6190', outroot), fill_value=MISSVAL) # noqa
return (dayVegStart5_ANN, dayVegEnd5_ANN, dayVegLen5_ANN, GDD5_ANN)
if __name__ == "__main__":
main()
if __name__ == "__test__":
main('test')
ani.py deleted 100644 → 0
View file @ ab984913
import numpy as np
import iris
import iris.plot as iplt
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.animation as animation
YEAR0 = 1961
def animate(year, cube):
ycube = cube[year - YEAR0]
hQuad.set_array(ycube.data.ravel())