from datetime import datetime, timedelta

import xarray as xr

import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.io import metar
from metpy.plots.declarative import (BarbPlot, ContourPlot, FilledContourPlot, MapPanel,
                                     PanelContainer, PlotObs)
from metpy.units import units

"""
Explanation: MetPy Declarative Syntax Tutorial
The declarative syntax that is a part of the MetPy packaged is designed to aid in simple
data exploration and analysis needs by simplifying the plotting context from typical verbose
Python code. The complexity of data wrangling and plotting are hidden behind the simplified
syntax to allow a lower barrier to investigating your data.
Imports
You'll note that the number of imports is smaller due to using the declarative syntax.
There is no need to import Matplotlib or Cartopy to your code as all of that is done
behind the scenes.
End of explanation
"""


# Open the netCDF file as a xarray Dataset and parse the full dataset
data = xr.open_dataset(get_test_data('GFS_test.nc', False)).metpy.parse_cf()

# View a summary of the Dataset
print(data)

"""
Explanation: Getting Data
Depending on what kind of data you are wanting to plot you'll use either Xarray (for gridded
data), Pandas (for CSV data), or the MetPy METAR parser (for METAR data).
We'll start this tutorial by reading in a gridded dataset using Xarray.
End of explanation
"""


plot_time = datetime(2010, 10, 26, 12)

"""
Explanation: Set Datetime
Set the date/time of that you desire to plot
End of explanation
"""


ds = data.metpy.sel(lat=slice(70, 10), lon=slice(360 - 150, 360 - 55))

"""
Explanation: Subsetting Data
MetPy provides wrappers for the usual xarray indexing and selection routines that can handle
quantities with units. For DataArrays, MetPy also allows using the coordinate axis types
mentioned above as aliases for the coordinates. And so, if we wanted data to be just over
the U.S. for plotting purposes
End of explanation
"""


ds['wind_speed'] = mpcalc.wind_speed(ds['u-component_of_wind_isobaric'],
                                     ds['v-component_of_wind_isobaric'])

"""
Explanation: For full details on xarray indexing/selection, see
xarray's documentation &lt;http://xarray.pydata.org/en/stable/indexing.html&gt;_.
Calculations
In MetPy 1.0 and later, calculation functions accept Xarray DataArray's as input and the
output a DataArray that can be easily added to an existing Dataset.
As an example, we calculate wind speed from the wind components and add it as a new variable
to our Dataset.
End of explanation
"""


# Set attributes for contours of Geopotential Heights at 300 hPa
cntr2 = ContourPlot()
cntr2.data = ds
cntr2.field = 'Geopotential_height_isobaric'
cntr2.level = 300 * units.hPa
cntr2.time = plot_time
cntr2.contours = list(range(0, 10000, 120))
cntr2.linecolor = 'black'
cntr2.linestyle = 'solid'
cntr2.clabels = True

"""
Explanation: Plotting
With that miniaml preparation, we are now ready to use the simplified plotting syntax to be
able to plot our data and analyze the meteorological situation.
General Structure


Set contour attributes


Set map characteristics and collect contours


Collect panels and plot


Show (or save) the results


Valid Plotting Types for Gridded Data:


ContourPlot()


FilledContourPlot()


ImagePlot()


PlotBarbs()


More complete descriptions of these and other plotting types, as well as the map panel and
panel container classes are at the end of this tutorial.
Let's plot a 300-hPa map with color-filled wind speed, which we calculated and added to
our Dataset above, and geopotential heights over the CONUS.
We'll start by setting attributes for contours of Geopotential Heights at 300 hPa.
We need to set at least the data, field, level, and time attributes. We'll set a few others
to have greater control over hour the data is plotted.
End of explanation
"""


# Set attributes for plotting color-filled contours of wind speed at 300 hPa
cfill = FilledContourPlot()
cfill.data = ds
cfill.field = 'wind_speed'
cfill.level = 300 * units.hPa
cfill.time = plot_time
cfill.contours = list(range(10, 201, 20))
cfill.colormap = 'BuPu'
cfill.colorbar = 'horizontal'
cfill.plot_units = 'knot'

"""
Explanation: Now we'll set the attributes for plotting color-filled contours of wind speed at 300 hPa.
Again, the attributes that must be set include data, field, level, and time. We'll also set
a colormap and colorbar to be purposeful for wind speed. Additionally, we'll set the
attribute to change the units from m/s to knots, which is the common plotting units for
wind speed.
End of explanation
"""


# Set the attributes for the map and add our data to the map
panel = MapPanel()
panel.area = [-125, -74, 20, 55]
panel.projection = 'lcc'
panel.layers = ['states', 'coastline', 'borders']
panel.title = f'{cfill.level.m}-hPa Heights and Wind Speed at {plot_time}'
panel.plots = [cfill, cntr2]

"""
Explanation: Once we have our contours (and any colorfill plots) set up, we will want to define the map
panel that we'll plot the data on. This is the place where we can set the view extent,
projection of our plot, add map lines like coastlines and states, set a plot title.
One of the key elements is to add the data to the map panel as a list with the plots
attribute.
End of explanation
"""


# Set the attributes for the panel and put the panel in the figure
pc = PanelContainer()
pc.size = (15, 15)
pc.panels = [panel]

