5975837719

Committed 25 Aug 2023 12:24PM UTC coverage: 85.547% (+0.05%) from 85.502%

Build # 5975837719

Build Type

push

github

Committed by

web-flow

Commit Message

Tiny edit on index page regarding youtube (#1628)

* Update index.rst

* Update index.rst

Run Details

11441 of 13374 relevant lines covered (85.55%)

3.76 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

70.27

/oggm/shop/cru.py

import logging
import warnings

# External libs
import numpy as np
import pandas as pd
import xarray as xr
from scipy import stats

# Optional libs
try:
    import salem
except ImportError:
    pass

# Locals
from oggm import cfg
from oggm import utils
from oggm import entity_task
from oggm.exceptions import MassBalanceCalibrationError, InvalidParamsError

# Module logger
log = logging.getLogger(__name__)

CRU_SERVER = ('https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.04/'
              'cruts.2004151855.v4.04/')

CRU_BASE = 'cru_ts4.04.1901.2019.{}.dat.nc'

CRU_CL = ('https://cluster.klima.uni-bremen.de/~oggm/climate/cru/'
          'cru_cl2.nc.zip')


def set_cru_url(url):
    """If you want to use a different server for CRU (for testing, etc)."""
    global CRU_SERVER
    CRU_SERVER = url


@utils.locked_func
def get_cru_cl_file():
    """Returns the path to the unpacked CRU CL file."""
    return utils.file_extractor(utils.file_downloader(CRU_CL))


@utils.locked_func
def get_cru_file(var=None):
    """Returns a path to the desired CRU baseline climate file.

    If the file is not present, download it.

    Parameters
    ----------
    var : str
        'tmp' for temperature
        'pre' for precipitation

    Returns
    -------
    str
        path to the CRU file
    """

    # Be sure input makes sense
    if var not in ['tmp', 'pre']:
        raise InvalidParamsError('CRU variable {} does not exist!'.format(var))

    # Download
    cru_filename = CRU_BASE.format(var)
    cru_url = CRU_SERVER + '{}/'.format(var) + cru_filename + '.gz'
    return utils.file_extractor(utils.file_downloader(cru_url))


@entity_task(log, writes=['climate_historical'])
def process_cru_data(gdir, tmp_file=None, pre_file=None, y0=None, y1=None,
                     output_filesuffix=None):
    """Processes and writes the CRU baseline climate data for this glacier.

    Interpolates the CRU TS data to the high-resolution CL2 climatologies
    (provided with OGGM) and writes everything to a NetCDF file.

    Parameters
    ----------
    gdir : :py:class:`oggm.GlacierDirectory`
        the glacier directory to process
    tmp_file : str
        path to the CRU temperature file (defaults to the current OGGM chosen
        CRU version)
    pre_file : str
        path to the CRU precip file (defaults to the current OGGM chosen
        CRU version)
    y0 : int
        the starting year of the timeseries to write. The default is to take
        the entire time period available in the file, but with this kwarg
        you can shorten it (to save space or to crop bad data)
    y1 : int
        the end year of the timeseries to write. The default is to take
        the entire time period available in the file, but with this kwarg
        you can shorten it (to save space or to crop bad data)
    output_filesuffix : str
        this add a suffix to the output file (useful to avoid overwriting
        previous experiments)
    """

    if cfg.PARAMS['baseline_climate'] != 'CRU':
        raise InvalidParamsError("cfg.PARAMS['baseline_climate'] should be "
                                 "set to CRU")

    # read the climatology
    ncclim = salem.GeoNetcdf(get_cru_cl_file())

    # and the TS data
    if tmp_file is None:
        tmp_file = get_cru_file('tmp')
    if pre_file is None:
        pre_file = get_cru_file('pre')
    nc_ts_tmp = salem.GeoNetcdf(tmp_file, monthbegin=True)
    nc_ts_pre = salem.GeoNetcdf(pre_file, monthbegin=True)

    # set temporal subset for the ts data (hydro years)
    yrs = nc_ts_pre.time.year
    y0 = yrs[0] if y0 is None else y0
    y1 = yrs[-1] if y1 is None else y1

    nc_ts_tmp.set_period(t0=f'{y0}-01-01', t1=f'{y1}-12-01')
    nc_ts_pre.set_period(t0=f'{y0}-01-01', t1=f'{y1}-12-01')
    time = nc_ts_pre.time
    ny, r = divmod(len(time), 12)
    assert r == 0

    lon = gdir.cenlon
    lat = gdir.cenlat

    # This is guaranteed to work because I prepared the file (I hope)
    ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=1)

    # get climatology data
    loc_hgt = ncclim.get_vardata('elev')
    loc_tmp = ncclim.get_vardata('temp')
    loc_pre = ncclim.get_vardata('prcp')
    loc_lon = ncclim.get_vardata('lon')
    loc_lat = ncclim.get_vardata('lat')

    # see if the center is ok
    if not np.isfinite(loc_hgt[1, 1]):
        # take another candidate where finite
        isok = np.isfinite(loc_hgt)

        # wait: some areas are entirely NaNs, make the subset larger
        _margin = 1
        while not np.any(isok):
            _margin += 1
            ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=_margin)
            loc_hgt = ncclim.get_vardata('elev')
            isok = np.isfinite(loc_hgt)
        if _margin > 1:
            log.debug('(%s) I had to look up for far climate pixels: %s',
                      gdir.rgi_id, _margin)

        # Take the first candidate (doesn't matter which)
        lon, lat = ncclim.grid.ll_coordinates
        lon = lon[isok][0]
        lat = lat[isok][0]
        # Resubset
        ncclim.set_subset()
        ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=1)
        loc_hgt = ncclim.get_vardata('elev')
        loc_tmp = ncclim.get_vardata('temp')
        loc_pre = ncclim.get_vardata('prcp')
        loc_lon = ncclim.get_vardata('lon')
        loc_lat = ncclim.get_vardata('lat')

    assert np.isfinite(loc_hgt[1, 1])
    isok = np.isfinite(loc_hgt)
    hgt_f = loc_hgt[isok].flatten()
    assert len(hgt_f) > 0.

    # maybe this will throw out of bounds warnings
    nc_ts_tmp.set_subset(corners=((lon, lat), (lon, lat)), margin=1)
    nc_ts_pre.set_subset(corners=((lon, lat), (lon, lat)), margin=1)

    # compute monthly anomalies
    # of temp
    ts_tmp = nc_ts_tmp.get_vardata('tmp', as_xarray=True)
    ts_tmp_avg = ts_tmp.sel(time=slice('1961-01-01', '1990-12-01'))
    ts_tmp_avg = ts_tmp_avg.groupby('time.month').mean(dim='time')
    ts_tmp = ts_tmp.groupby('time.month') - ts_tmp_avg
    # of precip
    ts_pre = nc_ts_pre.get_vardata('pre', as_xarray=True)
    ts_pre_avg = ts_pre.sel(time=slice('1961-01-01', '1990-12-01'))
    ts_pre_avg = ts_pre_avg.groupby('time.month').mean(dim='time')
    ts_pre_ano = ts_pre.groupby('time.month') - ts_pre_avg
    # scaled anomalies is the default. Standard anomalies above
    # are used later for where ts_pre_avg == 0
    ts_pre = ts_pre.groupby('time.month') / ts_pre_avg

    # interpolate to HR grid
    if np.any(~np.isfinite(ts_tmp[:, 1, 1])):
        # Extreme case, middle pix is not valid
        # take any valid pix from the 3*3 (and hope there's one)
        found_it = False
        for idi in range(2):
            for idj in range(2):
                if np.all(np.isfinite(ts_tmp[:, idj, idi])):
                    ts_tmp[:, 1, 1] = ts_tmp[:, idj, idi]
                    ts_pre[:, 1, 1] = ts_pre[:, idj, idi]
                    ts_pre_ano[:, 1, 1] = ts_pre_ano[:, idj, idi]
                    found_it = True
        if not found_it:
            msg = '({}) there is no climate data'.format(gdir.rgi_id)
            raise MassBalanceCalibrationError(msg)
    elif np.any(~np.isfinite(ts_tmp)):
        # maybe the side is nan, but we can do nearest
        ts_tmp = ncclim.grid.map_gridded_data(ts_tmp.values, nc_ts_tmp.grid,
                                              interp='nearest')
        ts_pre = ncclim.grid.map_gridded_data(ts_pre.values, nc_ts_pre.grid,
                                              interp='nearest')
        ts_pre_ano = ncclim.grid.map_gridded_data(ts_pre_ano.values,
                                                  nc_ts_pre.grid,
                                                  interp='nearest')
    else:
        # We can do bilinear
        ts_tmp = ncclim.grid.map_gridded_data(ts_tmp.values, nc_ts_tmp.grid,
                                              interp='linear')
        ts_pre = ncclim.grid.map_gridded_data(ts_pre.values, nc_ts_pre.grid,
                                              interp='linear')
        ts_pre_ano = ncclim.grid.map_gridded_data(ts_pre_ano.values,
                                                  nc_ts_pre.grid,
                                                  interp='linear')

    # take the center pixel and add it to the CRU CL clim
    # for temp
    loc_tmp = xr.DataArray(loc_tmp[:, 1, 1], dims=['month'],
                           coords={'month': ts_tmp_avg.month})
    ts_tmp = xr.DataArray(ts_tmp[:, 1, 1], dims=['time'],
                          coords={'time': time})
    ts_tmp = ts_tmp.groupby('time.month') + loc_tmp
    # for prcp
    loc_pre = xr.DataArray(loc_pre[:, 1, 1], dims=['month'],
                           coords={'month': ts_pre_avg.month})
    ts_pre = xr.DataArray(ts_pre[:, 1, 1], dims=['time'],
                          coords={'time': time})
    ts_pre_ano = xr.DataArray(ts_pre_ano[:, 1, 1], dims=['time'],
                              coords={'time': time})
    # scaled anomalies
    ts_pre = ts_pre.groupby('time.month') * loc_pre
    # standard anomalies
    ts_pre_ano = ts_pre_ano.groupby('time.month') + loc_pre
    # Correct infinite values with standard anomalies
    ts_pre.values = np.where(np.isfinite(ts_pre.values),
                             ts_pre.values,
                             ts_pre_ano.values)
    # The last step might create negative values (unlikely). Clip them
    ts_pre.values = utils.clip_min(ts_pre.values, 0)

    # done
    loc_hgt = loc_hgt[1, 1]
    loc_lon = loc_lon[1]
    loc_lat = loc_lat[1]
    assert np.isfinite(loc_hgt)
    assert np.all(np.isfinite(ts_pre.values))
    assert np.all(np.isfinite(ts_tmp.values))

    gdir.write_monthly_climate_file(time, ts_pre.values, ts_tmp.values,
                                    loc_hgt, loc_lon, loc_lat,
                                    filesuffix=output_filesuffix,
                                    source=nc_ts_tmp._nc.title[:10])

    ncclim._nc.close()
    nc_ts_tmp._nc.close()
    nc_ts_pre._nc.close()


