{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Station Plot\n",
    "\n",
    "This notebook is adapted from https://unidata.github.io/MetPy/latest/examples/plots/Station_Plot.html\n",
    "\n",
    "Make a station plot, complete with sky cover and weather symbols.\n",
    "\n",
    "The station plot itself is pretty straightforward, but there is a bit of code to perform the\n",
    "data-wrangling (hopefully that situation will improve in the future). Certainly, if you have\n",
    "existing point data in a format you can work with trivially, the station plot will be simple.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from metpy.calc import reduce_point_density\n",
    "from metpy.cbook import get_test_data\n",
    "from metpy.io import metar\n",
    "from metpy.plots import add_metpy_logo, current_weather, sky_cover, StationPlot\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get real-time metar data\n",
    "# Author: Yongjie Huang\n",
    "# Date: 2022-02-24\n",
    "from os.path import exists\n",
    "from urllib.request import urlretrieve\n",
    "\n",
    "def get_real_time_metar_data(fname):\n",
    "    base_url = 'https://thredds-test.unidata.ucar.edu/thredds/fileServer/noaaport/text/metar/'\n",
    "    data_dir = '/tmp/'\n",
    "    path_to_file = data_dir + fname\n",
    "    \n",
    "    if exists(path_to_file):\n",
    "        print(f'File exists: {path_to_file}')\n",
    "    else:\n",
    "        print(f\"Downloading file '{fname}' from '{base_url+fname}' to '{path_to_file}'.\")\n",
    "        urlretrieve(base_url+fname, path_to_file)\n",
    "\n",
    "    return path_to_file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The setup\n",
    "\n",
    "First download and read in the station data. \n",
    "\n",
    "We use the get_real_time_metar_data function defined above to download 'realtime' data from Unidata THREDDS server.\n",
    "the metar reader because it simplifies a lot of tasks,\n",
    "like dealing with separating text and assembling a pandas dataframe.\n",
    "\n",
    "https://thredds-test.unidata.ucar.edu/thredds/catalog/noaaport/text/metar/catalog.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Download metar data in text format.\n",
    "\n",
    "date =  datetime(2022, 3, 21, 12, 0)\n",
    "print( date )\n",
    "\n",
    "datestring = datetime.strftime( date, '%Y%m%d_%H%M' )\n",
    "print( datestring )\n",
    "\n",
    "filename = 'metar_'+datestring+'.txt'\n",
    "\n",
    "metar_data = get_real_time_metar_data(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "Parse a downloaded text file containing multiple METAR reports and/or text products, by calling metar.parse_metar_file which is part of the metpy.io package. The parser returns the data as pandas dataframe. See\n",
    "https://unidata.github.io/MetPy/latest/api/generated/metpy.io.parse_metar_file.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data = metar.parse_metar_file(metar_data)\n",
    "\n",
    "#data = metar.parse_metar_file(get_test_data('metar_20190701_1200.txt', as_file_obj=False))\n",
    "\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "print( data['station_id'].shape )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Drop rows with missing winds and sea level pressure\n",
    "data = data.dropna(how='any', subset=['wind_direction', 'wind_speed','air_pressure_at_sea_level','air_temperature'])\n",
    "print( data['station_id'].shape )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Set up the map projection\n",
    "\n",
    "\n",
    "proj = ccrs.LambertConformal(central_longitude=-95, central_latitude=35,\n",
    "                               standard_parallels=[35, 40])\n",
    "\n",
    "# Use the Cartopy map projection to transform station locations in lat/lon to the map coordinates\n",
    "point_locs = proj.transform_points(ccrs.PlateCarree(), data['longitude'].values,\n",
    "                                   data['latitude'].values)\n",
    "\n",
    "point_locs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "print( point_locs.shape, data['station_id'].values.shape )\n",
    "\n",
    "nstations = data['station_id'].values.size\n",
    "print( nstations )\n",
    "\n",
    "nstations\n",
    "data['longitude'][1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Screen observations and keep only those within the lat/lon bounds. Put kept observations and their (x,y) and (lat,lon) locations in new arrays. Rememeber how many observations are kept (n_stations)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xstation = np.zeros( nstations )\n",
    "ystation = np.zeros( nstations )\n",
    "lat_keep = np.zeros( nstations )\n",
    "lon_keep = np.zeros( nstations )\n",
    "pmsl_keep = np.zeros( nstations )\n",
    "t_keep = np.zeros( nstations )\n",
    "\n",
    "lonmin = -140\n",
    "lonmax = -50\n",
    "latmin = 20\n",
    "latmax = 60\n",
    "kk = 0\n",
    "print( nstations )\n",
    "\n",
    "for k in range(0, nstations ):\n",
    "    lonk = data['longitude'][k] \n",
    "    latk = data['latitude'][k] \n",
    "    \n",
    "    if ( (lonk-lonmax)*(lonk-lonmin) <= 0.0 and (latk-latmax)*(latk-latmin) <= 0.0):\n",
    "        xstation[kk] = point_locs[k,0]\n",
    "        ystation[kk] = point_locs[k,1]\n",
    "        lat_keep[kk] = latk\n",
    "        lon_keep[kk] = lonk\n",
    "        pmsl_keep[kk] = data['air_pressure_at_sea_level'].values[k]\n",
    "        t_keep[kk] = data['air_temperature'].values[k]\n",
    "        if np.isnan( pmsl_keep[kk]) or np.isnan( t_keep[kk]):  # check if pmsl is NaN/nan. If so, skip that station (don't count into kk)\n",
    "            print( k, pmsl_keep[kk] )\n",
    "        else:\n",
    "            kk += 1\n",
    "nstations_kept = kk + 1  # The total number of station kept should be the final count kk plus 1 (kk starts from 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "print(nstations_kept, np.min(xstation[0:nstations_kept]), np.max(xstation[0:nstations_kept]), \n",
    "      np.min(ystation[0:nstations_kept]), np.max(ystation[0:nstations_kept]) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Define  (x, y) coordinates of the objective analysis grid in the map projection space based on the \n",
    "# max/min x and y values found above.\n",
    "#\n",
    "xgrid = np.array( range( -2700000, 3100000, 50000 ))\n",
    "ygrid = np.array( range( -1600000, 3000000, 50000 ))\n",
    "nxgrid = xgrid.size\n",
    "nygrid = ygrid.size\n",
    "\n",
    "print( nxgrid, nygrid )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Calculate average station distance. \n",
    "\n",
    "Suggestions: Define a numpy array of size nstations_kept to store the distance to the nearest station for each station. Fill this array with distance beyond the largest possible, like 10,000,000 meter. For each station, find the distance of the closest station that is not itself. It turns out that there are many duplicate stations in the data so that the closest station distance for some stations is zero. So before averaging is done to find the average distance, make sure you skip the zero distances, either by using np.where function to find the non-zero distance indices to use for np.mean average fund (see Dr. Furtado's Notebook example), or write loop code yourself to average only the non-zero distances.  The average distance you find should tens of km. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "nearest_station_distance = np.zeros( nstations_kept )\n",
    "nearest_station_distance.fill(10000000)\n",
    "\n",
    "for k in range(0, nstations_kept ):\n",
    "    for kk in range(0, nstations_kept ):\n",
    "\n",
    "# write your codes here.\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# find the average station distance\n",
    "\n",
    "average_r = \n",
    "max_r = \n",
    "min_r = \n",
    "\n",
    "print( \"Average station distance =\", average_r )\n",
    "print( \"Max station distance =\", max_r )\n",
    "print( \"Min station distance =\", min_r )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Perform Cressman objective analysis for mean sea level pressure and surface temperature and store results in zgrid and tgrid respectively.\n",
    "\n",
    "## Please use 2x and 5x the average station spacing, and run the analysis twice and submit both results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "R2 = (average_r*5)**2\n",
    "\n",
    "zgrid = np.zeros((nygrid, nxgrid)) \n",
    "tgrid = np.zeros((nygrid, nxgrid)) \n",
    "\n",
    "mean_station_pmsl = \n",
    "mean_station_t = \n",
    "\n",
    "print( 'mean station psml and temperature:', mean_station_pmsl,  mean_station_t )\n",
    "\n",
    "for i in range( 0, nxgrid ):\n",
    "    for j in range(0, nygrid ): \n",
    "       \n",
    "# write your code here    \n",
    "    \n",
    "print('min and max of analyzed psml: ', np.min(zgrid), np.max(zgrid))\n",
    "print('min and max of analyzed temperature: ', np.min(tgrid), np.max(tgrid))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot analyzed pressure field together with station markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Change the DPI of the resulting figure. Higher DPI drastically improves the\n",
    "# look of the text rendering.\n",
    "plt.rcParams['savefig.dpi'] = 600\n",
    "\n",
    "# Create the figure and an axes set to the projection.\n",
    "fig = plt.figure(figsize=(20, 10))\n",
    "add_metpy_logo(fig, 1100, 300, size='large')\n",
    "ax = fig.add_subplot(1, 1, 1, projection=proj)\n",
    "\n",
    "# Add some various map elements to the plot to make it recognizable.\n",
    "ax.add_feature(cfeature.LAND)\n",
    "ax.add_feature(cfeature.OCEAN)\n",
    "ax.add_feature(cfeature.LAKES)\n",
    "ax.add_feature(cfeature.COASTLINE)\n",
    "ax.add_feature(cfeature.STATES)\n",
    "ax.add_feature(cfeature.BORDERS, linewidth=2)\n",
    "\n",
    "# Set plot bounds\n",
    "\n",
    "ax.set_extent((-125, -65, 20, 52))\n",
    "#\n",
    "# plot objectively analyzed mean sea level pressure field.\n",
    "#\n",
    "x2d, y2d = np.meshgrid( xgrid, ygrid )\n",
    "ctr = ax.contour(x2d, y2d, zgrid, colors='red', linewidths=2,transform=proj, levels=range(994,1040,2))\n",
    "ctr = ax.contourf(x2d, y2d, zgrid, transform=proj, levels=range(994,1040,2))\n",
    "ax.clabel(ctr, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)\n",
    "\n",
    "\n",
    "# Start the station plot by specifying the axes to draw on, as well as the\n",
    "# lon/lat of the stations (with transform). We also the fontsize to 12 pt.\n",
    "stationplot = StationPlot(ax, data['longitude'].values, data['latitude'].values,\n",
    "                          clip_on=True, transform=ccrs.PlateCarree(), fontsize=6)\n",
    "\n",
    "# Plot the cloud cover symbols in the center location. This uses the codes made above and\n",
    "# uses the `sky_cover` mapper to convert these values to font codes for the\n",
    "# weather symbol font.\n",
    "stationplot.plot_symbol('C', data['cloud_coverage'].values, sky_cover)\n",
    "\n",
    "\n",
    "Title = f\"Surface stations and mean sea level pressure contours for {date}\"\n",
    "ax.set_title(Title.title(), fontsize=16)\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Discusss the results obtained above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Plot analyzed surface temperature field together with station markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Change the DPI of the resulting figure. Higher DPI drastically improves the\n",
    "# look of the text rendering.\n",
    "plt.rcParams['savefig.dpi'] = 600\n",
    "\n",
    "# Create the figure and an axes set to the projection.\n",
    "fig = plt.figure(figsize=(20, 10))\n",
    "add_metpy_logo(fig, 1100, 300, size='large')\n",
    "ax = fig.add_subplot(1, 1, 1, projection=proj)\n",
    "\n",
    "# Add some various map elements to the plot to make it recognizable.\n",
    "ax.add_feature(cfeature.LAND)\n",
    "ax.add_feature(cfeature.OCEAN)\n",
    "ax.add_feature(cfeature.LAKES)\n",
    "ax.add_feature(cfeature.COASTLINE)\n",
    "ax.add_feature(cfeature.STATES)\n",
    "ax.add_feature(cfeature.BORDERS, linewidth=2)\n",
    "\n",
    "# Set plot bounds\n",
    "\n",
    "ax.set_extent((-125, -65, 20, 52))\n",
    "#\n",
    "# plot objectively analyzed mean sea level pressure field.\n",
    "#\n",
    "x2d, y2d = np.meshgrid( xgrid, ygrid )\n",
    "ctr = ax.contour(x2d, y2d, tgrid, colors='red', linewidths=2,transform=proj, levels=range(-20,30,1))\n",
    "ctr = ax.contourf(x2d, y2d, tgrid, transform=proj, levels=range(-20,30,1))\n",
    "ax.clabel(ctr, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)\n",
    "\n",
    "\n",
    "# Start the station plot by specifying the axes to draw on, as well as the\n",
    "# lon/lat of the stations (with transform). We also the fontsize to 12 pt.\n",
    "stationplot = StationPlot(ax, data['longitude'].values, data['latitude'].values,\n",
    "                          clip_on=True, transform=ccrs.PlateCarree(), fontsize=6)\n",
    "\n",
    "# Plot the cloud cover symbols in the center location. This uses the codes made above and\n",
    "# uses the `sky_cover` mapper to convert these values to font codes for the\n",
    "# weather symbol font.\n",
    "stationplot.plot_symbol('C', data['cloud_coverage'].values, sky_cover)\n",
    "\n",
    "\n",
    "Title = f\"Surface stations  and surface temperature (C) for {date}\"\n",
    "ax.set_title(Title.title(), fontsize=16)\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Discusss the results obtained above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}