"""
Explanation: Finally we'll collect all of the panels to plot on the figure, set the size of the figure,
and ultimately show or save the figure.
End of explanation
"""


# Show the image
pc.show()

"""
Explanation: All of our setting now produce the following map!
End of explanation
"""


# Set attributes for contours of Geopotential Heights at 300 hPa
cntr2 = ContourPlot()
cntr2.data = ds
cntr2.field = 'Geopotential_height_isobaric'
cntr2.level = 300 * units.hPa
cntr2.time = plot_time
cntr2.contours = list(range(0, 10000, 120))
cntr2.linecolor = 'black'
cntr2.linestyle = 'solid'
cntr2.clabels = True

# Set attributes for plotting color-filled contours of wind speed at 300 hPa
cfill = FilledContourPlot()
cfill.data = ds
cfill.field = 'wind_speed'
cfill.level = 300 * units.hPa
cfill.time = plot_time
cfill.contours = list(range(10, 201, 20))
cfill.colormap = 'BuPu'
cfill.colorbar = 'horizontal'
cfill.plot_units = 'knot'

"""
Explanation: That's it! What a nice looking map, with relatively simple set of code.
Adding Wind Barbs
We can easily add wind barbs to the plot we generated above by adding another plot type
and adding it to the panel. The plot type for wind barbs is PlotBarbs() and has its own
set of attributes to control plotting a vector quantity.
We start with setting the attributes that we had before for our 300 hPa plot including,
Geopotential Height contours, and color-filled wind speed.
End of explanation
"""


# Set attributes for plotting wind barbs
barbs = BarbPlot()
barbs.data = ds
barbs.time = plot_time
barbs.field = ['u-component_of_wind_isobaric', 'v-component_of_wind_isobaric']
barbs.level = 300 * units.hPa
barbs.skip = (3, 3)
barbs.plot_units = 'knot'

"""
Explanation: Now we'll set the attributes for plotting wind barbs, with the required attributes of data,
time, field, and level. The skip attribute is particularly useful for thining the number of
wind barbs that are plotted on the map and again we'll convert to units of knots.
End of explanation
"""


# Set the attributes for the map and add our data to the map
panel = MapPanel()
panel.area = [-125, -74, 20, 55]
panel.projection = 'lcc'
panel.layers = ['states', 'coastline', 'borders']
panel.title = f'{cfill.level.m}-hPa Heights and Wind Speed at {plot_time}'
panel.plots = [cfill, cntr2, barbs]

# Set the attributes for the panel and put the panel in the figure
pc = PanelContainer()
pc.size = (15, 15)
pc.panels = [panel]

# Show the figure
pc.show()

"""
Explanation: Add all of our plot types to the panel, don't forget to add in the new wind barbs to our plot
list!
End of explanation
"""


df = metar.parse_metar_file(get_test_data('metar_20190701_1200.txt', False), year=2019,
                            month=7)

# Let's take a look at the variables that we could plot coming from our METAR observations.
print(df.keys())

# Set the observation time
obs_time = datetime(2019, 7, 1, 12)

"""
Explanation: Plot Surface Obs
We can also plot surface (or upper-air) observations at point locations using the simplified
syntax. Whether it is surface or upper-air data, the PlotObs() class is what you would
want to use. Then you would add those observations to a map panel and then collect the panels
to plot the figure; similar to what you would do for a gridded plot.
End of explanation
"""


# Plot desired data
obs = PlotObs()
obs.data = df
obs.time = obs_time
obs.time_window = timedelta(minutes=15)
obs.level = None
obs.fields = ['cloud_coverage', 'air_temperature', 'dew_point_temperature',
              'air_pressure_at_sea_level', 'current_wx1_symbol']
obs.plot_units = [None, 'degF', 'degF', None, None]
obs.locations = ['C', 'NW', 'SW', 'NE', 'W']
obs.formats = ['sky_cover', None, None, lambda v: format(v * 10, '.0f')[-3:],
               'current_weather']
obs.reduce_points = 0.75
obs.vector_field = ['eastward_wind', 'northward_wind']

"""
Explanation: Setting of our attributes for plotting observations is pretty straignforward and just needs
to be lists for the variables, and a comparable number of items for plot characteristics that
are specific to the individual fields. For example, the locations around a station plot, the
plot units, and any plotting formats would all meed to have the same number of items as the
fields attribute.
Plotting wind bards is done through the vector_field attribute and you can reduce the number
of points plotted (especially important for surface observations) with the reduce points
attribute.
For a very basic plot of one field, the minimum required attributes are the data, time,
fields, and location attributes.
End of explanation
"""


# Panel for plot with Map features
panel = MapPanel()
panel.layout = (1, 1, 1)
panel.projection = 'lcc'
panel.area = 'in'
panel.layers = ['states']
panel.title = f'Surface plot for {obs_time}'
panel.plots = [obs]

# Bringing it all together
pc = PanelContainer()
pc.size = (10, 10)
pc.panels = [panel]

pc.show()

"""
Explanation: We use the same Classes for plotting our data on a map panel and collecting all of the
panels on the figure. In this case we'll focus in on the state of Indiana for plotting.
End of explanation
"""


