tornado-map/tornados_us.py

# -*- coding: utf-8 -*-

import csv
from dataclasses import dataclass
from datetime import datetime

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from cartopy import geodesic
import shapely.geometry as sgeom


proj = ccrs.AlbersEqualArea(
    central_longitude=-95,
    central_latitude=40,
)
water_blue = "#7ebfd4"


def phys(name, resolv):
    return cfeature.NaturalEarthFeature('physical', name, resolv, facecolor="none")


def cultural(name, resolv):
    return cfeature.NaturalEarthFeature('cultural', name, resolv, facecolor="none")


land = phys("land", "50m")
rivers = phys("rivers_lake_centerlines", "50m")
lakes = phys("lakes", "50m")
countries = cultural("admin_0_countries", "50m")
states = cultural("admin_1_states_provinces_lines", "50m")


@dataclass
class Tornado:
    om: str
    yr: str
    mo: str
    dy: str
    date: str
    time: str
    tz: str
    st: str
    stf: str
    stn: str
    mag: str
    inj: str
    fat: str
    loss: str
    closs: str
    slat: str
    slon: str
    elat: str
    elon: str
    len: str
    wid: str
    ns: str
    sn: str
    sg: str
    f1: str
    f2: str
    f3: str
    f4: str
    fc: str

    def __post_init__(self):
        self.id = int(self.om)
        self.dt = datetime.strptime(self.date, "%Y-%m-%d")
        self.slon = float(self.slon)
        self.slat = float(self.slat)
        self.elon = float(self.elon)
        self.elat = float(self.elat)

    def __repr__(self):
        return "<tornado {id} {date}>".format(id=self.id, date=self.dt)


def read_data():
    tornados = []
    with open("1950-2017_actual_tornadoes.csv", newline='') as csvfile:
        reader = csv.DictReader(csvfile)
        for i, row in enumerate(reader):
            tornados.append(Tornado(**row))
    return tornados


def draw_map(data):
    fig = plt.figure(frameon=False)
    ax = fig.add_subplot(111, projection=proj)
    ax.background_patch.set_facecolor(water_blue)
    ax.set_extent([-125, -65, 20, 50], ccrs.Geodetic())
    ax.add_feature(land, facecolor="#f0f0f0")
    ax.add_feature(rivers, edgecolor=water_blue)
    ax.add_feature(lakes, facecolor=water_blue)

    ax.add_feature(countries, edgecolor="grey", linewidth=0.2, alpha=0.6, dashes='--')
    ax.add_feature(states, edgecolor="grey", linewidth=0.2, alpha=0.4, dashes='--')

    for i, t in enumerate(data):
        print(i)
        if t.elon < 0:
            ax.plot([t.slon, t.elon], [t.slat, t.elat], 'r', lw=0.3, transform=ccrs.Geodetic())
        else:
            ax.plot([t.slon, t.slon], [t.slat, t.slat], 'r', lw=0.3, transform=ccrs.Geodetic())

    fig.tight_layout()
    plt.savefig("tornados_us.png", dpi=320)


if __name__ == '__main__':
    ts = read_data()
    draw_map(ts)