yes, definitely looks like there is no data.
I have attached another version of mine, in which the trend line is
disabled by default, but I suspect that would just delay the inevitable and
it would crash trying to do the plot.
I also fixed up a few plotting errors in my code to do with the mysteries
(to me) of layer ordering.
I also had a background bar showing either side of expected arrival - in
this version I have now changed that to start at the expected arrival and
stop 1 hour later.
On Friday, 21 January 2022 at 3:29:36 am UTC+10 [email protected] wrote:
> Hello,
> Not being a programmer, I probably shouldn't have messed with this, but
> being curious...
>
> I tried the code posted on github as well as the one by Cameron D. In both
> cases I got the following error:
> ```
> root@n4mrv:/home/bg/weewx_tonga_browse-main# python3 ./tonga.py [file
> from Cameron D]
>
> distance to eruption 12056.6 km arrival at 1642258360 (2022-01-15 09:52:39)
> opposite pulse arrival at 1642308921 (2022-01-15 23:55:21)
> second time around pulse arrival at 1642385471 (2022-01-16 21:11:10)
> Traceback (most recent call last):
> File "./tonga.py", line 178, in <module>
> plot_burst( cursor, arrival_time, hour_span, "primary" )
> File "./tonga.py", line 54, in plot_burst
> coeff = np.polyfit(xdata, ydata, background_order )
> File "<__array_function__ internals>", line 180, in polyfit
> File "/usr/local/lib/python3.8/dist-packages/numpy/lib/polynomial.py",
> line 638, in polyfit
> raise TypeError("expected non-empty vector for x")
> TypeError: expected non-empty vector for x
> ```
> I added my lat/lon information but may have missed something else I need
> to change. Python modules were installed as directed. Copy of weewx.sdb is
> in the same directory as the program.
> Thanks.
> Bob
> On Wednesday, January 19, 2022 at 10:32:49 AM UTC-5 [email protected]
> wrote:
>
>> On Wednesday, January 19, 2022 at 12:42:04 a.m. UTC-4 Cameron D wrote:
>>
>>>
>>> - as you get closer to the equator, tidal changes dominate the
>>> baseline in that timescale - I tried higher order polynomials, but they
>>> are
>>> next to useless.
>>>
>>> I also had little luck with higher order polynomials to remove the
>> general trend.
>>
>> I've put the script here: https://github.com/morrowwm/weewx_tonga_browse
>>
>
--
You received this message because you are subscribed to the Google Groups
"weewx-user" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/weewx-user/71655988-a5d6-4771-bb56-12fa455f24f2n%40googlegroups.com.
import time
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from geopy import distance
#DB="sqlite3"
DB="mysql"
if DB == "sqlite3":
import sqlite3
else:
import sys
import mysql.connector as MysqlConnector
from mysql.connector import errorcode as msqlerrcode
db_user="weewx"
db_pwd="insert-your-password"
db_name="weewx"
# change as needed - (latitude, longitude)
home = (44.8032, -63.6204)
# some barometers only change their reading every nth record in the database
# set oversampling = 1 if the barometer can change every archive record.
# set oversampling = 2 if the barometer updates every 10 minutes but you have a 5 minute sample period
oversampling = 1
# this will remove background if you have a smoothly varying long-term barometer trend
plot_trends=False
# order of polynomial fitted to background
background_order = 2
# you should not need to edit anything below here
# ===============================================
def plot_burst( cursor, pulse_time, span, order ):
start_time = 3600*(pulse_time // 3600 - span) # start at an even hour
xmax_hours = pulse_time + span*3600
query = "select datetime, barometer from archive where datetime > %.0f and datetime < %.0f order by dateTime;" %\
( start_time, xmax_hours )
cursor.execute(query)
result = cursor.fetchall()
xdata = []
ydata = []
tdata = []
rowcount = 0
for row in result:
rowcount += 1
if rowcount % oversampling == 0:
xdata.append(row[0])
tdata.append(mdates.datestr2num(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime(row[0]))))
ydata.append(row[1])
peakt = [mdates.datestr2num(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime( pulse_time+ 1800 )))]
plotx_range = (
mdates.datestr2num(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime( start_time ))),
mdates.datestr2num(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime( xmax_hours )))
)
bar_width = 1.0/24 # one hour
if plot_trends:
# remove any linear trend
coeff = np.polyfit(xdata, ydata, background_order )
trend = np.polyval(coeff, xdata)
fig, ax = plt.subplots(figsize=[12,10])
plt.text(tdata[1], np.max(ydata), "location: %0.2f, %0.2f" % home)
date_formatter = mdates.DateFormatter('%m-%d %H:%M')
ax.tick_params(axis="x", rotation=90)
ax.xaxis.set_major_locator(mdates.HourLocator(interval = 1))
ax.xaxis.set_major_formatter(date_formatter)
ax.xaxis.set_minor_locator(mdates.MinuteLocator(interval = 10))
fig.subplots_adjust(bottom=0.2)
ax3 = ax.twinx()
ax3.set_axis_off()
ax3.set_zorder(1)
ax3.bar( peakt, [1] , width=bar_width, color="lightgrey")
ax.set_xlim(left=plotx_range ) # start at a round hour
ax.set_zorder(2)
ax.patch.set_visible(False)
ax.plot(tdata, ydata, color="blue", marker="+", linewidth=1)
if plot_trends:
ax.plot(tdata, trend, color="lightgrey", linewidth=1)
ax2=ax.twinx()
ax2.plot(tdata, ydata - trend, color="lightblue", linewidth=3)
# save the plot as a file
filename = "./Hunga_Tonga_" + order + ".png"
fig.savefig( filename, format='png', dpi=100, bbox_inches='tight')
#plt.show()
def plot_all( cursor, e_time, pulse_time1, pulse_time2, pulse_time3, span ):
start_time = 3600*(e_time // 3600 - span) # start at an even hour
xmax_hours = pulse_time3 + span*3600
query = "select datetime, barometer from archive where datetime > %.0f and datetime < %.0f order by dateTime;" %\
( start_time, xmax_hours )
cursor.execute(query)
result = cursor.fetchall()
xdata = []
ydata = []
tdata = []
rowcount = 0
for row in result:
rowcount += 1
if rowcount % oversampling == 0:
xdata.append(row[0])
tdata.append(mdates.datestr2num(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime(row[0]))))
ydata.append(row[1])
peakt = [
mdates.datestr2num(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime( pulse_time1 + 1800 ))),
mdates.datestr2num(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime( pulse_time2 + 1800 ))),
mdates.datestr2num(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime( pulse_time3 + 1800 )))
]
plotx_range = (
mdates.datestr2num(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime( start_time ))),
mdates.datestr2num(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime( xmax_hours )))
)
bar_width = 1.0/24 # one hour
if plot_trends:
# remove any linear trend
coeff = np.polyfit(xdata, ydata, background_order )
trend = np.polyval(coeff, xdata)
fig, ax = plt.subplots(figsize=[12,10])
date_formatter = mdates.DateFormatter('%m-%d %H:%M')
ax.tick_params(axis="x", rotation=90)
ax.xaxis.set_major_locator(mdates.HourLocator(interval = 6))
ax.xaxis.set_major_formatter(date_formatter)
ax.xaxis.set_minor_locator(mdates.HourLocator(interval = 1))
fig.subplots_adjust(bottom=0.2)
ax3 = ax.twinx()
ax3.set_axis_off()
ax3.set_zorder(1)
ax3.bar( peakt, [1,1,1] , width=bar_width, color="lightgrey")
ax.set_xlim(left=plotx_range ) # start at a round hour
ax.set_zorder(2)
ax.patch.set_visible(False)
ax.plot(tdata, ydata, color="blue", marker="+", linewidth=1)
if plot_trends:
ax.plot(tdata, trend, color="lightgrey", linewidth=1)
ax2=ax.twinx()
ax2.set_zorder(3)
ax2.plot(tdata, ydata - trend, color="lightblue", linewidth=3)
# save the plot as a file
filename = "./Hunga_Tonga_" + "all" + ".png"
fig.savefig( filename, format='png', dpi=100, bbox_inches='tight')
#plt.show()
hunga_tonga = (-20.5452074472518, -175.38715105641674)
#speed_of_sound = 0.32 # km/s
speed_of_sound = 0.318 # km/s - derived value from fitting first and 2nd pulses to several points across Australia.
eruption_time = 1642220085 # 2022-01-15 04:14:45 UTC
earth_circumference=40040 # near enough
hour_span = 6
distance = distance.distance(hunga_tonga, home).km
travel_time = distance / speed_of_sound
arrival_time = eruption_time + travel_time
once_around = earth_circumference / speed_of_sound
opposite_wave_time = eruption_time + once_around - travel_time
return_wave_time = eruption_time + once_around + travel_time
print("distance to eruption %0.1f km arrival at %.0f (%s)" %
(distance, arrival_time, time.strftime("%Y-%m-%d %H:%M:%S", time.localtime(arrival_time))))
print(" opposite pulse arrival at %.0f (%s)" %
( opposite_wave_time , time.strftime("%Y-%m-%d %H:%M:%S", time.localtime(opposite_wave_time ))))
print(" second time around pulse arrival at %.0f (%s)" %
( return_wave_time , time.strftime("%Y-%m-%d %H:%M:%S", time.localtime(return_wave_time ))))
if DB == "sqlite3":
connection = sqlite3.connect('weewx.sdb')
else:
try:
connection = MysqlConnector.connect( user=db_user, password=db_pwd,
host='127.0.0.1', database=db_name)
except MysqlConnector.Error as err:
if err.errno == msqlerrcode.ER_ACCESS_DENIED_ERROR:
print("Bad user name or password")
elif err.errno == msqlerrcode.ER_BAD_DB_ERROR:
print("Database does not exist")
else:
print(err)
sys.exit( 2)
cursor = connection.cursor()
plot_burst( cursor, arrival_time, hour_span, "primary" )
plot_burst( cursor, opposite_wave_time, hour_span, "secondary" )
plot_burst( cursor, return_wave_time, hour_span, "return" )
plot_all( cursor, eruption_time, arrival_time, opposite_wave_time, return_wave_time, hour_span )
connection.close()