@entity_task(log, writes=['climate_historical'])
def process_dummy_cru_file(gdir, sigma_temp=2, sigma_prcp=0.5, seed=None,
                           y0=None, y1=None, output_filesuffix=None):
    """Create a simple baseline climate file for this glacier - for testing!

    This simply reproduces the climatology with a little randomness in it.

    TODO: extend the functionality by allowing a monthly varying sigma

    Parameters
    ----------
    gdir : GlacierDirectory
        the glacier directory
    sigma_temp : float
        the standard deviation of the random timeseries (set to 0 for constant
        ts)
    sigma_prcp : float
        the standard deviation of the random timeseries (set to 0 for constant
        ts)
    seed : int
        the RandomState seed
    y0 : int
        the starting year of the timeseries to write. The default is to take
        the entire time period available in the file, but with this kwarg
        you can shorten it (to save space or to crop bad data)
    y1 : int
        the starting year of the timeseries to write. The default is to take
        the entire time period available in the file, but with this kwarg
        you can shorten it (to save space or to crop bad data)
    output_filesuffix : str
        this add a suffix to the output file (useful to avoid overwriting
        previous experiments)
    """

    # read the climatology
    clfile = get_cru_cl_file()
    ncclim = salem.GeoNetcdf(clfile)

    y0 = 1901 if y0 is None else y0
    y1 = 2021 if y1 is None else y1

    time = pd.date_range(start=f'{y0}-01-01', end=f'{y1}-12-01', freq='MS')
    ny, r = divmod(len(time), 12)
    assert r == 0

    lon = gdir.cenlon
    lat = gdir.cenlat

    # This is guaranteed to work because I prepared the file (I hope)
    ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=1)

    # get climatology data
    loc_hgt = ncclim.get_vardata('elev')
    loc_tmp = ncclim.get_vardata('temp')
    loc_pre = ncclim.get_vardata('prcp')
    loc_lon = ncclim.get_vardata('lon')
    loc_lat = ncclim.get_vardata('lat')

    # see if the center is ok
    if not np.isfinite(loc_hgt[1, 1]):
        # take another candidate where finite
        isok = np.isfinite(loc_hgt)

        # wait: some areas are entirely NaNs, make the subset larger
        _margin = 1
        while not np.any(isok):
            _margin += 1
            ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=_margin)
            loc_hgt = ncclim.get_vardata('elev')
            isok = np.isfinite(loc_hgt)
        if _margin > 1:
            log.debug('(%s) I had to look up for far climate pixels: %s',
                      gdir.rgi_id, _margin)

        # Take the first candidate (doesn't matter which)
        lon, lat = ncclim.grid.ll_coordinates
        lon = lon[isok][0]
        lat = lat[isok][0]
        # Resubset
        ncclim.set_subset()
        ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=1)
        loc_hgt = ncclim.get_vardata('elev')
        loc_tmp = ncclim.get_vardata('temp')
        loc_pre = ncclim.get_vardata('prcp')
        loc_lon = ncclim.get_vardata('lon')
        loc_lat = ncclim.get_vardata('lat')

    assert np.isfinite(loc_hgt[1, 1])
    isok = np.isfinite(loc_hgt)
    hgt_f = loc_hgt[isok].flatten()
    assert len(hgt_f) > 0.


    # Make DataArrays
    rng = np.random.RandomState(seed)
    loc_tmp = xr.DataArray(loc_tmp[:, 1, 1], dims=['month'],
                           coords={'month': np.arange(1, 13)})
    ts_tmp = rng.randn(len(time)) * sigma_temp
    ts_tmp = xr.DataArray(ts_tmp, dims=['time'],
                          coords={'time': time})

    loc_pre = xr.DataArray(loc_pre[:, 1, 1], dims=['month'],
                           coords={'month': np.arange(1, 13)})
    ts_pre = utils.clip_min(rng.randn(len(time)) * sigma_prcp + 1, 0)
    ts_pre = xr.DataArray(ts_pre, dims=['time'],
                          coords={'time': time})

    # Create the time series
    ts_tmp = ts_tmp.groupby('time.month') + loc_tmp
    ts_pre = ts_pre.groupby('time.month') * loc_pre

    # done
    loc_hgt = loc_hgt[1, 1]
    loc_lon = loc_lon[1]
    loc_lat = loc_lat[1]
    assert np.isfinite(loc_hgt)

    gdir.write_monthly_climate_file(time, ts_pre.values, ts_tmp.values,
                                    loc_hgt, loc_lon, loc_lat,
                                    filesuffix=output_filesuffix,
                                    source='CRU CL2 and some randomness')
    ncclim._nc.close()

