import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters() # datetime converter for a matplotlib
import seaborn as sns
sns.set(style="ticks", font_scale=1.5)
from matplotlib.dates import DateFormatter
import matplotlib.dates as mdates
import matplotlib.ticker as ticker
import scipy as sp
import json
import requests
import os
import subprocess
from tqdm import tqdm
from scipy import signal
from scipy.signal import savgol_filter
# avoid "SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame."
pd.options.mode.chained_assignment = None # default='warn'71 savgol video
Import packages and stuff.
Download data from the IMS using an API.
# read token from file
with open('../archive/IMS-token.txt', 'r') as file:
TOKEN = file.readline()
# 28 = SHANI station
STATION_NUM = 28
start = "2022/01/01"
end = "2022/01/07"
filename = 'shani_2022_january.json'
# check if the JSON file already exists
# if so, then load file
if os.path.exists(filename):
with open(filename, 'r') as json_file:
data = json.load(json_file)
else:
# make the API request if the file doesn't exist
url = f"https://api.ims.gov.il/v1/envista/stations/{STATION_NUM}/data/?from={start}&to={end}"
headers = {'Authorization': f'ApiToken {TOKEN}'}
response = requests.get(url, headers=headers)
data = json.loads(response.text.encode('utf8'))
# save the JSON data to a file
with open(filename, 'w') as json_file:
json.dump(data, json_file)
# show data to see if it's alright
# dataLoad and process data.
df = pd.json_normalize(data['data'],record_path=['channels'], meta=['datetime'])
df['date'] = (pd.to_datetime(df['datetime'])
.dt.tz_localize(None) # ignores time zone information
)
df = df.pivot(index='date', columns='name', values='value')
df| name | Grad | RH | Rain | STDwd | TD | TDmax | TDmin | TG | TW | Time | WD | WDmax | WS | WS1mm | WSmax | Ws10mm |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| date | ||||||||||||||||
| 2022-01-01 00:00:00 | 0.0 | 77.0 | 0.0 | 10.3 | 11.2 | 11.2 | 11.1 | 10.7 | -9999.0 | 2354.0 | 75.0 | 64.0 | 5.0 | 6.0 | 7.0 | 5.5 |
| 2022-01-01 00:10:00 | 0.0 | 77.0 | 0.0 | 11.2 | 11.2 | 11.2 | 11.1 | 10.8 | -9999.0 | 1.0 | 77.0 | 84.0 | 4.7 | 5.5 | 6.6 | 4.9 |
| 2022-01-01 00:20:00 | 0.0 | 75.0 | 0.0 | 10.0 | 11.4 | 11.5 | 11.2 | 10.9 | -9999.0 | 20.0 | 80.0 | 83.0 | 5.1 | 6.2 | 8.0 | 5.1 |
| 2022-01-01 00:30:00 | 0.0 | 74.0 | 0.0 | 9.6 | 11.5 | 11.5 | 11.4 | 11.0 | -9999.0 | 22.0 | 76.0 | 74.0 | 4.8 | 5.8 | 7.3 | 5.0 |
| 2022-01-01 00:40:00 | 0.0 | 73.0 | 0.0 | 9.1 | 11.6 | 11.7 | 11.5 | 11.1 | -9999.0 | 34.0 | 74.0 | 64.0 | 4.8 | 5.7 | 7.2 | 5.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2022-01-06 23:10:00 | 0.0 | 36.0 | 0.0 | 16.1 | 11.6 | 12.0 | 11.1 | 6.8 | -9999.0 | 2310.0 | 144.0 | 126.0 | 0.7 | 1.6 | 2.0 | 0.7 |
| 2022-01-06 23:20:00 | 0.0 | 35.0 | 0.0 | 10.1 | 12.1 | 12.3 | 11.9 | 6.3 | -9999.0 | 2320.0 | 118.0 | 116.0 | 1.5 | 1.9 | 2.3 | 1.5 |
| 2022-01-06 23:30:00 | 0.0 | 36.0 | 0.0 | 7.1 | 12.4 | 12.6 | 11.9 | 7.3 | -9999.0 | 2330.0 | 113.0 | 116.0 | 2.5 | 3.0 | 3.3 | 2.5 |
| 2022-01-06 23:40:00 | 0.0 | 37.0 | 0.0 | 5.6 | 12.6 | 12.7 | 12.5 | 7.8 | -9999.0 | 2339.0 | 119.0 | 126.0 | 3.0 | 3.3 | 3.9 | 3.0 |
| 2022-01-06 23:50:00 | 0.0 | 39.0 | 0.0 | 11.5 | 11.9 | 12.6 | 11.5 | 7.1 | -9999.0 | 2341.0 | 102.0 | 108.0 | 1.9 | 2.5 | 2.8 | 2.9 |
864 rows × 16 columns
Define useful functions.
def concise(ax):
"""
Let python choose the best xtick labels for you
"""
locator = mdates.AutoDateLocator(minticks=3, maxticks=7)
formatter = mdates.ConciseDateFormatter(locator)
ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(formatter)
# dirty trick to have dates in the middle of the 24-hour period
# make minor ticks in the middle, put the labels there!
# from https://matplotlib.org/stable/gallery/ticks/centered_ticklabels.html
def center_dates(ax):
# show day of the month + month abbreviation. see full option list here:
# https://strftime.org
date_form = DateFormatter("%d %b")
# major ticks at midnight, every day
ax.xaxis.set_major_locator(mdates.DayLocator(interval=1))
ax.xaxis.set_major_formatter(date_form)
# minor ticks at noon, every day
ax.xaxis.set_minor_locator(mdates.HourLocator(byhour=[12]))
# erase major tick labels
ax.xaxis.set_major_formatter(ticker.NullFormatter())
# set minor tick labels as define above
ax.xaxis.set_minor_formatter(date_form)
# completely erase minor ticks, center tick labels
for tick in ax.xaxis.get_minor_ticks():
tick.tick1line.set_markersize(0)
tick.tick2line.set_markersize(0)
tick.label1.set_horizontalalignment('center')
def center_dates_two_panels(ax0, ax1):
# show day of the month + month abbreviation. see full option list here:
date_form = DateFormatter("%d %b")
# major ticks at midnight, every day
ax0.xaxis.set_major_locator(mdates.DayLocator(interval=1))
ax1.xaxis.set_major_locator(mdates.DayLocator(interval=1))
ax1.xaxis.set_major_formatter(date_form)
# minor ticks at noon, every day
ax1.xaxis.set_minor_locator(mdates.HourLocator(byhour=[12]))
# erase major tick labels
ax1.xaxis.set_major_formatter(ticker.NullFormatter())
# set minor tick labels as define above
ax1.xaxis.set_minor_formatter(date_form)
# completely erase minor ticks, center tick labels
for tick in ax0.xaxis.get_minor_ticks():
tick.tick1line.set_markersize(0)
tick.tick2line.set_markersize(0)
for tick in ax1.xaxis.get_minor_ticks():
tick.tick1line.set_markersize(0)
tick.tick2line.set_markersize(0)
tick.label1.set_horizontalalignment('center')We don’t need the full month, let’s cut the dataframe to fewer days.
We now redefine a narrower window, this will be the graph’s xlimits. We leave the dataframe as is, because we will need some data outside the graph’s limits.
fig, ax = plt.subplots(figsize=(8,5))
ax.plot(df.loc[start:end, 'TD'], color='black')
ax.set(ylim=[5, 17.5],
xlim=[start, end],
ylabel="Temperature (°C)",
title="Yatir Forest, 2022",
yticks=[5,10,15])
center_dates(ax)
# fig.savefig("sliding_YF_temperature_2022.png")Looks good. Let’s move on.
71.1 Savgol filter
%matplotlib widget
fig, ax = plt.subplots(figsize=(8,5))
ax.plot(df.loc[start:end, 'TD'], color='black')
sg = savgol_filter(df['TD'], 13, 2)
i = 500
ax.plot(df.index[:i], sg[:i], color='xkcd:hot pink')
window_pts = 31
p_order = 3
window_x = np.arange(i - window_pts // 2, i + window_pts // 2)
window_y = df['TD'][i - window_pts // 2:i + window_pts // 2]
# Fit and plot polynomial inside the window
fitted_y, coeffs = fit_polynomial(window_x, window_y, p_order)
whole_x = np.arange(len(df))
whole_y = df['TD'].values
poly = np.poly1d(coeffs)
whole_poly = poly(whole_x)
ax.plot(df.index, whole_poly, color='xkcd:sun yellow', lw=2)
# ax.plot(df.index[window_x], fitted_y, color='0.8', lw=3)
ax.plot(df.index[window_x], fitted_y, color='xkcd:mustard', lw=2)
ax.set(ylim=[5, 17.5],
xlim=[start, end],
ylabel="Temperature (°C)",
title="Yatir Forest, 2022",
yticks=[5,10,15])
center_dates(ax)
# fig.savefig("sliding_YF_temperature_2022.png")
%matplotlib widget
fig, ax = plt.subplots(figsize=(8,5))
class Lines:
"""
empty class, later will be populated with graph objects.
this is useful to draw and erase lines on demand.
"""
pass
lines = Lines()
# set graph y limits
ylim = [3, 22]
# choose here windown width in minutes
window_width_min = 500.0
window_width_min_integer = int(window_width_min) # same but integer
window_width_int = int(window_width_min // 10 + 1) # window width in points
N = len(df) # df length
t_swipe = pd.date_range(start=pd.to_datetime(start) - pd.Timedelta(minutes=window_width_min) - pd.Timedelta(minutes=30),
end=pd.to_datetime(end) + pd.Timedelta(minutes=60),
freq="10min")
# starting time
t0 = t_swipe[0]
ind0 = df.index.get_loc(t0) + window_width_int//2 + 1
# show sliding window on the top panel as a light blue shade
lines.fill_bet = ax.fill_between([t0, t0 + pd.Timedelta(minutes=window_width_min)],
y1=ylim[0], y2=ylim[1], alpha=0.1, zorder=-1)
sg = savgol_filter(df['TD'], window_width_int, p_order)
df.loc[:, 'sg'] = sg
# plot temperature
ax.plot(df.loc[start:end, 'TD'], color="black")
# define x,y data inside window to execute polyfit on
window_x = np.arange(ind0 - window_width_int // 2, ind0 + window_width_int // 2)
window_y = df['TD'][ind0 - window_width_int // 2:ind0 + window_width_int // 2].values
# fit and plot polynomial inside the window
fitted_y, coeffs = fit_polynomial(window_x, window_y, p_order)
# get x,y data for the whole array
whole_x = np.arange(len(df))
whole_y = df['TD'].values
poly = np.poly1d(coeffs)
whole_poly = poly(whole_x)
# calculate the middle of the sliding window
window_middle = t0 + pd.Timedelta(minutes=window_width_min/2)
# plot a pink line showing the result of the moving average
# from the beginning to the middle of the sliding window
lines.pink_line, = ax.plot(df.loc[start:window_middle, 'sg'], color="xkcd:hot pink", lw=3)
lines.poly_all, = ax.plot(df.index, whole_poly, color='xkcd:sun yellow', lw=2)
lines.poly_window, = ax.plot(df.index[window_x], fitted_y, color='xkcd:mustard', lw=2)
# emphasize the location of the middle on the window with a circle
lines.pink_circle, = ax.plot([window_middle], [df.loc[window_middle, 'sg']],
marker='o', markerfacecolor="None", markeredgecolor="xkcd:dark pink", markeredgewidth=2,
markersize=8)
# some explanation
ax.text(0.99, 0.97, f"savitzky-golay\nwidth = {window_width_int:.0f} pts\npoly order = {p_order}", transform=ax.transAxes,
horizontalalignment='right', verticalalignment='top',
fontsize=14)
# axis tweaking
ax.set(ylim=ylim,
xlim=[start, end],
ylabel="Temperature (°C)",
yticks=[5,10,15,20],
title="Yatir Forest, 2022")
# adjust dates on both panels as defined before
center_dates(ax)
def update_swipe(k, lines):
"""
updates both panels, given the index k along which the window is sliding
"""
# left side of the sliding window
t0 = t_swipe[k]
# middle position
window_middle = t0 + pd.Timedelta(minutes=window_width_min/2)
ind0 = df.index.get_loc(window_middle)
# erase previous blue shade on the top graph
lines.fill_bet.remove()
# fill again the blue shade in the updated window position
lines.fill_bet = ax.fill_between([t0, t0 + pd.Timedelta(minutes=window_width_min)],
y1=ylim[0], y2=ylim[1], alpha=0.1, zorder=-1, color="tab:blue")
# update pink curve
lines.pink_line.set_data(df[start:window_middle].index,
df.loc[start:window_middle, 'sg'].values)
# update pink circle
lines.pink_circle.set_data([window_middle], [df.loc[window_middle, 'sg']])
# define x,y data inside window to execute polyfit on
window_x = np.arange(ind0 - window_width_int // 2, ind0 + window_width_int // 2)
window_y = df['TD'][ind0 - window_width_int // 2:ind0 + window_width_int // 2]
# fit and plot polynomial inside the window
fitted_y, coeffs = fit_polynomial(window_x, window_y, p_order)
poly = np.poly1d(coeffs)
whole_poly = poly(whole_x)
lines.poly_all.set_data(df.index, whole_poly)
lines.poly_window.set_data(df.index[window_x], fitted_y)
fig.savefig(f"pngs/savgol{window_width_int}/savgol_zero.png", dpi=600)
# create a tqdm progress bar
progress_bar = tqdm(total=len(t_swipe), unit="iteration")
# loop over all sliding indices, update graph and then save it
for fignum, i in enumerate(np.arange(0, len(t_swipe)-1, 1)):
update_swipe(i, lines)
fig.savefig(f"pngs/savgol{window_width_int}/savgol_{window_width_int}_{fignum:03}.png", dpi=600)
# update the progress bar
progress_bar.update(1)
# close the progress bar
progress_bar.close()100%|█████████▉| 634/635 [13:07<00:01, 1.24s/iteration]Combine all saved images into one mp4 video.
# Define the path to your PNG images
pngs_path = f"pngs/savgol51"
pngs_name = f"savgol_51_%03d.png"
# Define the output video file path
video_output = f"output_savgol51.mp4"
# Use ffmpeg to create a video from PNG images
# desired framerate. choose 24 if you don't know what to do
fr = 12
# run command
ffmpeg_cmd = f"ffmpeg -framerate {fr} -i {pngs_path}/{pngs_name} -c:v libx264 -vf fps={fr} {video_output}"
subprocess.run(ffmpeg_cmd, shell=True)ffmpeg version 6.1.1 Copyright (c) 2000-2023 the FFmpeg developers
built with Apple clang version 15.0.0 (clang-1500.1.0.2.5)
configuration: --prefix=/usr/local/Cellar/ffmpeg/6.1.1_2 --enable-shared --enable-pthreads --enable-version3 --cc=clang --host-cflags= --host-ldflags='-Wl,-ld_classic' --enable-ffplay --enable-gnutls --enable-gpl --enable-libaom --enable-libaribb24 --enable-libbluray --enable-libdav1d --enable-libharfbuzz --enable-libjxl --enable-libmp3lame --enable-libopus --enable-librav1e --enable-librist --enable-librubberband --enable-libsnappy --enable-libsrt --enable-libsvtav1 --enable-libtesseract --enable-libtheora --enable-libvidstab --enable-libvmaf --enable-libvorbis --enable-libvpx --enable-libwebp --enable-libx264 --enable-libx265 --enable-libxml2 --enable-libxvid --enable-lzma --enable-libfontconfig --enable-libfreetype --enable-frei0r --enable-libass --enable-libopencore-amrnb --enable-libopencore-amrwb --enable-libopenjpeg --enable-libspeex --enable-libsoxr --enable-libzmq --enable-libzimg --disable-libjack --disable-indev=jack --enable-videotoolbox --enable-audiotoolbox
libavutil 58. 29.100 / 58. 29.100
libavcodec 60. 31.102 / 60. 31.102
libavformat 60. 16.100 / 60. 16.100
libavdevice 60. 3.100 / 60. 3.100
libavfilter 9. 12.100 / 9. 12.100
libswscale 7. 5.100 / 7. 5.100
libswresample 4. 12.100 / 4. 12.100
libpostproc 57. 3.100 / 57. 3.100
Input #0, image2, from 'pngs/savgol51/savgol_51_%03d.png':
Duration: 00:00:52.83, start: 0.000000, bitrate: N/A
Stream #0:0: Video: png, rgba(pc, gbr/unknown/unknown), 4800x3000 [SAR 23622:23622 DAR 8:5], 12 fps, 12 tbr, 12 tbn
Stream mapping:
Stream #0:0 -> #0:0 (png (native) -> h264 (libx264))
Press [q] to stop, [?] for help
[libx264 @ 0x7f9cb7906dc0] using SAR=1/1
[libx264 @ 0x7f9cb7906dc0] using cpu capabilities: MMX2 SSE2Fast SSSE3 SSE4.2 AVX FMA3 BMI2 AVX2
[libx264 @ 0x7f9cb7906dc0] profile High 4:4:4 Predictive, level 6.0, 4:4:4, 8-bit
[libx264 @ 0x7f9cb7906dc0] 264 - core 164 r3108 31e19f9 - H.264/MPEG-4 AVC codec - Copyleft 2003-2023 - http://www.videolan.org/x264.html - options: cabac=1 ref=3 deblock=1:0:0 analyse=0x3:0x113 me=hex subme=7 psy=1 psy_rd=1.00:0.00 mixed_ref=1 me_range=16 chroma_me=1 trellis=1 8x8dct=1 cqm=0 deadzone=21,11 fast_pskip=1 chroma_qp_offset=4 threads=6 lookahead_threads=1 sliced_threads=0 nr=0 decimate=1 interlaced=0 bluray_compat=0 constrained_intra=0 bframes=3 b_pyramid=2 b_adapt=1 b_bias=0 direct=1 weightb=1 open_gop=0 weightp=2 keyint=250 keyint_min=12 scenecut=40 intra_refresh=0 rc_lookahead=40 rc=crf mbtree=1 crf=23.0 qcomp=0.60 qpmin=0 qpmax=69 qpstep=4 ip_ratio=1.40 aq=1:1.00
Output #0, mp4, to 'output_savgol51.mp4':
Metadata:
encoder : Lavf60.16.100
Stream #0:0: Video: h264 (avc1 / 0x31637661), yuv444p(tv, progressive), 4800x3000 [SAR 1:1 DAR 8:5], q=2-31, 12 fps, 12288 tbn
Metadata:
encoder : Lavc60.31.102 libx264
Side data:
cpb: bitrate max/min/avg: 0/0/0 buffer size: 0 vbv_delay: N/A
[out#0/mp4 @ 0x7f9cb7806000] video:3513kB audio:0kB subtitle:0kB other streams:0kB global headers:0kB muxing overhead: 0.225481%
frame= 634 fps=3.2 q=-1.0 Lsize= 3521kB time=00:00:52.58 bitrate= 548.5kbits/s speed=0.27x
[libx264 @ 0x7f9cb7906dc0] frame I:3 Avg QP:13.70 size:146882
[libx264 @ 0x7f9cb7906dc0] frame P:232 Avg QP:19.02 size: 7856
[libx264 @ 0x7f9cb7906dc0] frame B:399 Avg QP:24.04 size: 3341
[libx264 @ 0x7f9cb7906dc0] consecutive B-frames: 11.4% 10.7% 10.4% 67.5%
[libx264 @ 0x7f9cb7906dc0] mb I I16..4: 32.8% 61.0% 6.2%
[libx264 @ 0x7f9cb7906dc0] mb P I16..4: 0.5% 0.7% 0.3% P16..4: 0.4% 0.2% 0.1% 0.0% 0.0% skip:97.7%
[libx264 @ 0x7f9cb7906dc0] mb B I16..4: 0.1% 0.0% 0.0% B16..8: 1.9% 0.3% 0.0% direct: 0.0% skip:97.6% L0:50.4% L1:47.8% BI: 1.7%
[libx264 @ 0x7f9cb7906dc0] 8x8 transform intra:51.3% inter:35.8%
[libx264 @ 0x7f9cb7906dc0] coded y,u,v intra: 7.2% 6.5% 4.4% inter: 0.1% 0.1% 0.0%
[libx264 @ 0x7f9cb7906dc0] i16 v,h,dc,p: 89% 10% 1% 0%
[libx264 @ 0x7f9cb7906dc0] i8 v,h,dc,ddl,ddr,vr,hd,vl,hu: 31% 3% 65% 0% 0% 0% 0% 0% 0%
[libx264 @ 0x7f9cb7906dc0] i4 v,h,dc,ddl,ddr,vr,hd,vl,hu: 54% 7% 23% 3% 1% 4% 1% 5% 1%
[libx264 @ 0x7f9cb7906dc0] Weighted P-Frames: Y:0.0% UV:0.0%
[libx264 @ 0x7f9cb7906dc0] ref P L0: 58.7% 5.7% 25.2% 10.4%
[libx264 @ 0x7f9cb7906dc0] ref B L0: 85.9% 11.8% 2.4%
[libx264 @ 0x7f9cb7906dc0] ref B L1: 96.9% 3.1%
[libx264 @ 0x7f9cb7906dc0] kb/s:544.58CompletedProcess(args='ffmpeg -framerate 12 -i pngs/savgol51/savgol_51_%03d.png -c:v libx264 -vf fps=12 output_savgol51.mp4', returncode=0)