In [163]:
%matplotlib inline
%config InlineBackend.figure_format='retina'
from datacharm import *
from enerpi.api import enerpi_data_catalog
from enerpi.enerplot import plot_tile_last_24h, plot_power_consumption_hourly
import random
import math
import pysolar
import datetime as dt
import pytz
import matplotlib.dates as mpd
LAT, LONG = 38.631463, -0.866402
TZ = pytz.timezone('Europe/Madrid')
def _tslim(ax, h0=12, hf=22, m0=0, mf=0):
xlim = [mpd.num2date(x, tz=TZ) for x in ax.get_xlim()]
new = [mpd.date2num(d.replace(hour=h, minute=m)) for d, h, m in zip(xlim, [h0, hf], [m0, mf])]
ax.set_xlim(new)
return ax
# Catálogo y lectura de todos los datos.
cat = enerpi_data_catalog()
data, data_s = cat.get_all_data()
print_info(data_s.tail())
LDR = pd.DataFrame(data.ldr).tz_localize(TZ)
print_cyan(LDR.describe())
LDR.hist(bins=(LDR.max() - LDR.min() + 1).values[0] // 5, log=True, figsize=(18, 8));
In [ ]:
In [ ]:
In [81]:
# Plot day
def _get_rday(data):
t0 = data.index[0].date()
tf = data.index[-1].date()
delta = tf - t0
rd = random.randrange(0, delta.days + 1, 1)
day = t0 + dt.timedelta(days=rd)
return data.loc[day.strftime('%Y-%m-%d'):day.strftime('%Y-%m-%d')]
def _alt_azi(d):
return pysolar.solar.get_altitude(LAT, LONG, d), pysolar.solar.get_azimuth(LONG, LAT, d)
def get_solar_day(day=dt.datetime.today(), step_minutes=5, tz=TZ, step_calc_seg=30):
if (type(day) is pd.DataFrame) or (type(day) is pd.Series):
df_sol = pd.DataFrame(day)
sunrise, sunset = pysolar.util.get_sunrise_sunset(LAT, LONG, df_sol.index[0])
df_sol.loc[(sunrise > df_sol.index) | (sunset < df_sol.index), 'altitud'] = 0
idx_calc_solar = df_sol.index[df_sol['altitud'].isnull()][::step_calc_seg]
alts = [pysolar.solar.get_altitude(LAT, LONG, d) for d in idx_calc_solar]
azims = [pysolar.solar.get_azimuth(LAT, LONG, d) for d in idx_calc_solar]
irrads = [pysolar.radiation.get_radiation_direct(d, alt) if alt > 0 else 0 for d, alt in zip(idx_calc_solar, alts)]
df_sol.loc[idx_calc_solar, 'altitud'] = alts
df_sol.loc[idx_calc_solar, 'azimut'] = azims
df_sol.loc[idx_calc_solar, 'irradiacion_cs'] = irrads
#df_sol.altitud = df_sol.altitud.where(df_sol.altitud > 0, other=0)
#df_sol['irradiacion_cs'] = df_sol.apply(lambda x: pysolar.radiation.get_radiation_direct(x.name, x['altitud']), axis=1)
df_sol = df_sol.interpolate().fillna(0)
else:
day = tz.localize(day).replace(hour=0, minute=0, second=0, microsecond=0)
tt = [day + dt.timedelta(minutes=m) for m in range(0, 24 * 60, step_minutes)]
df_sol = pd.DataFrame([_alt_azi(d) for d in tt], columns=['altitud', 'azimut'], index=tt)
df_sol.altitud = df_sol.altitud.where(df_sol.altitud > 0, other=0)
df_sol['irradiacion_cs'] = df_sol.apply(lambda x: pysolar.radiation.get_radiation_direct(x.name, x['altitud']), axis=1)
return df_sol
df_sol = get_solar_day(day=dt.datetime.today(), step_minutes=5)
data_day = get_solar_day(_get_rday(LDR))
data_day.head()
Out[81]:
In [82]:
df_sol_p = data_day.copy()
print_red(df_sol_p.describe())
df_sol_p.altitud *= 10.
#df_sol_p.azimut += 180.
#df_sol_p.azimut *= -1.
df_sol_p[['ldr', 'altitud']].plot(figsize=(18, 6));
print_magenta('AZIMUT RANGE: 0 --> -360 sentido horario, comienza en SUR. W=-90, N=-180, E=-270, S=0/-360')
In [298]:
idx_calc_solar[0]
#[pysolar.solar.get_azimuth(LAT, LONG, d) for d in idx_calc_solar]
Out[298]:
In [126]:
d_rs['d_ldr'] = d_rs.ldr.diff().fillna(0)
d_rs['on'] = d_rs['d_ldr'] > 100
d_rs['off'] = d_rs['d_ldr'] < -100
ax = d_rs[['ldr', 'altitud']].plot(figsize=(18, 6));
d_rs[d_rs['on']].ldr.plot(lw=0, marker='o', markersize=10, color='red')
d_rs[d_rs['off']].ldr.plot(lw=0, marker='o', markersize=10, color='blue')
_tslim(ax, 20, 21)
Out[126]:
In [315]:
def _azimut_S(azi):
if azi + 180 < 0:
return -(360 + azi)
return -azi
azimuts = pd.Series([_azimut_S(pysolar.solar.get_azimuth(LAT, LONG, d)) for d in idx_calc_solar], index=idx_calc_solar)
azimuts.plot()
Out[315]:
In [ ]:
In [ ]:
In [ ]:
In [266]:
def _f_roll(x):
delta_pos = np.sum(x[x > 0])
delta_neg = np.sum(x[x < 0])
return delta_pos + delta_neg
#def _f_state(delta):
# if delta > 100:
# return 1
# elif delta < -100:
# return -1
# return 0
def _min_max_rel(x):
mid = x[len(x) // 2]
if mid == np.max(x):
return 1
elif mid == np.min(x):
return -1
return 0
def _first_on(x):
if (np.sum(x) == 1.) & (x[-1] == 1.):
return True
return False
def _last_on(x):
if (np.sum(x) == 1.) & (x[-1] == 1.):
return True
return False
# TODO Resolver NaN's
resample_inicial = '2s'
delta_roll_threshold = 100
data_analysis = data_day[['ldr', 'altitud', 'azimut', 'irradiacion_cs']].resample(resample_inicial).mean()
data_analysis['delta_roll'] = data_analysis.ldr.diff().fillna(0).rolling(3, center=True).apply(_f_roll).fillna(0)
data_analysis['on'] = data_analysis['delta_roll'] > delta_roll_threshold
data_analysis['off'] = data_analysis['delta_roll'] < -delta_roll_threshold
data_analysis['on'].rolling?
data_analysis['ch_state'] = 0
data_analysis.loc[data_analysis[data_analysis['on'].rolling(3).apply(_first_on).fillna(False) > 0].index, 'ch_state'] = 1
data_analysis.loc[data_analysis[data_analysis['off'].rolling(3).apply(_last_on).fillna(False) > 0].index, 'ch_state'] = -1
#print_red(data_analysis[data_analysis['on']])
#data_analysis['light_state'] = data_analysis['delta_roll'].apply(_f_state)
#data_analysis['state_ch'] = data_analysis['delta_roll'].rolling(3, center=True).apply(_min_max_rel)
#encendidos = data_analysis[data_analysis['on'].rolling(3).apply(_first_on).fillna(False) > 0]
#data_analysis[data_analysis['on'] & (data_analysis['state_ch'] == 1)]
#apagados = data_analysis[data_analysis['off'].rolling(3).apply(_last_on).fillna(False) > 0]
#data_analysis[data_analysis['off'] & (data_analysis['state_ch'] == -1)]
#print_info(encendidos[['ldr','delta_roll', 'on', 'off', 'state_ch']])
#apagados[['ldr','delta_roll', 'on', 'off', 'state_ch']]
#[['ldr','delta_roll', 'on', 'off', 'state_ch']]
cambios = data_analysis[data_analysis['ch_state'] != 0]
cambios
Out[266]:
In [ ]:
In [ ]:
In [256]:
In [ ]:
In [271]:
encendidos = cambios[cambios.ch_state == 1].index
apagados = cambios[cambios.ch_state == -1].index
data_analysis['artif_light'] = False
data_analysis['artif_level_max'] = 0
data_analysis['artif_level'] = 0
for i_enc in encendidos:
apagado = apagados[apagados > i_enc]
if len(apagado) > 0:
apagado = apagado[0]
if apagado:
data_analysis.loc[i_enc:apagado, 'artif_light'] = True
data_analysis.loc[i_enc:apagado, 'artif_level_max'] = data_analysis.ldr.loc[i_enc:apagado].max()
#df_on = data_analysis.loc[i_enc: apagado]
print_red('Encendido L={:.0f}, SEGS={:.0f}'.format(data_analysis.ldr.loc[i_enc:apagado].max(), (apagado - i_enc).total_seconds()))
#print_cyan(df_on.ldr.describe().T)
#ax = data_analysis.loc[i_enc - pd.Timedelta('1min'): apagado + pd.Timedelta('1min')].ldr.plot(figsize=(18, 5))
#ax.scatter(x=i_enc, y=data_analysis.loc[i_enc, 'ldr'], color='red', marker='o', lw=0)
#ax.scatter(x=apagado, y=data_analysis.loc[apagado, 'ldr'], color='blue', marker='o', lw=0)
#plt.show()
ax = data_analysis['artif_level_max'].plot(figsize=(18, 10))
_tslim(ax, 0, 2)
plt.show()
ax = data_analysis['artif_level_max'].plot(figsize=(18, 10))
_tslim(ax, 20, 21)
Out[271]:
In [273]:
#data_analysis[apagado, 'ldr']
data_analysis[['ldr', 'artif_level_max', 'altitud']].plot()
Out[273]:
In [283]:
# Reconstrucción LDR natural:
#data_analysis[data_analysis.artif_level_max > 0].index
index_art = data_analysis[(data_analysis.artif_level_max > 0).shift(-2) | (data_analysis.artif_level_max > 0).shift(2)].index
data_analysis.ldr.drop(index_art).resample('1min').max().interpolate().plot(figsize=(18, 8))
Out[283]:
In [291]:
ax = data_analysis.ldr.drop(index_art).resample('2min').max().interpolate().plot(figsize=(18, 8), color='darkorange')
data_analysis.artif_level_max.plot(ax=ax, color='red')
_tslim(ax, 0, 2)
plt.show()
ax = data_analysis.ldr.drop(index_art).resample('2min').max().interpolate().plot(figsize=(18, 8), color='darkorange')
data_analysis.artif_level_max.plot(ax=ax, color='red')
_tslim(ax, 20, 22)
ax = data_analysis.ldr.drop(index_art).resample('2min').max().interpolate().plot(figsize=(18, 8), color='darkorange')
data_analysis.artif_level_max.plot(ax=ax, color='red')
Out[291]:
In [239]:
ax = data_day.iloc[::5].ldr.plot(figsize=(18, 8), color='darkorange')
encendidos.ldr.plot(marker='d', color='red', lw=0, ax=ax)
apagados.ldr.plot(marker='o', color='blue', lw=0, ax=ax)
Out[239]:
In [296]:
data_analysis.join(data_analysis.ldr.drop(index_art).resample('2min').max().interpolate().rename('natural')).interpolate()
Out[296]:
In [293]:
data_analysis[['altitud', 'azimut', ]]
Out[293]:
In [226]:
g_horas_enc = encendidos.groupby(lambda x: x.hour).groups
for hora, idx_enc in g_horas_enc.items():
ax = data_day.ldr.plot(figsize=(18, 6), color='darkorange')
encendidos.loc[idx_enc].ldr.plot(marker='d', color='red', lw=0, ax=ax)
_tslim(ax, hora, hora + 1)
plt.show()
break
In [225]:
def _min_max_rel_m(x):
print(type(x), len(x), x)
#assert()
mid = x[len(x) // 2]
if mid == np.max(x):
return 1
elif mid == np.min(x):
return -1
return 0
data_analysis.iloc[:20][['ldr', 'delta_roll']].rolling(3, center=True).apply(_min_max_rel_m)
Out[225]:
In [152]:
dday = _get_rday(LDR)
#.resample('5min').fillna(method='bfill').round(0).astype(int)
print_info(dday.head())
print_cyan(len(dday))
dday['delta_v'] = (dday['ldr'].shift(1) - dday['ldr']).fillna(0)
only_ch = dday[dday.delta_v != 0]
print_cyan(len(only_ch))
only_ch.plot()
Out[152]:
In [157]:
# Limpieza de LDR:
# * ROLLING
dday['event'] = dday.delta_v.abs() > 50
dday.loc[dday['event'], 'event_s'] = dday[dday['event']].delta_v > 0
dday.sum()
Out[157]:
In [159]:
dday[dday.event & dday.event_s]
Out[159]:
In [174]:
rs = dday.resample('1min')
ldr_dayrs = rs.ldr.apply(lambda x: pd.Series({'v_min': np.min(x),
'v_max': np.max(x),
'v_median': np.median(x)}) #.dropna().astype(int)
).unstack()
#ldr_dayrs = ldr_dayrs.join((rs.event.sum() > 0) * 500)
ax = ldr_dayrs.plot(figsize=(18, 8), alpha=.8, lw=1)
(rs.event.sum() > 0) * 700
Out[174]:
In [181]:
import matplotlib.dates as mpd
def _xlim_time(ax, t0, tf):
xlim = ax.get_xlim()
mpd.num.xlim[0]
dday.ldr.rolling(5).max().plot()
plt.gca().get_xlim()
Out[181]:
In [ ]: