Map noaa average wind speed data given a cdf4 file

noaa_avg_windspd_map.py

from netCDF4 import Dataset
from datetime import datetime, timedelta
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap


class MonthMean(object):

    def __init__(self, data, month, year):
        self.data = data
        self.month = month
        self.year = year

    def __repr__(self):
        return self.month + ',' + str(self.year)


class NOAAData(object):

    def __init__(self, filep, variable=None):
        self._month_names = [
            'January', 'February', 'March', 'April', 'May', 'June', 'July',
            'August', 'September', 'October', 'November', 'December']
        dataset = Dataset(filep)
        date = dataset.variables['time'].units
        self.date = datetime.strptime(date[12:-2], '%Y-%m-%d %H:%M:%S')
        start, end = dataset.variables['time'].actual_range
        self.start = self.date + timedelta(hours=start)
        self.end = self.date + timedelta(hours=end)
        self.lon = dataset.variables['lon'][:]
        if self.lon.max() > 180:
            self.lon -= 180
        self.lat = dataset.variables['lat'][:]
        variable = variable if variable else dataset.variables.keys()[-1]
        self.data = dataset[variable]
        self.dataset = dataset
        self.months = self.sort_data(self.data)
        self.months_dict = self._make_months_dict()

    def sort_data(self, data):
        months = self._month_names
        year = self.start.year
        m = 0
        months_data = []
        for month in data:
            months_data.append(MonthMean(month, months[m], year))
            if m == 11:
                m = 0
                year += 1
            else:
                m += 1
        return months_data

    def _make_months_dict(self):
        months = self._month_names
        months_dict = {}
        for month in months:
            months_dict[month] = [
                month_data for month_data in self.months if
                month_data.month == month]
        return months_dict

    def _filter_months(self, month_names, start=None, end=None):
        start = start if start else self.start.year
        end = end if end else self.end.year
        months = [self.months_dict[month] for month in month_names]
        filtered_months = []
        years = range(start, end+1)
        for month in months:
            filtered_months += [year.data for year in month if year.year in
                                years]
        return np.array(filtered_months)

    def before_data_base_map(
            self, llcrnrlon=None, llcrnrlat=None, urcrnrlon=None,
            urcrnrlat=None, llcrnrx=None, llcrnry=None, urcrnrx=None,
            urcrnry=None, width=None, height=None, projection='cyl',
            resolution='c', area_thresh=None, rsphere=6370997.0, ellps=None,
            lat_ts=None, lat_1=None, lat_2=None, lat_0=None, lon_0=None,
            lon_1=None, lon_2=None, o_lon_p=None, o_lat_p=None, k_0=None,
            no_rot=False, suppress_ticks=True, satellite_height=35786000,
            boundinglat=None, fix_aspect=True, anchor='C', celestial=False,
            round=False, epsg=None, ax=None):
        base_map = Basemap(
            llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat, llcrnrx, llcrnry,
            urcrnrx, urcrnry, width, height, projection, resolution,
            area_thresh, rsphere, ellps, lat_ts, lat_1, lat_2, lat_0, lon_0,
            lon_1, lon_2, o_lon_p, o_lat_p, k_0, no_rot, suppress_ticks,
            satellite_height, boundinglat, fix_aspect, anchor, celestial,
            round, epsg, ax)
        fig = plt.figure()
        # Adjust Map to be in Center of Map
        fig.subplots_adjust(top=0.85)
        ax = fig.add_subplot(111)
        # Make the Map have Coastlines
        base_map.drawcoastlines()

        return base_map, ax

    def after_data_base_map(self, base_map, mer_labels=None, par_labels=None,
                            mer_range=np.arange(-180, 180, 20),
                            par_range=np.arange(-90, 90, 5)):
        base_map.colorbar()
        # Plot Meridian/Longitude Lines 20 Degrees Apart
        # Label on the Left and Bottom of Graph
        base_map.drawmeridians(mer_range, labels=mer_labels)
        # Plot Latitude Lines Starting From 60 Degrees to 90, 5 Degrees Apart
        # Label on the Left
        base_map.drawparallels(par_range, labels=par_labels)

    def make_contour_map(
            self, months=None, start=None, end=None, lon=None, lat=None,
            llcrnrlon=None, llcrnrlat=None, urcrnrlon=None,
            urcrnrlat=None, llcrnrx=None, llcrnry=None, urcrnrx=None,
            urcrnry=None, width=None, height=None, projection='cyl',
            resolution='c', area_thresh=None, rsphere=6370997.0, ellps=None,
            lat_ts=None, lat_1=None, lat_2=None, lat_0=None, lon_0=None,
            lon_1=None, lon_2=None, o_lon_p=None, o_lat_p=None, k_0=None,
            no_rot=False, suppress_ticks=True, satellite_height=35786000,
            boundinglat=None, fix_aspect=True, anchor='C', celestial=False,
            round=False, epsg=None, ax=None,
            mer_labels=[False, False, False, False],
            par_labels=[False, False, False, False],
            mer_range=np.arange(-180, 180, 20),
            par_range=np.arange(-90, 90, 5), fill=True
            ):
        start = start if start else self.start.year
        end = end if end else self.end.year
        if not months:
            months = self._month_names
        filtered_months = self._filter_months(months, start, end)
        data = reduce(lambda a, b: a + b, filtered_months)/len(filtered_months)
        # Set Up Latitude and Longitude into Arrays that can be Mapped
        lon = lon if lon is not None else self.lon
        lat = lat if lat is not None else self.lat
        lon, lat = np.meshgrid(lon, lat)
        # Set Latitiude and Longitude to Mappable
        base_map, ax = self.before_data_base_map(
            llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat, llcrnrx, llcrnry,
            urcrnrx, urcrnry, width, height, projection, resolution,
            area_thresh, rsphere, ellps, lat_ts, lat_1, lat_2, lat_0, lon_0,
            lon_1, lon_2, o_lon_p, o_lat_p, k_0, no_rot, suppress_ticks,
            satellite_height, boundinglat, fix_aspect, anchor, celestial,
            round, epsg, ax)
        x, y = base_map(lon, lat)
        # Plot Data as Filled in Contours
        if fill:
            base_map.contourf(x, y, data)
        else:
            base_map.contour(x, y, data)
        self.after_data_base_map(base_map, mer_labels, par_labels, mer_range,
                                 par_range)
        # Title Plot
        if len(months) > 1:
            months_title = months[0] + '-' + months[-1]
        else:
            months[0]
        years = range(start, end + 1)
        years_title = str(years[0]) + '-' + str(years[-1])
        ax.set_title(
            '%s, %s in %s, %s' % (
                self.data.var_desc, self.data.units, months_title, years_title)
            )
        plt.show()