2024 Solar Eclipse Weather Report

#!/usr/bin/python3
#eclipse 2024 wether predictior - alnwlsn 2024
#using NWS forecast data api
#eclipse path extracted manually from KML at http://xjubier.free.fr/en/site_pages/SolarEclipsesGoogleEarth.html
import simplekml
import requests
import json
import csv
import time
import os.path
from datetime import datetime
import random
import subprocess

times = [
    "2024-04-07T16:00:00-04:00",
    "2024-04-08T13:00:00-04:00",
    "2024-04-08T14:00:00-04:00",
    "2024-04-08T15:00:00-04:00",
    "2024-04-08T16:00:00-04:00"
]

timesInterest = {}
for i in times:
    timesInterest[(datetime.fromisoformat(i).timestamp())] = i

def jj(j): #(debugging) write a json file to look at
    with open("jj.json", "w") as f:
        f.write(json.dumps(j, indent=2, sort_keys=True))

def get_grid_coordinates(x, y):
    url = f"https://api.weather.gov/points/{x},{y}"
    # print(url)
    response = requests.get(url)
    
    if response.status_code == 200:
        data = response.json()
        gridX = data['properties']['gridX']
        gridY = data['properties']['gridY']
        gridId = data['properties']['gridId']
        return gridId,gridX, gridY
    else:
        print("Error:", response.status_code)
        quit()
        return None, None
    
def get_forecast(grid_tuple): #gets hour-by hour forecast
    file = f'responses/{grid_tuple[0]},{grid_tuple[1]},{grid_tuple[2]}.json'
    if os.path.isfile(file): #do not get existing files
        return True
    url = f"https://api.weather.gov/gridpoints/{grid_tuple[0]}/{grid_tuple[1]},{grid_tuple[2]}/forecast/hourly"
    print(url)
    response = requests.get(url)
    
    if response.status_code == 200:
        data = response.json()
        with open(file, "w") as f:
            f.write(json.dumps(data, indent=2, sort_keys=True))
        return data
    else:
        print("Error:", response.status_code)
        if response.status_code == 500: #rate limiting timeout
            time.sleep(2)
        return False


def interpolate_points(points):
    interpolated_points = []
    for i in range(len(points) - 1):
        x1, y1 = points[i]
        x2, y2 = points[i + 1]
        interpolated_points.append((x1, y1))
        for j in range(1, 8):
            x_interpolated = x1 + (x2 - x1) * (j / 8)
            y_interpolated = y1 + (y2 - y1) * (j / 8)
            interpolated_points.append((x_interpolated, y_interpolated))
    interpolated_points.append(points[-1])  # Adding the last point
    return interpolated_points

def step1(): #get NWS grid coodinates of all points along the path
    path = []
    with open('path.csv') as f:
        index = 0
        rows = csv.reader(f)
        for row in rows:
            path.append((float(row[0]),float(row[1])))
            index += 1

    path = path[622:685] #trim this as needed
    path = interpolate_points(path)

    grid = []
    n = 0
    for i in path:
        g = get_grid_coordinates(i[1],i[0])
        # print(g)
        grid.append(g)
        # if n >= 4:
        #     break
        print(n,g)
        n+=1
        # time.sleep(1)
    with open("grid1.json", "w") as f:
        f.write(json.dumps(grid, indent=2, sort_keys=True))

def step2(): #extend the grid by some in up down grid direction (randomly)
    ex_grid = set()
    with open('grid1.json') as f:
        j = json.load(f)
        for i in j:
            # ex_grid.add((i[0], i[1]-4, i[2]))
            # ex_grid.add((i[0], i[1]-3, i[2]))
            # ex_grid.add((i[0], i[1]-2, i[2]))
            # ex_grid.add((i[0], i[1]-1, i[2]))
            # ex_grid.add((i[0], i[1], i[2]))
            # ex_grid.add((i[0], i[1]+1, i[2]))
            # ex_grid.add((i[0], i[1]+2, i[2]))
            # ex_grid.add((i[0], i[1]+3, i[2]))
            # ex_grid.add((i[0], i[1]+4, i[2]))
            for l in range(1,4):
                ex_grid.add((i[0], i[1]+random.randint(-30, 30), i[2]))

    ex_grid = list(ex_grid)
    # jj(ex_grid)
    print(ex_grid)
    print(len(ex_grid))
    with open("grid2.json", "w") as f:
        f.write(json.dumps(ex_grid, indent=2, sort_keys=True))

def step3(): #get the forecast for all the grid entries (store as a bunch of files)
    with open('grid2.json') as f:
    # with open('grid3.json') as f: #or this if you did step4 earlier
        j = json.load(f)
        # print(len(j))
        # time.sleep(4)
        for i in j:
            print(i)
            if False==get_forecast(i): #with 3 retries
                if False==get_forecast(i):
                    get_forecast(i)

def step4(): #reduce grid to ones I know are accessable by api (because I got the file in step3)
     g = []
     with open('grid2.json') as f:
        j = json.load(f)
        for i in j:
            file = f'responses/{i[0]},{i[1]},{i[2]}.json'
            if os.path.isfile(file):
                g.append(i)
        with open("grid3.json", "w") as f:
            f.write(json.dumps(g, indent=2, sort_keys=True))

def colorer(weather):
    #in BGR format
    opacity = "60"
    if "Mostly Cloudy" in weather:
        return opacity+'0080ff'
    if "Partly Cloudy" in weather:
        return opacity+'00bfff'
    if "Partly Sunny" in weather:
        return opacity+'00ffff'
    if "Mostly Sunny" in weather:
        return opacity+'00ffbf'
    if "Sunny" in weather:
        return opacity+'00ff00'
    return opacity+'0000ff' #unknown weather

def step5():
    kml = simplekml.Kml()
    subfolders = {}
    for i in timesInterest.keys():
        subfolders[i] = kml.newfolder(name=timesInterest[i])

    with open('grid2.json') as f:
        j = json.load(f)
        for i in j:
            file = f'responses/{i[0]},{i[1]},{i[2]}.json'
            if os.path.isfile(file):
                with open(file) as g:
                    color = '99ffac59'
                    forecast = json.load(g)
                    #reformat cordinates list into [(lon, lat),(lon, lat)...]
                    geo_raw = forecast['geometry']['coordinates'][0]
                    geo = []
                    for point in geo_raw:
                        geo.append((point[0],point[1]))
                    # quit()
                    hourly = forecast['properties']['periods']
                    for h in hourly:
                        startt = datetime.fromisoformat(h['startTime']).timestamp()
                        if startt in timesInterest.keys():
                            conditions = h['shortForecast']
                            color = colorer(conditions)
                            # print(timesInterest[startt],color,conditions)                       
                            pol = subfolders[startt].newpolygon(name=f'{i[0]},{i[1]},{i[2]},{conditions}', outerboundaryis=geo)
                            pol.style.linestyle.color = color
                            pol.style.polystyle.color = simplekml.Color.changealphaint(100, color)
                    # break
        # index += 1
        fstamp = 'eclipse-weather-'+str(int(time.time()))
        kml.save('/home/alnwlsn/temp/'+fstamp+".kml")

# step1() #first find the NWS grid coordinates for all points along the eclipse
# step2() #extend the collection with grid coords +1/-1 from the ones we have already
# step3() #grab forecast data for all grid points we generated. Try running this one a couple times
# step4() #store a list of all areas we know we can get, so I don't have to check against the API guess later (it can't do marine areas, for one)
step5() #make the KML