1	import logging	9✔
2	import warnings	9✔
3
4	# External libs
5	import numpy as np	9✔
6	import pandas as pd	9✔
7	import xarray as xr	9✔
8	from scipy import stats	9✔
9
10	# Optional libs
11	try:	9✔
12	import salem	9✔
13	except ImportError:
14	pass
15
16	# Locals
17	from oggm import cfg	9✔
18	from oggm import utils	9✔
19	from oggm import entity_task	9✔
20	from oggm.exceptions import MassBalanceCalibrationError, InvalidParamsError	9✔
21
22	# Module logger
23	log = logging.getLogger(__name__)	9✔
24
25	CRU_SERVER = ('https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.04/'	9✔
26	'cruts.2004151855.v4.04/')
27
28	CRU_BASE = 'cru_ts4.04.1901.2019.{}.dat.nc'	9✔
29
30	CRU_CL = ('https://cluster.klima.uni-bremen.de/~oggm/climate/cru/'	9✔
31	'cru_cl2.nc.zip')
32
33
34	def set_cru_url(url):	9✔
35	"""If you want to use a different server for CRU (for testing, etc)."""
36	global CRU_SERVER
37	CRU_SERVER = url	×
38
39
40	@utils.locked_func	9✔
41	def get_cru_cl_file():	9✔
42	"""Returns the path to the unpacked CRU CL file."""
43	return utils.file_extractor(utils.file_downloader(CRU_CL))	6✔
44
45
46	@utils.locked_func	9✔
47	def get_cru_file(var=None):	9✔
48	"""Returns a path to the desired CRU baseline climate file.
49
50	If the file is not present, download it.
51
52	Parameters
53	----------
54	var : str
55	'tmp' for temperature
56	'pre' for precipitation
57
58	Returns
59	-------
60	str
61	path to the CRU file
62	"""
63
64	# Be sure input makes sense
65	if var not in ['tmp', 'pre']:	2✔
66	raise InvalidParamsError('CRU variable {} does not exist!'.format(var))	×
67
68	# Download
69	cru_filename = CRU_BASE.format(var)	2✔
70	cru_url = CRU_SERVER + '{}/'.format(var) + cru_filename + '.gz'	2✔
71	return utils.file_extractor(utils.file_downloader(cru_url))	2✔
72
73
74	@entity_task(log, writes=['climate_historical'])	9✔
75	def process_cru_data(gdir, tmp_file=None, pre_file=None, y0=None, y1=None,	9✔
76	output_filesuffix=None):
77	"""Processes and writes the CRU baseline climate data for this glacier.
78
79	Interpolates the CRU TS data to the high-resolution CL2 climatologies
80	(provided with OGGM) and writes everything to a NetCDF file.
81
82	Parameters
83	----------
84	gdir : :py:class:`oggm.GlacierDirectory`
85	the glacier directory to process
86	tmp_file : str
87	path to the CRU temperature file (defaults to the current OGGM chosen
88	CRU version)
89	pre_file : str
90	path to the CRU precip file (defaults to the current OGGM chosen
91	CRU version)
92	y0 : int
93	the starting year of the timeseries to write. The default is to take
94	the entire time period available in the file, but with this kwarg
95	you can shorten it (to save space or to crop bad data)
96	y1 : int
97	the end year of the timeseries to write. The default is to take
98	the entire time period available in the file, but with this kwarg
99	you can shorten it (to save space or to crop bad data)
100	output_filesuffix : str
101	this add a suffix to the output file (useful to avoid overwriting
102	previous experiments)
103	"""
104
105	if cfg.PARAMS['baseline_climate'] != 'CRU':	3✔
106	raise InvalidParamsError("cfg.PARAMS['baseline_climate'] should be "	×
107	"set to CRU")
108
109	# read the climatology
110	ncclim = salem.GeoNetcdf(get_cru_cl_file())	3✔
111
112	# and the TS data
113	if tmp_file is None:	3✔
114	tmp_file = get_cru_file('tmp')	2✔
115	if pre_file is None:	3✔
116	pre_file = get_cru_file('pre')	2✔
117	nc_ts_tmp = salem.GeoNetcdf(tmp_file, monthbegin=True)	3✔
118	nc_ts_pre = salem.GeoNetcdf(pre_file, monthbegin=True)	3✔
119
120	# set temporal subset for the ts data (hydro years)
121	yrs = nc_ts_pre.time.year	3✔
122	y0 = yrs[0] if y0 is None else y0	3✔
123	y1 = yrs[-1] if y1 is None else y1	3✔
124
125	nc_ts_tmp.set_period(t0=f'{y0}-01-01', t1=f'{y1}-12-01')	3✔
126	nc_ts_pre.set_period(t0=f'{y0}-01-01', t1=f'{y1}-12-01')	3✔
127	time = nc_ts_pre.time	3✔
128	ny, r = divmod(len(time), 12)	3✔
129	assert r == 0	3✔
130
131	lon = gdir.cenlon	3✔
132	lat = gdir.cenlat	3✔
133
134	# This is guaranteed to work because I prepared the file (I hope)
135	ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=1)	3✔
136
137	# get climatology data
138	loc_hgt = ncclim.get_vardata('elev')	3✔
139	loc_tmp = ncclim.get_vardata('temp')	3✔
140	loc_pre = ncclim.get_vardata('prcp')	3✔
141	loc_lon = ncclim.get_vardata('lon')	3✔
142	loc_lat = ncclim.get_vardata('lat')	3✔
143
144	# see if the center is ok
145	if not np.isfinite(loc_hgt[1, 1]):	3✔
146	# take another candidate where finite
147	isok = np.isfinite(loc_hgt)	×
148
149	# wait: some areas are entirely NaNs, make the subset larger
150	_margin = 1	×
151	while not np.any(isok):	×
152	_margin += 1	×
153	ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=_margin)	×
154	loc_hgt = ncclim.get_vardata('elev')	×
155	isok = np.isfinite(loc_hgt)	×
156	if _margin > 1:	×
157	log.debug('(%s) I had to look up for far climate pixels: %s',	×
158	gdir.rgi_id, _margin)
159
160	# Take the first candidate (doesn't matter which)
161	lon, lat = ncclim.grid.ll_coordinates	×
162	lon = lon[isok][0]	×
163	lat = lat[isok][0]	×
164	# Resubset
165	ncclim.set_subset()	×
166	ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=1)	×
167	loc_hgt = ncclim.get_vardata('elev')	×
168	loc_tmp = ncclim.get_vardata('temp')	×
169	loc_pre = ncclim.get_vardata('prcp')	×
170	loc_lon = ncclim.get_vardata('lon')	×
171	loc_lat = ncclim.get_vardata('lat')	×
172
173	assert np.isfinite(loc_hgt[1, 1])	3✔
174	isok = np.isfinite(loc_hgt)	3✔
175	hgt_f = loc_hgt[isok].flatten()	3✔
176	assert len(hgt_f) > 0.	3✔
177
178	# maybe this will throw out of bounds warnings
179	nc_ts_tmp.set_subset(corners=((lon, lat), (lon, lat)), margin=1)	3✔
180	nc_ts_pre.set_subset(corners=((lon, lat), (lon, lat)), margin=1)	3✔
181
182	# compute monthly anomalies
183	# of temp
184	ts_tmp = nc_ts_tmp.get_vardata('tmp', as_xarray=True)	3✔
185	ts_tmp_avg = ts_tmp.sel(time=slice('1961-01-01', '1990-12-01'))	3✔
186	ts_tmp_avg = ts_tmp_avg.groupby('time.month').mean(dim='time')	3✔
187	ts_tmp = ts_tmp.groupby('time.month') - ts_tmp_avg	3✔
188	# of precip
189	ts_pre = nc_ts_pre.get_vardata('pre', as_xarray=True)	3✔
190	ts_pre_avg = ts_pre.sel(time=slice('1961-01-01', '1990-12-01'))	3✔
191	ts_pre_avg = ts_pre_avg.groupby('time.month').mean(dim='time')	3✔
192	ts_pre_ano = ts_pre.groupby('time.month') - ts_pre_avg	3✔
193	# scaled anomalies is the default. Standard anomalies above
194	# are used later for where ts_pre_avg == 0
195	ts_pre = ts_pre.groupby('time.month') / ts_pre_avg	3✔
196
197	# interpolate to HR grid
198	if np.any(~np.isfinite(ts_tmp[:, 1, 1])):	3✔
199	# Extreme case, middle pix is not valid
200	# take any valid pix from the 3*3 (and hope there's one)
201	found_it = False	×
202	for idi in range(2):	×
203	for idj in range(2):	×
204	if np.all(np.isfinite(ts_tmp[:, idj, idi])):	×
205	ts_tmp[:, 1, 1] = ts_tmp[:, idj, idi]	×
206	ts_pre[:, 1, 1] = ts_pre[:, idj, idi]	×
207	ts_pre_ano[:, 1, 1] = ts_pre_ano[:, idj, idi]	×
208	found_it = True	×
209	if not found_it:	×
210	msg = '({}) there is no climate data'.format(gdir.rgi_id)	×
211	raise MassBalanceCalibrationError(msg)	×
212	elif np.any(~np.isfinite(ts_tmp)):	3✔
213	# maybe the side is nan, but we can do nearest
214	ts_tmp = ncclim.grid.map_gridded_data(ts_tmp.values, nc_ts_tmp.grid,	×
215	interp='nearest')
216	ts_pre = ncclim.grid.map_gridded_data(ts_pre.values, nc_ts_pre.grid,	×
217	interp='nearest')
218	ts_pre_ano = ncclim.grid.map_gridded_data(ts_pre_ano.values,	×
219	nc_ts_pre.grid,
220	interp='nearest')
221	else:
222	# We can do bilinear
223	ts_tmp = ncclim.grid.map_gridded_data(ts_tmp.values, nc_ts_tmp.grid,	3✔
224	interp='linear')
225	ts_pre = ncclim.grid.map_gridded_data(ts_pre.values, nc_ts_pre.grid,	3✔
226	interp='linear')
227	ts_pre_ano = ncclim.grid.map_gridded_data(ts_pre_ano.values,	3✔
228	nc_ts_pre.grid,
229	interp='linear')
230
231	# take the center pixel and add it to the CRU CL clim
232	# for temp
233	loc_tmp = xr.DataArray(loc_tmp[:, 1, 1], dims=['month'],	3✔
234	coords={'month': ts_tmp_avg.month})
235	ts_tmp = xr.DataArray(ts_tmp[:, 1, 1], dims=['time'],	3✔
236	coords={'time': time})
237	ts_tmp = ts_tmp.groupby('time.month') + loc_tmp	3✔
238	# for prcp
239	loc_pre = xr.DataArray(loc_pre[:, 1, 1], dims=['month'],	3✔
240	coords={'month': ts_pre_avg.month})
241	ts_pre = xr.DataArray(ts_pre[:, 1, 1], dims=['time'],	3✔
242	coords={'time': time})
243	ts_pre_ano = xr.DataArray(ts_pre_ano[:, 1, 1], dims=['time'],	3✔
244	coords={'time': time})
245	# scaled anomalies
246	ts_pre = ts_pre.groupby('time.month') * loc_pre	3✔
247	# standard anomalies
248	ts_pre_ano = ts_pre_ano.groupby('time.month') + loc_pre	3✔
249	# Correct infinite values with standard anomalies
250	ts_pre.values = np.where(np.isfinite(ts_pre.values),	3✔
251	ts_pre.values,
252	ts_pre_ano.values)
253	# The last step might create negative values (unlikely). Clip them
254	ts_pre.values = utils.clip_min(ts_pre.values, 0)	3✔
255
256	# done
257	loc_hgt = loc_hgt[1, 1]	3✔
258	loc_lon = loc_lon[1]	3✔
259	loc_lat = loc_lat[1]	3✔
260	assert np.isfinite(loc_hgt)	3✔
261	assert np.all(np.isfinite(ts_pre.values))	3✔
262	assert np.all(np.isfinite(ts_tmp.values))	3✔
263
264	gdir.write_monthly_climate_file(time, ts_pre.values, ts_tmp.values,	3✔
265	loc_hgt, loc_lon, loc_lat,
266	filesuffix=output_filesuffix,
267	source=nc_ts_tmp._nc.title[:10])
268
269	ncclim._nc.close()	3✔
270	nc_ts_tmp._nc.close()	3✔
271	nc_ts_pre._nc.close()	3✔
272
273
274	@entity_task(log, writes=['climate_historical'])	9✔
275	def process_dummy_cru_file(gdir, sigma_temp=2, sigma_prcp=0.5, seed=None,	9✔
276	y0=None, y1=None, output_filesuffix=None):
277	"""Create a simple baseline climate file for this glacier - for testing!
278
279	This simply reproduces the climatology with a little randomness in it.
280
281	TODO: extend the functionality by allowing a monthly varying sigma
282
283	Parameters
284	----------
285	gdir : GlacierDirectory
286	the glacier directory
287	sigma_temp : float
288	the standard deviation of the random timeseries (set to 0 for constant
289	ts)
290	sigma_prcp : float
291	the standard deviation of the random timeseries (set to 0 for constant
292	ts)
293	seed : int
294	the RandomState seed
295	y0 : int
296	the starting year of the timeseries to write. The default is to take
297	the entire time period available in the file, but with this kwarg
298	you can shorten it (to save space or to crop bad data)
299	y1 : int
300	the starting year of the timeseries to write. The default is to take
301	the entire time period available in the file, but with this kwarg
302	you can shorten it (to save space or to crop bad data)
303	output_filesuffix : str
304	this add a suffix to the output file (useful to avoid overwriting
305	previous experiments)
306	"""
307
308	# read the climatology
309	clfile = get_cru_cl_file()	3✔
310	ncclim = salem.GeoNetcdf(clfile)	3✔
311
312	y0 = 1901 if y0 is None else y0	3✔
313	y1 = 2021 if y1 is None else y1	3✔
314
315	time = pd.date_range(start=f'{y0}-01-01', end=f'{y1}-12-01', freq='MS')	3✔
316	ny, r = divmod(len(time), 12)	3✔
317	assert r == 0	3✔
318
319	lon = gdir.cenlon	3✔
320	lat = gdir.cenlat	3✔
321
322	# This is guaranteed to work because I prepared the file (I hope)
323	ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=1)	3✔
324
325	# get climatology data
326	loc_hgt = ncclim.get_vardata('elev')	3✔
327	loc_tmp = ncclim.get_vardata('temp')	3✔
328	loc_pre = ncclim.get_vardata('prcp')	3✔
329	loc_lon = ncclim.get_vardata('lon')	3✔
330	loc_lat = ncclim.get_vardata('lat')	3✔
331
332	# see if the center is ok
333	if not np.isfinite(loc_hgt[1, 1]):	3✔
334	# take another candidate where finite
335	isok = np.isfinite(loc_hgt)	×
336
337	# wait: some areas are entirely NaNs, make the subset larger
338	_margin = 1	×
339	while not np.any(isok):	×
340	_margin += 1	×
341	ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=_margin)	×
342	loc_hgt = ncclim.get_vardata('elev')	×
343	isok = np.isfinite(loc_hgt)	×
344	if _margin > 1:	×
345	log.debug('(%s) I had to look up for far climate pixels: %s',	×
346	gdir.rgi_id, _margin)
347
348	# Take the first candidate (doesn't matter which)
349	lon, lat = ncclim.grid.ll_coordinates	×
350	lon = lon[isok][0]	×
351	lat = lat[isok][0]	×
352	# Resubset
353	ncclim.set_subset()	×
354	ncclim.set_subset(corners=((lon, lat), (lon, lat)), margin=1)	×
355	loc_hgt = ncclim.get_vardata('elev')	×
356	loc_tmp = ncclim.get_vardata('temp')	×
357	loc_pre = ncclim.get_vardata('prcp')	×
358	loc_lon = ncclim.get_vardata('lon')	×
359	loc_lat = ncclim.get_vardata('lat')	×
360
361	assert np.isfinite(loc_hgt[1, 1])	3✔
362	isok = np.isfinite(loc_hgt)	3✔
363	hgt_f = loc_hgt[isok].flatten()	3✔
364	assert len(hgt_f) > 0.	3✔
365
366
367	# Make DataArrays
368	rng = np.random.RandomState(seed)	3✔
369	loc_tmp = xr.DataArray(loc_tmp[:, 1, 1], dims=['month'],	3✔
370	coords={'month': np.arange(1, 13)})
371	ts_tmp = rng.randn(len(time)) * sigma_temp	3✔
372	ts_tmp = xr.DataArray(ts_tmp, dims=['time'],	3✔
373	coords={'time': time})
374
375	loc_pre = xr.DataArray(loc_pre[:, 1, 1], dims=['month'],	3✔
376	coords={'month': np.arange(1, 13)})
377	ts_pre = utils.clip_min(rng.randn(len(time)) * sigma_prcp + 1, 0)	3✔
378	ts_pre = xr.DataArray(ts_pre, dims=['time'],	3✔
379	coords={'time': time})
380
381	# Create the time series
382	ts_tmp = ts_tmp.groupby('time.month') + loc_tmp	3✔
383	ts_pre = ts_pre.groupby('time.month') * loc_pre	3✔
384
385	# done
386	loc_hgt = loc_hgt[1, 1]	3✔
387	loc_lon = loc_lon[1]	3✔
388	loc_lat = loc_lat[1]	3✔
389	assert np.isfinite(loc_hgt)	3✔
390
391	gdir.write_monthly_climate_file(time, ts_pre.values, ts_tmp.values,	3✔
392	loc_hgt, loc_lon, loc_lat,
393	filesuffix=output_filesuffix,
394	source='CRU CL2 and some randomness')
395	ncclim._nc.close()	3✔

OGGM / oggm / 5975837719

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous