Preprocess forcing for esemble simulation
Change forcing RCA3
import numpy as np
import xarray as xr
import pylab as plt
import multiprocessing as mp
import sys
import argparse
import os
from scipy.optimize import curve_fit
import glob
from cdo import Cdo
from subprocess import call
cdo = Cdo()
import warnings
warnings.filterwarnings("ignore")
Find fit for MCA and LIA in RCA3
ds = xr.open_dataset("/gfs1/work/mvkfbiow/postprocessing_baltic_past" +
"/data/rca_temp_full.nc")
temp1d = ds.temp.mean(["lon", "lat"])
temp1d_m = temp1d.rolling(time=365, min_periods=1).mean()
x = np.arange(len(temp1d))
def do_fit(x,temp1d_m):
"""
Polyfit function
"""
fit_params = np.polyfit(x, temp1d_m.values, 3)
z = np.poly1d(fit_params)
predict = z(x)
return predict
predict = do_fit(np.arange(0, len(temp1d_m)), temp1d_m)
polynom_fit = xr.DataArray(predict, coords=[temp1d.time])
def do_sin_fit(x, freq, amplitude):
return np.sin(x * freq ) * amplitude + 275.50742
popt, pcov = curve_fit(do_sin_fit, x, predict, bounds=(1/3000, [1/1900, 0.4]))
original_sin_fit = do_sin_fit(x, *popt)
original_sin_fit = xr.DataArray(original_sin_fit, coords=[temp1d.time])
amplitudes = [0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
plt.figure(figsize=(12, 8))
ax = plt.subplot(111)
ax.spines["top"].set_visible(False)
ax.spines["bottom"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["left"].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.yaxis.grid()
temp1d.rolling(time=365).mean().plot(ax=ax)
for i in amplitudes:
original_sin_fit = do_sin_fit(x, popt[0], i)
original_sin_fit = xr.DataArray(original_sin_fit, coords=[temp1d.time])
original_sin_fit.plot(ax = ax, alpha = 0.9, color = "green", label="Sinus fit {}".format(i), )
polynom_fit.plot(ax = ax, alpha = 0.3, color = "red", label="Polynom", )
ax.legend()
plt.title("Spatial Averaged Temperature", fontsize=17, ha="center")
plt.text(-0.1, -0.2,"Data: Regional Climate Model (RCA3), Schimanke et. al, 2012",
horizontalalignment='left',
verticalalignment='center',
transform = ax.transAxes)
Text(-0.1, -0.2, 'Data: Regional Climate Model (RCA3), Schimanke et. al, 2012')
Construct new data with fit
ds3 = xr.Dataset({'data': (('time'), temp1d.values)},
{'time': temp1d_m.time})
for i in amplitudes:
original_sin_fit = do_sin_fit(x, popt[0], i)
original_sin_fit = xr.DataArray(original_sin_fit, coords=[temp1d.time])
var = temp1d - temp1d.mean() # get variance in data
new = var + original_sin_fit # estimate new data with variance and fit
diffs = new - temp1d # diff between new and old
new = new.where(diffs > 0, other=temp1d)
diffs = new - temp1d
ds3['new_data_{}'.format(i)] = new
ds3['diffs_{}'.format(i)] = diffs
plt.figure(figsize=(12, 8))
ax = plt.subplot(211)
ax.spines["top"].set_visible(False)
ax.spines["bottom"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["left"].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.yaxis.grid()
jet= plt.get_cmap('Spectral_r')
colors = iter(jet(np.linspace(0,1,6)))
ds3.data.rolling(time=365).mean().plot(ax = ax, alpha = 0.3, color = "blue", label="Original data", )
for i in amplitudes:
ds3['new_data_{}'.format(i)].rolling(time=365).mean().plot(ax = ax, alpha = 0.8, color = next(colors), label="New data {}".format(i))
ax.legend()
plt.title("New temperature Data for MCA", fontsize=17, ha="center")
ax = plt.subplot(212)
ax.spines["top"].set_visible(False)
ax.spines["bottom"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["left"].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.yaxis.grid()
colors = iter(jet(np.linspace(0,1,6)))
for i in amplitudes:
ds3['diffs_{}'.format(i)].plot(ax = ax, alpha = 0.9, color = next(colors), label="Difference {}".format(i), )
ax.legend()
plt.title("Spatial Averaged Temperature", fontsize=17, ha="center")
plt.text(-0.1, -0.2,"Data: Regional Climate Model (RCA3), Schimanke et. al, 2012",
horizontalalignment='left',
verticalalignment='center',
transform = ax.transAxes)
plt.tight_layout()
Store Sin fits
for i in amplitudes:
original_sin_fit = do_sin_fit(x, popt[0], i)
original_sin_fit = xr.DataArray(original_sin_fit, coords=[temp1d.time]).rename("data")
original_sin_fit.to_netcdf("/gfs1/work/mvkfbiow/postprocessing_baltic_past/data/"
+ "ref_for_temp_change_{}.nc".format(i))
Apply Change
Changes air temperature by applying sinus fit + variations
Note: ref_for_temp_change.nc
estimated on local machine in /work/publications/AMO_ERGOM/prep…
exp_dict = ({"original": "rca_temp_full.nc"})
for i in amplitudes:
exp_dict[str(i)] = "tairK_amplitude_{}.mom.data.nc".format(i)
plot = False
for count, year in enumerate(range(961, 1451)):
print("\r" + str(year) , end="")
ds = xr.open_mfdataset("/gfs2/work/mviowmod/DATABASE/BALTIC_MOM/forcing_950-1800/mom/{}/tairK.mom.dta.nc".format(year)).load()
ds2 = xr.open_dataset("/gfs1/work/mvkfbiow/postprocessing_baltic_past/data/ref_for_temp_change_{}.nc".format(i)).load()
ref_temp = 275.50696
ds3 = ds.copy(deep=True)
for index, month in enumerate(range(1, 13)):
current_selection = ds.sel(TIME='{}-{}'.format(str(year).zfill(4),str(month).zfill(2))).tairK.copy()
new = (current_selection - ref_temp) + ds2.sel(time=str(year+50).zfill(4)).isel(time=index).data.values
ds3.tairK.loc[dict(TIME='{}-{}'.format(str(year).zfill(4),str(month).zfill(2)))] = new.values
ds3.to_netcdf("/gfs1/work/mvkfbiow/forcing_950-1800/mom/{}/{}".format(year, exp_dict[str(i)]), unlimited_dims=['TIME'])
# if plot == True:
# diff = ds3 - ds#
#
# f, (ax,ax2,ax3) = plt.subplots(3)
# ds3.mean(["XLON", "YLAT"]).tairK.plot(ax=ax)
# ds.mean(["XLON", "YLAT"]).tairK.plot(ax=ax)
# ds3.isel(TIME=1).tairK.plot.pcolormesh(ax=ax2, vmin=260, vmax=290)
# ds.isel(TIME=1).tairK.plot.pcolormesh(ax=ax3, vmin=260, vmax=290)
Merge data for validation
print(exp_dict)
{'original': '/gfs1/work/mvkfbiow/postprocessing_baltic_past/data/rca_temp_full.nc', '0.4': 'tairK_amplitude_0.4.mom.data.nc', '0.5': 'tairK_amplitude_0.5.mom.data.nc', '0.6': 'tairK_amplitude_0.6.mom.data.nc', '0.7': 'tairK_amplitude_0.7.mom.data.nc', '0.8': 'tairK_amplitude_0.8.mom.data.nc', '0.9': 'tairK_amplitude_0.9.mom.data.nc'}
path_forcing = "/gfs1/work/mvkfbiow/forcing_950-1800/mom"
BEGIN = 951
END = 961
for count in amplitudes:
fname = exp_dict[str(count)]
for i in range(BEGIN, END):
cdo.monmean(input=path_forcing + "/" + str(i) + "/" + str(fname), output="/gfs1/work/mvkfbiow/postprocessing_baltic_past/data/" +
str(i).zfill(4) + "_" + fname)
cdo.mergetime(input="/gfs1/work/mvkfbiow/postprocessing_baltic_past/data/????_" + fname,
output="/gfs1/work/mvkfbiow/postprocessing_baltic_past/data/" + fname)
files = glob.glob("/gfs1/work/mvkfbiow/postprocessing_baltic_past/data/????_" + fname)
for f in files:
os.remove(f)
plt.figure(figsize=(12, 8))
# Remove the plot frame lines. They are unnecessary chartjunk.
ax = plt.subplot(111)
ax.spines["top"].set_visible(False)
ax.spines["bottom"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["left"].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.yaxis.grid()
for keys, name in exp_dict.items():
DS = xr.open_dataset("/gfs1/work/mvkfbiow/postprocessing_baltic_past/data/" + name)
try:
volume_flow = DS.tairK.mean(["XLON", "YLAT"])
except:
volume_flow = DS.temp.mean(["lon", "lat"])
volume_flow.name = name
volume_flow.plot(ax = ax, alpha = 0.3, color = "blue", label="", )
plt.title("tairK RCA3; shifted by 50 years", fontsize=17, ha="center")
plt.text(-0.1, -0.2,"Data: Regional Climate Model (RCA3), Schimanke et. al, 2012",
horizontalalignment='left',
verticalalignment='center',
transform = ax.transAxes)
plt.savefig("../tmp.png", dpi = 300)
Correct Humidity for new temperature
class specific_humidity():
"""
Calculates original RH and vapour pressure (e),
and saturation vapour pressure (es)
using temperature, pressure, humidity
Recalculates SH for new temperature using old RH and new
es to calculate e. e is used to calculate new SH.
es is calculated using Clausis-Clapeyron
"""
def __init__(self, p, q, t):
self.rd = 287.058
self.rv = 461.52
self.t = t
self.p = p
self.q = q
self.e = self.q * self.p / (0.378 * self.q + 0.622)
self.es = 611 * np.exp(17.67*(self.t-273.15)/(self.t-29.65))
self.rh = self.e / self.es
def calc_q(self, t):
es = 611 * np.exp(17.67*(t-273.15)/(t-29.65))
e = self.rh * es
self.q2 = ( self.rd * e / (self.rv*(self.p-e))).rename("shumi")
return self.q2
for i in amplitudes:
fname = exp_dict[str(i)]
print(fname)
for year in range(951, 961):
print("\r" + str(year) , end="")
t_new = xr.open_dataset("/gfs1/work/mvkfbiow/forcing_950-1800/mom/" + str(year) + "/" + fname)
p = xr.open_dataset("/gfs2/work/mviowmod/DATABASE/BALTIC_MOM/forcing_950-1800/mom/" + str(year) + "/pair.mom.dta.nc")
t = xr.open_dataset("/gfs2/work/mviowmod/DATABASE/BALTIC_MOM/forcing_950-1800/mom/" + str(year) + "/tairK.mom.dta.nc")
q = xr.open_dataset("/gfs2/work/mviowmod/DATABASE/BALTIC_MOM/forcing_950-1800/mom/" + str(year) + "/shumi.mom.dta.nc")
calc_hum = specific_humidity(p.pair, q.shumi, t.tairK)
s1 = calc_hum.calc_q(t_new.tairK).to_dataset(name = 'shumi')
s1.to_netcdf("/gfs1/work/mvkfbiow/forcing_950-1800/mom/" + str(year) + "/shumi" + fname.strip("tairK"), unlimited_dims=['TIME'])
tairK_amplitude_0.4.mom.data.nc
951
952
953
954
955
956
957
958
959
960
tairK_amplitude_0.5.mom.data.nc
951
952
953
954
955
956
957
958
959
960
tairK_amplitude_0.6.mom.data.nc
951
952
953
954
955
956
957
958
959
960
tairK_amplitude_0.7.mom.data.nc
951
952
953
954
955
956
957
958
959
960
tairK_amplitude_0.8.mom.data.nc
951
952
953
954
955
956
957
958
959
960
tairK_amplitude_0.9.mom.data.nc
951
952
953
954
955
956
957
958
959
960
Issues with grid?
had issues with saturation pressure overflow, cdo setgrid
resolved this
/gfs1/work/mvkfbiow/forcing_950-1800/mom/ contains .sh script with cdo setgrid
mygrid = "/gfs1/work/mvkfbiow/forcing_950-1800/mom/951/swdn.rco.dta.nc"
for i in amplitudes:
ifile = "shumi" + exp_dict[str(i)].strip("tairK")
print(ifile)
for year in range(952, 961):
cdo.setgrid(mygrid,
input="/gfs1/work/mvkfbiow/forcing_950-1800/mom/{}/".format(year) + ifile,
output="/gfs1/work/mvkfbiow/forcing_950-1800/mom/{}/".format(year) + "tmp.nc")
call("module load nco && ncatted -a calendar,TIME,o,c,julian /gfs1/work/mvkfbiow/forcing_950-1800/mom/{}/".format(year) + "tmp.nc", shell=True)
os.rename("/gfs1/work/mvkfbiow/forcing_950-1800/mom/{}/".format(year) + "tmp.nc",
"/gfs1/work/mvkfbiow/forcing_950-1800/mom/{}/".format(year) + ifile)
shumi_amplitude_0.4.mom.data.nc
shumi_amplitude_0.5.mom.data.nc
shumi_amplitude_0.6.mom.data.nc
shumi_amplitude_0.7.mom.data.nc
shumi_amplitude_0.8.mom.data.nc
shumi_amplitude_0.9.mom.data.nc
Validate new humidity
open_mfdataset( concat new dimension )
"""
enter year of interest
"""
year = 960
orig = xr.open_dataset("/gfs2/work/mviowmod/DATABASE/BALTIC_MOM/forcing_950-1800/mom/{}/shumi.mom.dta.nc".format(year))
orig_t = xr.open_dataset("/gfs2/work/mviowmod/DATABASE/BALTIC_MOM/forcing_950-1800/mom/{}/tairK.mom.dta.nc".format(year))
orig_gotland = orig.mean(["XLON", "YLAT"]).load()
orig_gotland_t = orig_t.mean(["XLON", "YLAT"]).load()
plt.figure(figsize=(12, 8))
jet= plt.get_cmap('Spectral_r')
colors = iter(jet(np.linspace(0,1,6)))
# Remove the plot frame lines. They are unnecessary chartjunk.
ax = plt.subplot(211)
ax.spines["top"].set_visible(False)
ax.spines["bottom"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["left"].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.yaxis.grid()
plt.title("Specific Humidity", fontsize=17, ha="center")
ax2 = plt.subplot(212)
ax2.spines["top"].set_visible(False)
ax2.spines["bottom"].set_visible(False)
ax2.spines["right"].set_visible(False)
ax2.spines["left"].set_visible(False)
ax2.get_xaxis().tick_bottom()
ax2.get_yaxis().tick_left()
ax2.yaxis.grid()
orig_gotland.shumi.plot(ax = ax, alpha = 0.3, color = "blue", label="Original", )
for i in amplitudes:
col = next(colors)
ifiles = "shumi" + exp_dict[str(i)].strip("tairK")
ifilet = exp_dict[str(i)]
new = xr.open_dataset("/gfs1/work/mvkfbiow/forcing_950-1800/mom/{}/{}".format(year, ifiles))
new_t = xr.open_dataset("/gfs1/work/mvkfbiow/forcing_950-1800/mom/{}/{}".format(year, ifilet))
new_gotland = new.mean(["lon", "lat"]).load()
new_gotland_t = new_t.mean(["XLON", "YLAT"]).load()
new_gotland.shumi.plot(ax = ax, alpha = 0.3, color = col, label="Shumi amp {}".format(i), )
(- orig_gotland_t.tairK + new_gotland_t.tairK).plot(ax = ax2, alpha = 0.9, color = col, label="Diff temp amp {}".format(i))
ax.legend()
ax2.legend()
plt.title("Temperature", fontsize=17, ha="center")
plt.text(-0.1, -0.2,"Data: Regional Climate Model (RCA3), Schimanke et. al, 2012",
horizontalalignment='left',
verticalalignment='center',
transform = ax.transAxes)
plt.tight_layout()
ds = xr.open_dataset("/gfs1/work/mvkfbiow/forcing_950-1800/mom/951/shumi.mom.dta.nc")
ds_s = xr.open_dataset("/gfs1/work/mvkfbiow/forcing_950-1800/mom/951/shumi_amplitude_05.mom.data.nc")
print(ds_s)
print(ds)
<xarray.Dataset>
Dimensions: (TIME: 3416, XLON: 129, YLAT: 89)
Coordinates:
* XLON (XLON) float64 7.5 7.75 8.0 8.25 8.5 ... 38.5 38.75 39.0 39.25 39.5
* YLAT (YLAT) float64 49.5 49.75 50.0 50.25 50.5 ... 70.75 71.0 71.25 71.5
* TIME (TIME) object 0950-12-01 00:00:00 ... 0952-01-31 21:00:00
Data variables:
shumi (TIME, YLAT, XLON) float32 ...
<xarray.Dataset>
Dimensions: (TIME: 3416, XLON: 129, YLAT: 89)
Coordinates:
* XLON (XLON) float64 7.5 7.75 8.0 8.25 8.5 ... 38.5 38.75 39.0 39.25 39.5
* YLAT (YLAT) float64 49.5 49.75 50.0 50.25 50.5 ... 70.75 71.0 71.25 71.5
* TIME (TIME) object 0950-12-01 00:00:00 ... 0952-01-31 21:00:00
Data variables:
shumi (TIME, YLAT, XLON) float32 ...
Attributes:
CDI: Climate Data Interface version 1.9.0 (http://mpimet.mpg.de/...
Conventions: CF-1.0
CDO: Climate Data Operators version 1.9.0 (http://mpimet.mpg.de/...
history: Thu Nov 23 11:56:52 2017: ncatted -a calendar,TIME,o,c,juli...
NCO: 4.6.7