Weather Sonde automated prediction
So, I've gotten the bug to chase weather balloons. The balloons in my region are launched by the National Weather Service from a location adjoining the airport in Buffalo, New York. Only a small percentage of the balloons make it southeast enough to be within reasonable chase distance from my home near Ithaca, New York, 110 air miles away. I ran manual predictions every day, to see how the balloons might progress, but that got tedious. I decided some automation was in order.
CUFS Prediction Software
To begin, I picked up the CUFS Prediction Wrapper for python. It provided tools to download the latest NWS prediction data as well as tools to run predictions. I run a Raspberry Pi at home, so I set up cron jobs to load the data each morning, and to run predictions for the next 7 days.Downloading Wind Data
I set up a cron job to download the wind data at 3am each day:0 3 * * * /home/mqh1/bin/wind_grabber.sh > /home/mqh1/balloons/cusf/logs/windgrabber.log 2>& 1
The contents of the "wind_grabber.sh' script are as follows:
python \
/home/mqh1/balloons/cusf/cusf_predictor_wrapper/apps/get_wind_data.py \
--lat=42 --lon=-77 --latdelta=2 --londelta=4 -f 180 \
-m 0p25_1hr -o /tmp/gfs/
The lat/long as well as latedelta/londelta parameters are sufficient to include most of the likely flight locations of a weather balloon out of Buffalo. I cache the results in /tmp/gfs on the pi. "-f 180" is looking for the next 180 hours of predictions (7 1/2 days).
The result is that I have current wind data automatically downloaded on a daily basis that should cover the areas of interest for any NWS flights. Since I do my own High Altitude Balloon launches, as well as pico flights, I don't feel too guilty about downloading this data daily.
Todo: I'd like to find out when the NWS uploads the data, and tune my download and daily prediction schedule to use the latest data.
Running predictions
Typical flights
I was curious what the best parameters would be to choose for ascent rate / descent rate / and burst altitude for the NWS balloon. Fortunately, two of the three are provided on a regular basis by the NWS. They have an index with the codes for all the launch locations. I was able to download a big file with all of the flight data from the last two years of flights. I stuck it in a spreadsheet, and did some number crunching. I had to eliminate a few flights that were clearly outliers. They had really screwy ascent rates or burst altitudes compared to the rest. I removed the anomalies that had glaring differences, and ran averages on the remainder. I got the following averages.Average burst altitude: 31,420 meters
Average ascent rate: 5.28 m/sec
Unfortunately, the NWS doesn't log the descent data, so I don't have a good measure of the average descent rate of the sondes. Based on conversations with folks who do this, it seems like most people go with 5m/s. It can vary widely however.
Todo: If I get a fully automated station running, I'd like to get better data about descent rates and decide how to handle it. I might use an average, or I might do predictions on each end of the range.
Automated predictions
The CUSF software makes it really easy to run predictions. I run them after I download the GFS data. I haven't measured the typical download rate yet, so for now, I run them 90 minutes after the wind data download.30 4 * * * /home/mqh1/bin/sonde_predict.py > /home/mqh1/balloons/cusf/logs/sonde_predict.log 2>& 1
The prediction script is a bit more complicated than the download script. I started from the sonde_predict.py script provided by the CUSF wrapper and made three principal mods.
- Determining which predictions are coming in "near' Ithaac
- Generating two separate KML files. One for "nearby" and one for "all" predictions.
- Emailing results.
Determining nearby predictions
After the prediction runs, I take the final landing coordinates and determine whether they are "nearby" Ithaca. I do this with a Point in Polygon (PIP) algorithm.
I used Google Earth to draw a bounding box covering an area around Ithaca in which I was willing to search for balloons (Add / Polygon). I then exported the polygon to a KML (right-click, Save Place As), I was then able to grab the coordinates from the KML file, and just copy/paste them to create an array in the code. Note, I had edit the data to reverse the lat/long coming out of the KML file to make them more friendly. In retrospect, I could have left them reversed, and just reversed the order in the code below. Still, I think this is more "human friendly" to see them in the typical order.
I used Google Earth to draw a bounding box covering an area around Ithaca in which I was willing to search for balloons (Add / Polygon). I then exported the polygon to a KML (right-click, Save Place As), I was then able to grab the coordinates from the KML file, and just copy/paste them to create an array in the code. Note, I had edit the data to reverse the lat/long coming out of the KML file to make them more friendly. In retrospect, I could have left them reversed, and just reversed the order in the code below. Still, I think this is more "human friendly" to see them in the typical order.
# Coordinate around Ithaca (larger grid)
nearby_coords = [
(42.5894692373013, -77.13750089643533),
(42.13873798563936, -77.11683620792473),
(42.1355764571338, -75.91998110781731),
(42.62563890419801, -75.89559117314741)
]
The coords of the prediction of the payload are in the script in the flight_path[] array, so we grab the last position.
finalcoords = flight_path[-1]
latlong = (finalcoords[1], finalcoords[2])
landingsite = Point(latlong)
Then we can determine whether it's in the "nearby_coords" polygon we've established:
if landingsite.within(nearby_poly):
Creating a second KML file
If we're within the "nearby_poly", I add the KML plot to a separate list. I modified the provided script to generate two KML files. One for "all' predictions, and one for "nearby" predictions. I write both of the predictions into an HTML directory for the web server I run on the Pi. That makes it trivial to download them to my laptop to inspect them with Google Earth.
Emailing results
Finally, I modified the prediction script to generate an email summarizing the number of nearby predictions, and containing HTML links to both the nearby and full KML files. Each morning, I receive an email that tells me if I have nearby predictions. If I do, I can click the link to see the predictions for the coming week that land near Ithaca.
In the image below, I highlighted the bounding box for my "nearby" predictions.
Todo: I'd like to find a way to provide a link into google maps which will show the 2D representation of the prediction by just clicking a link, rather than having to download a KML file and launch Google Earth.
Sample Code
Disclaimers: Shipped by weight not volume. Some settling may occur during shipping. Not actual size. Consult your computer professional before running. If this code causes and erection lasting more than 4 hours, consult a physician.. no, seriously!
#!/usr/bin/env python
#
# Project Horus
# CUSF Standalone Predictor Python Wrapper
# Radiosonde Launch Predictor
# Copyright 2017 Mark Jessop <vk5qi@rfhead.net>
#
# In this example we run predictions for a radiosonde launch from Adelaide Airport
# for every launch for the next 7 days, and write the tracks out to a KML file.
#
# For this script to work correctly, we need a weeks worth of GFS data available within the gfs directory.
# This can be gathered using get_wind_data.py (or scripted with wind_grabber.sh), using :
# python get_wind_data.py --lat=-33 --lon=139 --latdelta=10 --londelta=10 -f 168 -m 0p25_1hr -o gfs
# with lat/lon parameters chaged as appropriate.
#
# Mods to this script written by Mike Hojnowski / KD2EAT. It is provided for informational purposes. It will
# not run in a generic environment without modification. In particular, look for the email addresses and for
# the home directories sprinkled in paths all over.
#
import smtplib
import time
import fastkml
import datetime
from dateutil.parser import parse
from cusfpredict.predict import Predictor
from pygeoif.geometry import Point, LineString
from cusfpredict.utils import *
from shapely.geometry import Point, Polygon
#Email Variables
SMTP_SERVER = 'smtp.gmail.com' #Email Server (don't change!)
SMTP_PORT = 587 #Server Port (don't change!)
GMAIL_USERNAME = 'somebody@gmail.com' #change this to match your gmail account
GMAIL_PASSWORD = 'ImDumbButNotTHATDumb' #change this to match your gmail password
SENDTO = 'ThisCouldBeYou@gmail.com' #change to the recipient of the emails
## Coordinate around Ithaca (smaller grid)
#nearby_coords = [
# (42.44387782645627, -77.01612142816262),
# (42.2569086411454,-76.99969683890008),
# (42.3432446683892, -76.17009222210787),
# (42.57580936472674, -76.28671583838107),
# ]
#
# Coordinate around Ithaca (larger grid)
nearby_coords = [
(42.5894692373013, -77.13750089643533),
(42.13873798563936, -77.11683620792473),
(42.1355764571338, -75.91998110781731),
(42.62563890419801, -75.89559117314741)
]
nearby_poly = Polygon(nearby_coords)
# BIG polygon for testing. No predictions were near Ithaca the day I was coding, so I made a bigger polygon.
huge_coords = [
(43.44387782645627, -78.01612142816262),
(41.2569086411454,-78.99969683890008),
(41.3432446683892, -75.17009222210787),
(43.57580936472674, -75.28671583838107),
]
# Just substitude "huge_poly" for "nearby_poly" in the code below, if you need to test during unfavorable predictions.
huge_poly = Polygon(huge_coords)
# Predictor Parameters
PRED_BINARY = "/home/mqh1/balloons/cusf/cusf_predictor_wrapper/apps/pred"
GFS_PATH = "/tmp/gfs"
NEARBY_OUTPUT_FILE = "/home/mqh1/public_html/nearby_predictions.kml"
ALL_OUTPUT_FILE = "/home/mqh1/public_html/all_predictions.kml"
# Launch Parameters - Update as appropriate for your launch site
# Launch parameters based on NWS data from 01/01/2018 - 12/22/2019 averages
LAUNCH_LAT = 42.9414
LAUNCH_LON = -78.7193
LAUNCH_ALT = (218+500)/2.0; # Barometric and GPS varied. I took the average.
ASCENT_RATE = 5.26
DESCENT_RATE = 5.0
BURST_ALT = 31420
# Set the launch time to the current UTC day, but set the hours to the 12Z sonde
current_day = datetime.datetime.utcnow()
LAUNCH_TIME = datetime.datetime(current_day.year, current_day.month, current_day.day, 11, 15)
# Parameter Variations
# These can all be left at zero, or you can add a range of delta values
launch_time_variations = range(0,168,12) # Every 12 hours from now until 7 days time.
# A list to store prediction results
nearby_predictions = []
all_predictions = []
nearby = 0
# Create the predictor object.
pred = Predictor(bin_path=PRED_BINARY, gfs_path=GFS_PATH)
# Iterate through the range of launch times set above
for _delta_time in launch_time_variations:
# Calculate the launch time for the current prediction.
_launch_time = LAUNCH_TIME + datetime.timedelta(seconds=_delta_time*3600)
print ('Running prediction for ' + str(_launch_time))
# Run the prediction
flight_path = pred.predict(
launch_lat=LAUNCH_LAT,
launch_lon=LAUNCH_LON,
launch_alt=LAUNCH_ALT,
ascent_rate=ASCENT_RATE,
descent_rate=DESCENT_RATE,
burst_alt=BURST_ALT,
launch_time=_launch_time)
# If we only get a single entry in the output array, it means we don't have any wind data for this time
# Continue on to the next launch time.
if len(flight_path) == 1:
continue
#print(flight_path[-1])
coords = flight_path[-1]
latlong = (coords[1], coords[2])
#print latlong
landingsite = Point(latlong)
#print landingsite.within(nearby_poly)
#print landingsite.within(huge_poly)
# Generate a descriptive comment for the track and placemark.
pred_time_string = _launch_time.strftime("%Y%m%d-%H%M")
pred_comment = "%s %.1f/%.1f/%.1f" % (pred_time_string, ASCENT_RATE, BURST_ALT, DESCENT_RATE)
# Add the track and placemark to our list of predictions
all_predictions.append(flight_path_to_geometry(flight_path, comment=pred_comment))
all_predictions.append(flight_path_landing_placemark(flight_path, comment=pred_comment))
if landingsite.within(nearby_poly):
# If nearby, add the placemark to our nearby predictions list
nearby_predictions.append(flight_path_to_geometry(flight_path, comment=pred_comment))
nearby_predictions.append(flight_path_landing_placemark(flight_path, comment=pred_comment))
print("Prediction Succeeds: %s" % pred_comment)
nearby += 1
# Write out the prediction data to the KML file
kml_comment = "Sonde Predictions - %s" % gfs_model_age()
write_flight_path_kml(all_predictions, filename=ALL_OUTPUT_FILE, comment=kml_comment)
write_flight_path_kml(nearby_predictions, filename=NEARBY_OUTPUT_FILE, comment=kml_comment)
###################
# Email portion
class Emailer:
def sendmail(self, recipient, subject, content):
#Create Headers
headers = ["From: " + GMAIL_USERNAME, "Subject: " + subject, "To: " + recipient,
"MIME-Version: 1.0", "Content-Type: text/html"]
headers = "\r\n".join(headers)
#Connect to Gmail Server
session = smtplib.SMTP(SMTP_SERVER, SMTP_PORT)
session.ehlo()
session.starttls()
session.ehlo()
#Login to Gmail
session.login(GMAIL_USERNAME, GMAIL_PASSWORD)
#Send Email & Exit
session.sendmail(GMAIL_USERNAME, recipient, headers + "\r\n\r\n" + content)
session.quit
sender = Emailer()
sendTo = SENDTO
emailSubject = "Sonde predictions nearby: " + str(nearby)
emailContent = str(nearby) + ' predictions have been detected near Ithaca.<br><br>' +\
'\nNearby Predictions: <a href=http://thehojos.com/~mqh1/nearby_predictions.kml>http://thehojos.com/~mqh1/nearby_predictions.kml</a><br><br>' +\
'\nAll Predictions: <a href=http://thehojos.com/~mqh1/all_predictions.kml>http://thehojos.com/~mqh1/all_predictions.kml</a>'
sender.sendmail(sendTo, emailSubject, emailContent)
print("Email Sent")
No comments:
Post a Comment