Sunday, December 29, 2019

Weather sonde automated prediction

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. 

  1. Determining which predictions are coming in "near' Ithaac
  2. Generating two separate KML files.  One for "nearby" and one for "all" predictions.
  3. 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.

# 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