Files
weather-station/sensors/main.py

578 lines
24 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
import time, math, board, smbus2, serial, psutil, pynmea2, bme680
import adafruit_ltr390, threading, os, statistics, requests
from datetime import datetime, timezone
from collections import deque
from dotenv import load_dotenv
from pathlib import Path
BASE_DIR = Path(__file__).resolve().parent.parent
def cfg():
load_dotenv(BASE_DIR / '.env', override=True)
load_dotenv(BASE_DIR / '.env.local', override=True)
return {
'LATITUDE': float(os.getenv('LATITUDE', '0.0')),
'LONGITUDE': float(os.getenv('LONGITUDE', '0.0')),
'ALTITUDE': float(os.getenv('ALTITUDE', '0.0')),
'DB_HOST': os.getenv('DB_HOST', 'http://localhost:8428'),
'GPS_PORT': os.getenv('GPS_PORT', '/dev/ttyS0'),
'GPS_BAUD': int( os.getenv('GPS_BAUD', '9600')),
'GAS_REFERENCE': int( os.getenv('GAS_REFERENCE', '250000')),
'PRESSURE_TREND_WINDOW': int( os.getenv('PRESSURE_TREND_WINDOW', '60')),
'LUNA_CAL_SAMPLES': int( os.getenv('LUNA_CAL_SAMPLES', '20')),
'LUNA_MAX_DIST_CM': int( os.getenv('LUNA_MAX_DIST_CM', '800')),
'LUNA_MIN_DEPTH_CM': int( os.getenv('LUNA_MIN_DEPTH_CM', '2')),
'AS3935_ADDR': int( os.getenv('AS3935_ADDR', '0x03'), 16),
'AS3935_NOISE_FLOOR': int( os.getenv('AS3935_NOISE_FLOOR', '5')),
'AS3935_WATCHDOG': int( os.getenv('AS3935_WATCHDOG', '3')),
'AS3935_SPIKE_REJ': int( os.getenv('AS3935_SPIKE_REJ', '7')),
'LOOP_INTERVAL': float( os.getenv('LOOP_INTERVAL', '1.0')),
'FLOOD_TEMP_THRESHOLD': float( os.getenv('FLOOD_TEMP_THRESHOLD', '10.0')),
'SLUSH_TEMP_THRESHOLD': float( os.getenv('SLUSH_TEMP_THRESHOLD', '0.0')),
'SNOW_STRENGTH_MIN': int( os.getenv('SNOW_STRENGTH_MIN', '700')),
'SEISMIC_NOISE_FLOOR': float( os.getenv('SEISMIC_NOISE_FLOOR', '0.02')),
'TEMP_CALIBRATION': float( os.getenv('TEMP_CALIBRATION', '0.0')),
'VISIBILITY_MAX': float( os.getenv('VISIBILITY_MAX', '50.0')),
'CLOUD_SUN_MIN_ELEV': float( os.getenv('CLOUD_SUN_MIN_ELEV', '5.0')),
'CLOUD_LUX_CLEAR_SKY': float( os.getenv('CLOUD_LUX_CLEAR_SKY', '100000.0')),
'CLOUD_HUMIDITY_WEIGHT': float( os.getenv('CLOUD_HUMIDITY_WEIGHT', '0.3')),
}
# ── VictoriaMetrics ───────────────────────────────────────────────────────────
def write(c, fields, collection='weather'):
field_str = ','.join(
f"{k}={'true' if v is True else 'false' if v is False else v}"
for k, v in fields.items() if v is not None
)
if not field_str: return
requests.post(f"{c['DB_HOST']}/write", data=f"{collection} {field_str}", timeout=2)
# ── Helpers ───────────────────────────────────────────────────────────────────
def try_init(fn, name):
try:
r = fn()
print(f"{name}")
return r
except Exception as e:
print(f"{name}: {e}")
return None
def try_read(fn):
try: return fn()
except: return None
def s16(v): return v - 65536 if v > 32767 else v
try: bus = smbus2.SMBus(1)
except: bus = None; print("[I2C] SMBus init failed")
try: i2c = board.I2C()
except: i2c = None; print("[I2C] board.I2C init failed")
# ── BME680 ────────────────────────────────────────────────────────────────────
pressure_buffer = deque(maxlen=60)
def init_bme():
sensor = bme680.BME680(bme680.I2C_ADDR_SECONDARY)
sensor.set_humidity_oversample(bme680.OS_2X)
sensor.set_pressure_oversample(bme680.OS_4X)
sensor.set_temperature_oversample(bme680.OS_8X)
sensor.set_filter(bme680.FILTER_SIZE_3)
sensor.set_gas_status(bme680.ENABLE_GAS_MEAS)
sensor.set_gas_heater_temperature(320)
sensor.set_gas_heater_duration(150)
sensor.select_gas_heater_profile(0)
return sensor
def heat_index(tc, rh):
tf = tc * 9/5 + 32
if tf < 80: return round(tc, 2)
hi = (-42.379 + 2.04901523*tf + 10.14333127*rh
- 0.22475541*tf*rh - 0.00683783*tf**2
- 0.05481717*rh**2 + 0.00122874*tf**2*rh
+ 0.00085282*tf*rh**2 - 0.00000199*tf**2*rh**2)
return round((hi - 32) * 5/9, 2)
def absolute_humidity(tc, rh):
return round((6.112 * math.exp((17.67*tc)/(tc+243.5)) * rh * 2.1674) / (273.15+tc), 3)
def vpd(tc, rh):
es = 0.6108 * math.exp((17.27*tc)/(tc+237.3))
return round(es - (rh/100)*es, 3)
def gas_to_aqi(gas_ohms, humidity, gas_reference):
if gas_ohms is None: return None
gas_score = min(gas_ohms / gas_reference, 1.0) * 75
hum_score = 25 - abs(humidity - 40) * 0.5
quality = min(max(gas_score + hum_score, 0), 100)
aqi = round((1 - quality / 100) * 500)
return aqi
def sea_level_pressure(pressure_hpa, alt_m):
if alt_m is None: return None
return round(pressure_hpa / math.pow(1 - (alt_m / 44330.0), 5.255), 2)
def get_pressure_trend(window):
if pressure_buffer.maxlen != window:
pressure_buffer.__init__(window)
if len(pressure_buffer) < 2: return None
delta = pressure_buffer[-1] - pressure_buffer[0]
return round(delta, 2)
def read_bme(c, sensor, alt_m=None):
if not sensor: return None, None, None
try:
if not sensor.get_sensor_data(): return None, None, None
t, rh, p = sensor.data.temperature, sensor.data.humidity, sensor.data.pressure
gas = round(sensor.data.gas_resistance, 0) if sensor.data.heat_stable else None
aqi_score = gas_to_aqi(gas, rh, c['GAS_REFERENCE'])
pressure_buffer.append(p)
p_rate = get_pressure_trend(c['PRESSURE_TREND_WINDOW'])
dp = round(t - (100 - rh) / 5.0, 2)
write(c, {
'temperature': round(t + c['TEMP_CALIBRATION'], 2),
'humidity': round(rh, 2),
'dew_point': dp,
'pressure_hpa': round(p, 2),
'pressure_slp': sea_level_pressure(p, alt_m),
'pressure_rate': p_rate,
'humidity_abs': absolute_humidity(t, rh),
'vapor_pressure_deficit': vpd(t, rh),
'heat_index': heat_index(t, rh),
'air_quality_ohms': gas,
'air_quality': aqi_score,
})
return t, rh, aqi_score
except Exception as e:
print(f"[BME680] {e}")
return None, None, None
# ── MPU6050 ───────────────────────────────────────────────────────────────────
# ── MPU6050 ───────────────────────────────────────────────────────────────────
MPU_ADDR = 0x68
_mpu_peak = None
_mpu_lock = threading.Lock()
_mpu_grav = None # slow-tracked gravity vector [x, y, z]
GRAV_ALPHA = 0.05 # EMA rate while still
STILL_SPREAD = 0.03 # g, max peak-to-peak per axis to count as "still"
_still_buf = deque(maxlen=50) # ~0.5s @ 100Hz
def read_accel():
d = try_read(lambda: bus.read_i2c_block_data(MPU_ADDR, 0x3B, 6))
if not d: return None
return [
s16((d[0] << 8) | d[1]) / 16384.0,
s16((d[2] << 8) | d[3]) / 16384.0,
s16((d[4] << 8) | d[5]) / 16384.0,
]
def rotation_to_z(g):
"""Rodrigues rotation mapping gravity vector g -> +Z. Returns (R, |g|)."""
mag = math.sqrt(sum(v * v for v in g))
u = [v / mag for v in g]
kx, ky = u[1], -u[0] # axis = u × z
s = math.sqrt(kx * kx + ky * ky) # sin(theta)
c = u[2] # cos(theta)
if s < 1e-8:
return [[1, 0, 0], [0, c, 0], [0, 0, c]], mag
kx, ky = kx / s, ky / s
v = 1 - c
return [
[c + kx * kx * v, kx * ky * v, ky * s],
[ky * kx * v, c + ky * ky * v, -kx * s],
[-ky * s, kx * s, c ],
], mag
def is_still(buf):
if len(buf) < buf.maxlen: return False
return all(
max(s[i] for s in buf) - min(s[i] for s in buf) < STILL_SPREAD
for i in range(3)
)
def calibrate_mpu(samples=200):
if not bus: return None
print(" Calibrating MPU6050, keep still...")
sums, n = [0.0, 0.0, 0.0], 0
for _ in range(samples):
a = read_accel()
if a:
sums = [s + v for s, v in zip(sums, a)]
n += 1
time.sleep(0.005)
if n == 0: return None
g = [s / n for s in sums]
mag = math.sqrt(sum(v * v for v in g))
if not 0.5 < mag < 1.5:
print(f" ❌ MPU6050: bad gravity magnitude {mag:.2f}g")
return None
print(f" ✅ MPU6050: gravity {mag:.3f}g on ({g[0]:.2f}, {g[1]:.2f}, {g[2]:.2f}) → mapped to Z")
return g
def mpu_loop():
global _mpu_peak, _mpu_grav
while True:
if bus and _mpu_grav:
a = read_accel()
if a:
_still_buf.append(a)
R, g = rotation_to_z(_mpu_grav)
ax = R[0][0] * a[0] + R[0][1] * a[1] + R[0][2] * a[2]
ay = R[1][0] * a[0] + R[1][1] * a[1] + R[1][2] * a[2]
az = R[2][0] * a[0] + R[2][1] * a[1] + R[2][2] * a[2] - g
mag = math.sqrt(ax**2 + ay**2 + az**2)
# Re-track gravity whenever the sensor is *still* —
# absorbs bias shifts of any size, blocks shakes
if is_still(_still_buf):
_mpu_grav = [p + (v - p) * GRAV_ALPHA for p, v in zip(_mpu_grav, a)]
with _mpu_lock:
if _mpu_peak is None or mag > _mpu_peak['seismic_magnitude']:
_mpu_peak = {
'seismic_x': round(ax, 3),
'seismic_y': round(ay, 3),
'seismic_z': round(az, 3),
'seismic_magnitude': round(mag, 4),
}
time.sleep(0.01)
def flush_mpu(c):
global _mpu_peak
with _mpu_lock:
data, _mpu_peak = _mpu_peak, None
if data and data['seismic_magnitude'] > c['SEISMIC_NOISE_FLOOR']:
write(c, data)
else:
write(c, {'seismic_x': 0.0, 'seismic_y': 0.0, 'seismic_z': 0.0, 'seismic_magnitude': 0.0})
# ── QMC5883L ──────────────────────────────────────────────────────────────────
def read_compass(c):
if not bus: return
try:
addr = 0x1E
try_read(lambda: bus.write_byte_data(addr, 0x09, 0x1D))
d = try_read(lambda: bus.read_i2c_block_data(addr, 0x00, 6))
if not d: return
x = s16((d[1] << 8) | d[0])
y = s16((d[3] << 8) | d[2])
z = s16((d[5] << 8) | d[4])
heading = math.degrees(math.atan2(y, x))
if heading < 0: heading += 360
write(c, {'compass_x': x, 'compass_y': y, 'compass_z': z, 'heading': round(heading, 1)})
except Exception as e:
print(f"[QMC5883L] {e}")
# ── LTR390 ───────────────────────────────────────────────────────────────────
uv_dose_mj = 0.0
uv_dose_date = datetime.now().date()
last_uv_time = time.time()
dli_lux_acc = 0.0
dli_date = datetime.now().date()
last_lux_time = time.time()
daylight_start = None
DAYLIGHT_LUX_THRESHOLD = 50
def sun_elevation(lat, lon):
now = datetime.now(timezone.utc)
doy = now.timetuple().tm_yday
hour_ut = now.hour + now.minute / 60 + now.second / 3600
decl = math.radians(23.45 * math.sin(math.radians((360 / 365) * (doy - 81))))
lstm = 15 * round(lon / 15)
eot = (9.87 * math.sin(math.radians(2 * (360/365) * (doy - 81)))
- 7.53 * math.cos(math.radians((360/365) * (doy - 81)))
- 1.5 * math.sin(math.radians((360/365) * (doy - 81))))
solar_time = hour_ut + (lon - lstm) / 15 + eot / 60
ha = math.radians(15 * (solar_time - 12))
lat_r = math.radians(lat)
elev = math.degrees(math.asin(
math.sin(lat_r) * math.sin(decl) +
math.cos(lat_r) * math.cos(decl) * math.cos(ha)
))
return elev
def estimate_cloud_cover(c, lux, temp, dew_point, lat, lon):
elev = sun_elevation(lat, lon)
temp_dp_spread = temp - dew_point
hum_cloud = max(0.0, min(100.0, (1 - temp_dp_spread / 20.0) * 100))
if elev < c['CLOUD_SUN_MIN_ELEV']:
return round(hum_cloud, 1), round(elev, 2)
theoretical = max(c['CLOUD_LUX_CLEAR_SKY'] * math.sin(math.radians(elev)), 1.0)
solar_ratio = min(lux / theoretical, 1.0)
solar_cloud = (1.0 - solar_ratio) * 100
w = c['CLOUD_HUMIDITY_WEIGHT']
blended = (1 - w) * solar_cloud + w * hum_cloud
return round(min(max(blended, 0.0), 100.0), 1), round(elev, 2)
def estimate_visibility(c, rh, aqi):
rh_clamped = min(rh, 99.5)
beta_rh = 0.000146 * math.exp(0.06 * rh_clamped)
aqi_factor = 1.0 + ((aqi or 0) / 500.0) * 0.5
extinction = beta_rh * aqi_factor
return round(min(3.912 / extinction, c['VISIBILITY_MAX']), 2)
def read_ltr(c, sensor, temp=None, dew_point=None, rh=None, aqi=None, lat=0.0, lon=0.0):
global uv_dose_mj, uv_dose_date, last_uv_time
global dli_lux_acc, dli_date, last_lux_time, daylight_start
if not sensor: return
try:
lux = round(sensor.lux, 2)
uvi = round(sensor.uvi, 2)
except Exception as e:
print(f"[LTR390] {e}")
return
now = time.time()
today = datetime.now().date()
if today != uv_dose_date: uv_dose_mj, uv_dose_date = 0.0, today
if today != dli_date: dli_lux_acc, dli_date, daylight_start = 0.0, today, None
uv_dose_mj += uvi * 25 * (now - last_uv_time)
dli_lux_acc += lux * (now - last_lux_time)
last_uv_time = now
last_lux_time = now
if lux >= DAYLIGHT_LUX_THRESHOLD and daylight_start is None:
daylight_start = now
fields = {
'lux': lux,
'uv_index': uvi,
'solar_wm2': round(lux / 120, 2),
'uv_dose': round(uv_dose_mj, 2),
'daily_light_integral': round(dli_lux_acc / 54 / 1_000_000, 4),
}
if temp is not None and dew_point is not None:
cloud, sun_elev = estimate_cloud_cover(c, lux, temp, dew_point, lat, lon)
fields['clouds'] = cloud
fields['sun_elevation'] = sun_elev
if rh is not None:
fields['visibility'] = estimate_visibility(c, rh, aqi)
write(c, fields)
# ── TF-Luna ───────────────────────────────────────────────────────────────────
luna_baseline_cm = None
def calibrate_luna(c):
global luna_baseline_cm
if not bus: print(" ❌ TF-Luna: no I2C bus"); return
print(" Calibrating TF-Luna...")
readings = []
for _ in range(c['LUNA_CAL_SAMPLES']):
d = try_read(lambda: bus.read_i2c_block_data(0x10, 0x00, 6))
if d:
dist, strength = d[0] | (d[1] << 8), d[2] | (d[3] << 8)
if strength >= 100 and 0 < dist <= c['LUNA_MAX_DIST_CM']:
readings.append(dist)
time.sleep(0.1)
if readings:
luna_baseline_cm = round(sum(readings) / len(readings), 1)
print(f" ✅ TF-Luna baseline: {luna_baseline_cm} cm")
else:
print(" ❌ TF-Luna: calibration failed")
def read_luna(c):
if not bus: return
try:
d = try_read(lambda: bus.read_i2c_block_data(0x10, 0x00, 6))
if not d: return
dist = d[0] | (d[1] << 8)
strength = d[2] | (d[3] << 8)
if strength < 100 or dist <= 0 or dist > c['LUNA_MAX_DIST_CM']: return
fields = {'ground_distance': dist, 'ground_calibration': luna_baseline_cm, 'lidar_strength': strength}
if luna_baseline_cm and (luna_baseline_cm - dist) >= c['LUNA_MIN_DEPTH_CM']:
depth = round(luna_baseline_cm - dist, 1)
fields['accumulation'] = depth
else:
fields['accumulation'] = 0
write(c, fields)
except Exception as e:
print(f"[TF-Luna] {e}")
# ── AS3935 ────────────────────────────────────────────────────────────────────
lightning_distances = deque(maxlen=10)
lightning_count = 0
lightning_window_start = time.time()
def init_as3935(c):
if not bus: return
try:
addr = c['AS3935_ADDR']
bus.write_byte_data(addr, 0x01, (c['AS3935_NOISE_FLOOR'] << 4) | c['AS3935_WATCHDOG'])
bus.write_byte_data(addr, 0x02, c['AS3935_SPIKE_REJ'])
bus.write_byte_data(addr, 0x00, 0x24)
print(" ✅ AS3935")
except Exception as e:
print(f" ❌ AS3935: {e}")
def read_as3935(c):
global lightning_count, lightning_window_start
if not bus: return
try:
addr = c['AS3935_ADDR']
data = try_read(lambda: bus.read_i2c_block_data(addr, 0x00, 9))
if not data: return
interrupt = data[3] & 0x0F
energy = ((data[5] & 0x1F) << 16) | (data[4] << 8) | data[3]
try_read(lambda: bus.read_byte_data(addr, 0x03))
if interrupt not in (0x04, 0x08) or energy < 1000: return
dist = data[6] & 0x3F
lightning_distances.append(dist)
lightning_count += 1
elapsed = time.time() - lightning_window_start
if elapsed >= 3600:
lightning_count, lightning_window_start, elapsed = 1, time.time(), 1
trend = 'Stationary'
if len(lightning_distances) >= 3:
avg_early = sum(list(lightning_distances)[:3]) / 3
avg_late = sum(list(lightning_distances)[-3:]) / 3
if avg_late < avg_early - 2: trend = 'Approaching'
elif avg_late > avg_early + 2: trend = 'Retreating'
write(c, {
'lightning_distance': dist,
'lightning_rate': round(lightning_count / (elapsed / 3600), 1),
'storm_trend': trend,
})
except Exception as e:
print(f"[AS3935] {e}")
# ── GPS ───────────────────────────────────────────────────────────────────────
gps_cache = {}
def read_gps(c):
try:
ser = serial.Serial(c['GPS_PORT'], c['GPS_BAUD'], timeout=1)
line = ser.readline().decode('ascii', errors='replace').strip()
ser.close()
if 'GGA' in line:
msg = pynmea2.parse(line)
gps_cache.update({
'latitude': msg.latitude,
'longitude': msg.longitude,
'altitude': msg.altitude,
'gps_satellites': int(msg.num_sats),
})
if gps_cache:
write(c, {k: v for k, v in gps_cache.items() if v is not None})
except Exception as e:
print(f"[GPS] {e}")
return gps_cache
# ── Rain ──────────────────────────────────────────────────────────────────────
def read_rain(c):
write(c, {'precipitation': 0, 'raining': False})
# ── Wind ──────────────────────────────────────────────────────────────────────
def read_wind(c):
write(c, {'wind_direction': 0, 'wind_speed': 0, 'wind_gusts': 0})
# ── System ────────────────────────────────────────────────────────────────────
def read_system(c):
try:
cpu_temp = None
for key in ('cpu_thermal', 'coretemp', 'k10temp', 'acpitz'):
if key in (temps := psutil.sensors_temperatures()):
cpu_temp = round(temps[key][0].current, 1)
break
cpu_freq = psutil.cpu_freq()
mem = psutil.virtual_memory()
disk = psutil.disk_usage('/')
write(c, {
'cpu_temp': cpu_temp,
'cpu_usage': round(psutil.cpu_percent(interval=None), 1),
'cpu_freq_mhz': round(cpu_freq.current, 1) if cpu_freq else None,
'mem_used_mb': round(mem.used / 1024**2, 1),
'mem_total_mb': round(mem.total / 1024**2, 1),
'mem_percent': mem.percent,
'disk_used_gb': round(disk.used / 1024**3, 2),
'disk_total_gb': round(disk.total / 1024**3, 2),
'disk_percent': disk.percent,
}, collection='system')
except Exception as e:
print(f"[System] {e}")
# ── Main ──────────────────────────────────────────────────────────────────────
if __name__ == '__main__':
print("\n🌦 Weather Station — Initialising...\n")
c = cfg()
bme_sensor = try_init(init_bme, 'BME680')
ltr_sensor = try_init(lambda: adafruit_ltr390.LTR390(i2c) if i2c else (_ for _ in ()).throw(Exception('no i2c')), 'LTR390')
if bus:
try_read(lambda: bus.write_byte_data(MPU_ADDR, 0x6B, 0))
time.sleep(0.1)
_mpu_grav = try_init(calibrate_mpu, 'MPU6050')
else:
_mpu_grav = None
print(" ❌ MPU6050: no I2C bus")
threading.Thread(target=mpu_loop, daemon=True).start()
init_as3935(c)
calibrate_luna(c)
print("\n🚀 Starting sensor loop...\n")
last_temp = None
last_rh = None
last_aqi = None
while True:
c = cfg()
gps = read_gps(c)
lat = gps.get('latitude') or c['LATITUDE']
lon = gps.get('longitude') or c['LONGITUDE']
alt = gps.get('altitude') or c['ALTITUDE']
result = read_bme(c, bme_sensor, alt_m=alt)
if result and result[0] is not None:
last_temp, last_rh, last_aqi = result
dp = None
if last_temp is not None and last_rh is not None:
dp = round(last_temp - (100 - last_rh) / 5.0, 2)
flush_mpu(c)
read_compass(c)
read_ltr(c, ltr_sensor,
temp=last_temp, dew_point=dp,
rh=last_rh, aqi=last_aqi,
lat=lat, lon=lon)
read_luna(c)
read_as3935(c)
read_rain(c)
read_wind(c)
read_system(c)
time.sleep(c['LOOP_INTERVAL'])