weather.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Get weather data\n",
    "\n",
    "Weather can impact the number of bike riders on the road. We want to get weather data from 2013-2022 from [NOAA](https://www.ncei.noaa.gov/data/global-hourly/archive/csv/)\n",
    "* Bulk download the dataset for specified years from https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.ncdc:C00516\n",
    "* Untar and extract\n",
    "* Keep only weather station data for the KNYC (Central Park, Manhattan) weather station\n",
    "\n",
    "Weather stations and can be found [here](https://www.ncei.noaa.gov/access/search/data-search/global-summary-of-the-day?startDate=2013-01-01T00:00:00&endDate=2022-03-20T23:59:59&pageNum=1&stations=72505394728). Data can by downloaded manually for specific stations and years in lieu of running this script\n",
    "* NY CITY CENTRAL PARK, NY US (72505394728.csv)\n",
    "\n",
    "### Data\n",
    "We are interested in temperature, precipitation, and whether a day had snow\n",
    "* Global Hourly (downloaded using this script)\n",
    "* Global Summary of Day (manually downloaded from [here](https://www.ncei.noaa.gov/access/search/data-search/global-summary-of-the-day?startDate=2013-01-01T00:00:00&endDate=2022-03-20T23:59:59&pageNum=1&stations=72505394728))\n",
    "* Daily Summary (This is newer and replaces GSOD above, so actually using this for analysis)\n",
    "\n",
    "### Data dictionaries\n",
    "* [Global Hourly CSV Help](https://www.ncei.noaa.gov/data/global-hourly/doc/CSV_HELP.pdf)\n",
    "* [Global Hourly ISD Format](https://www.ncei.noaa.gov/data/global-hourly/doc/isd-format-document.pdf)\n",
    "* [Global Summary of Day](https://www.ncdc.noaa.gov/cdo-web/datasets#global-summary-of-the-day)\n",
    "* [Daily Summary](https://www.ncdc.noaa.gov/cdo-web/datasets#GHCND)\n",
    "\n",
    "# WARNING: Running the notebook will download many gigabytes of data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import requests\n",
    "import logging\n",
    "import sys\n",
    "import tarfile\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "YEARS = range(2013, 2023)\n",
    "NOAA_URL = \"https://www.ncei.noaa.gov/data/global-hourly/archive/csv/\"\n",
    "FILE_EXTENSION = \".tar.gz\"\n",
    "WEATHER_DIR = \"data/weather/\"\n",
    "HOURLY_DIR = WEATHER_DIR + \"hourly/\"\n",
    "DAILY_DIR = WEATHER_DIR + \"daily/\"\n",
    "KNYC_FILE = \"72505394728.csv\"\n",
    "GHCN_FILE = \"GHCN-Daily.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# create directory to save weather data to\n",
    "if not os.path.exists(WEATHER_DIR):\n",
    "    os.makedirs(os.path.dirname(WEATHER_DIR))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Global Hourly Weather Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# download all weather *tar.gz files\n",
    "for year in YEARS:\n",
    "    noaa_url = NOAA_URL + str(year) + FILE_EXTENSION\n",
    "    with requests.get(noaa_url, stream=True) as r:\n",
    "        r.raise_for_status()\n",
    "\n",
    "        local_filename = WEATHER_DIR + str(year) + FILE_EXTENSION\n",
    "        logging.debug(f\"downloaded {noaa_url}. saving to {local_filename}...\")\n",
    "        with open(local_filename, \"w+b\") as f:\n",
    "            f.write(r.raw.read(decode_content=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# extract (only KNYC weather data)\n",
    "for year in YEARS:\n",
    "    tar = tarfile.open(WEATHER_DIR + str(year) + FILE_EXTENSION)\n",
    "    YEAR_DIR = WEATHER_DIR + str(year) + \"/\"\n",
    "\n",
    "    if not os.path.exists(YEAR_DIR):\n",
    "        os.makedirs(os.path.dirname(YEAR_DIR))\n",
    "\n",
    "    tar.extract(member=KNYC_FILE, path=YEAR_DIR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def c2f(celsius):\n",
    "    return celsius * 9 / 5 + 32"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# create hourly weather dir\n",
    "if not os.path.exists(HOURLY_DIR):\n",
    "    os.makedirs(os.path.dirname(HOURLY_DIR))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# process and save hourly weather data\n",
    "# \"hourly\" weather for NYC from KNYC weather station\n",
    "for year in YEARS:\n",
    "    hourly_file = WEATHER_DIR + str(year) + \"/\" + KNYC_FILE\n",
    "    hourly = pd.read_csv(hourly_file, usecols=[\"DATE\", \"TMP\", \"AA1\"])\n",
    "    hourly.dropna(axis=0, inplace=True)\n",
    "\n",
    "    # TMP -> keep digits before comma -> divide by 10, keep sign (in celcius -> convert to F)\n",
    "    hourly[\"TMP\"] = (\n",
    "        hourly[\"TMP\"]\n",
    "        .apply(lambda s: int(s.split(\",\")[0]) / 10)\n",
    "        .replace(999.9, np.NAN)\n",
    "        .map(c2f)\n",
    "    )\n",
    "    # AA1 -> precipitation. <?, amount precip in mm/10, ?, ?>\n",
    "    hourly[\"AA1\"] = (\n",
    "        hourly[\"AA1\"].apply(lambda a: int(a.split(\",\")[1]) / 10).replace(999.9, np.NAN)\n",
    "    )\n",
    "    # DATE -> need to round to nearest hour, drop duplicate hours, if no hour impute\n",
    "\n",
    "    # drop the 999.9->NAs as time series model needs data\n",
    "    hourly.dropna(axis=0, inplace=True)\n",
    "    hourly.rename(\n",
    "        columns={\"DATE\": \"date\", \"TMP\": \"temp_f\", \"AA1\": \"precip_mm\"}, inplace=True\n",
    "    )\n",
    "\n",
    "    hourly.to_csv(HOURLY_DIR + str(year) + \".csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Global Summary of Daily Weather Data\n",
    "NOTE: This data was downloaded manually as specified above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# create daily weather dir\n",
    "if not os.path.exists(DAILY_DIR):\n",
    "    os.makedirs(os.path.dirname(DAILY_DIR))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# process and save daily weather data\n",
    "# daily weather for NYC from KNYC weather station\n",
    "for year in YEARS:\n",
    "    daily_file = WEATHER_DIR + str(year) + \".csv\"\n",
    "    # daily weather from KNYC weather station: NY CITY CENTRAL PARK, NY US\n",
    "    daily = pd.read_csv(daily_file, usecols=[\"DATE\", \"TEMP\", \"SNDP\", \"PRCP\"])\n",
    "\n",
    "    # TEMP -> leave as is (in Farenheit)\n",
    "    # SNDP in cm -> 0 for 999.9 and 1 for any snow\n",
    "    daily[\"SNDP\"].replace(999.9, 0, inplace=True)\n",
    "    daily[\"SNDP\"] = np.sign(daily[\"SNDP\"])\n",
    "    daily[\"SNDP\"] = daily[\"SNDP\"].astype(\"int\")  # 0/1 instead of boolean?\n",
    "\n",
    "    daily.rename(\n",
    "        columns={\n",
    "            \"DATE\": \"date\",\n",
    "            \"TEMP\": \"temp_f\",\n",
    "            \"SNDP\": \"is_snowday\",\n",
    "            \"PRCP\": \"precip_mm\",\n",
    "        },\n",
    "        inplace=True,\n",
    "    )\n",
    "\n",
    "    daily.to_csv(DAILY_DIR + str(year) + \".csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GHCN-Daily (Daily Summary)\n",
    "https://www.ncdc.noaa.gov/cdo-web/datasets#GHCND\n",
    "\n",
    "This dataset has more information than GSOD. The file `GHCN-Daily.csv` is for data from KNYC from Jan 1 2013 to March 20 2022"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "ghcn = pd.read_csv(\n",
    "    WEATHER_DIR + GHCN_FILE,\n",
    "    usecols=[\n",
    "        \"DATE\",\n",
    "        \"PRCP\",\n",
    "        \"SNOW\",\n",
    "        \"SNWD\",\n",
    "        \"TMAX\",\n",
    "        \"TMIN\",\n",
    "        \"WT01\",\n",
    "        \"WT02\",\n",
    "        \"WT22\",\n",
    "        \"WT03\",\n",
    "        \"WT08\",\n",
    "        \"WT16\",\n",
    "    ],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# See: https://www1.ncdc.noaa.gov/pub/data/cdo/documentation/GHCND_documentation.pdf\n",
    "# units are metric (mm, celsius)\n",
    "# PRCP (precipitation mm); SNOW (snowfall mm); SNWD (snow depth mm); TMAX and TMIN (celsius)\n",
    "# Weather attributes:\n",
    "# WT01 -> fog\n",
    "# WT02 -> heavy fog\n",
    "# WT22 -> ice fog\n",
    "# WT03 -> thunder\n",
    "# WT08 -> haze, smoke\n",
    "# WT16 -> rain\n",
    "ghcn[\"TAVG_F\"] = round(c2f((ghcn[\"TMAX\"] + ghcn[\"TMIN\"]) / 2), 1)\n",
    "ghcn[\"Normal Fog\"] = np.invert(np.isnan(ghcn[\"WT01\"]))\n",
    "ghcn[\"Heavy Fog\"] = np.invert(np.isnan(ghcn[\"WT02\"]))\n",
    "ghcn[\"Ice Fog\"] = np.invert(np.isnan(ghcn[\"WT22\"]))\n",
    "ghcn[\"Fog\"] = ghcn[\"Normal Fog\"] | ghcn[\"Heavy Fog\"] | ghcn[\"Ice Fog\"]\n",
    "ghcn[\"Thunder\"] = np.invert(np.isnan(ghcn[\"WT03\"]))\n",
    "ghcn[\"Haze\"] = np.invert(np.isnan(ghcn[\"WT08\"]))\n",
    "ghcn[\"Rain\"] = np.invert(np.isnan(ghcn[\"WT16\"]))\n",
    "\n",
    "# convert boolean to int\n",
    "ghcn[[\"Normal Fog\", \"Heavy Fog\", \"Ice Fog\", \"Fog\", \"Thunder\", \"Haze\", \"Rain\"]] = ghcn[\n",
    "    [\"Normal Fog\", \"Heavy Fog\", \"Ice Fog\", \"Fog\", \"Thunder\", \"Haze\", \"Rain\"]\n",
    "].replace({True: 1, False: 0})\n",
    "# drop converted columns\n",
    "ghcn.drop([\"TMAX\", \"TMIN\"] + list(ghcn.filter(regex=\"WT\")), axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "ghcn.to_csv(WEATHER_DIR + \"GHCN-Daily-Cleaned.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}