Source code for pipeline.infrastructure.displays.plotweather

import math

import matplotlib
import matplotlib.pyplot as plt
import numpy as np

import pipeline.infrastructure as infrastructure
from pipeline.infrastructure import casa_tools
from pipeline.infrastructure.utils.conversion import mjd_seconds_to_datetime

LOG = infrastructure.get_logger(__name__)


[docs]def plot_weather(vis='', figfile='', station=[], help=False): """ Compiles and plots the major weather parameters for the specified ms. Station can be a single integer or integer string, or a list of two integers. The default empty list means to plot all data from up to 2 of the stations present in the data. The default plot file name will be 'vis'.weather.png. """ if help: LOG.info("plot_weather(vis='', figfile='', station=[])") LOG.info(" Plots pressure, temperature, relative humidity, wind speed and direction.") LOG.info("Station can be a single integer or integer string, or a list of two integers.") LOG.info("The default empty list means to plot the data form up to 2 of the stations") LOG.info("present in the data. The default plot file name will be 'vis'.weather.png.") return myfontsize = 8 try: # Fetch data from weather table in MS. with casa_tools.TableReader(vis+"/WEATHER") as table: available_cols = table.colnames() mjdsec = table.getcol('TIME') pressure = table.getcol('PRESSURE') relative_humidity = table.getcol('REL_HUMIDITY') temperature = table.getcol('TEMPERATURE') # Nobeyama does not have DEW_POINT and NS_WX_STATION_ID dew_point = table.getcol('DEW_POINT') if 'DEW_POINT' in available_cols else None wind_direction = (180 / math.pi) * table.getcol('WIND_DIRECTION') wind_speed = table.getcol('WIND_SPEED') stations = table.getcol('NS_WX_STATION_ID') if 'NS_WX_STATION_ID' in available_cols else [] except: LOG.info("Could not open WEATHER table. Did you importasdm with asis='*'?") return mjdsec1 = mjdsec vis = vis.split('/')[-1] unique_stations = np.unique(stations) try: with casa_tools.TableReader(vis + '/ASDM_STATION') as table: station_names = table.getcol('name') except: LOG.info("Could not open ASDM_STATION table. The Station IDs (instead of Names) will be used.") station_names = None unique_station_names = [] for station_id in unique_stations: station_name = str(station_id) if station_names is not None: if any([wx_prefix.lower() in station_names[station_id].lower() for wx_prefix in ['WSTB', 'Meteo', 'OSF']]): station_name = station_names[station_id].replace('Meteo', '') unique_station_names.append(station_name) if station: if isinstance(station, int): if station not in unique_stations: LOG.info("Station %d is not in the data. Present are: " % station, unique_stations) return unique_stations = [station] elif isinstance(station, list): if len(station) > 2: LOG.info("Only 2 stations can be overlaid.") return if station[0] not in unique_stations: LOG.info("Station %d is not in the data. Present are: " % station[0], unique_stations) return if station[1] not in unique_stations: LOG.info("Station %d is not in the data. Present are: " % station[1], unique_stations) return unique_stations = station elif isinstance(station, str): if station.isdigit(): if int(station) not in unique_stations: LOG.info("Station %s is not in the data. Present are: " % station, unique_stations) return unique_stations = [int(station)] else: LOG.info("Invalid station ID, it must be an integer, or list of integers.") return if len(unique_stations) > 1: first_station_rows = np.where(stations == unique_stations[0])[0] second_station_rows = np.where(stations == unique_stations[1])[0] pressure2 = pressure[second_station_rows] relative_humidity2 = relative_humidity[second_station_rows] temperature2 = temperature[second_station_rows] dew_point2 = dew_point[second_station_rows] if dew_point is not None else None wind_direction2 = wind_direction[second_station_rows] wind_speed2 = wind_speed[second_station_rows] mjdsec2 = mjdsec[second_station_rows] pressure = pressure[first_station_rows] relative_humidity = relative_humidity[first_station_rows] temperature = temperature[first_station_rows] dew_point = dew_point[first_station_rows] if dew_point is not None else None wind_direction = wind_direction[first_station_rows] wind_speed = wind_speed[first_station_rows] mjdsec1 = mjdsec[first_station_rows] if np.mean(temperature2) > 100: # convert to Celsius temperature2 -= 273.15 if dew_point2 is not None and np.mean(dew_point2) > 100: dew_point2 -= 273.15 if np.mean(temperature) > 100: # convert to Celsius temperature -= 273.15 if dew_point is not None and np.mean(dew_point) > 100: dew_point -= 273.15 if dew_point is not None and np.mean(dew_point) == 0: # assume it is not measured and use NOAA formula to compute from humidity: dew_point = ComputeDewPointCFromRHAndTempC(relative_humidity, temperature) if np.mean(relative_humidity) < 0.001: if dew_point is None or np.count_nonzero(dew_point) == 0: # dew point is all zero so it was not measured, so cap the rH at small non-zero value relative_humidity = 0.001 * np.ones(len(relative_humidity)) else: LOG.info("Replacing zeros in relative humidity with value computed from dew point and temperature.") dew_point_wvp = computeWVP(dew_point) ambient_wvp = computeWVP(temperature) LOG.info("dWVP=%f, aWVP=%f" % (dew_point_wvp[0], ambient_wvp[0])) relative_humidity = 100*(dew_point_wvp/ambient_wvp) # take timerange from OBSERVATION table if there is only one unique timestamp if len(np.unique(mjdsec)) == 1: with casa_tools.TableReader(vis+"/OBSERVATION") as table: timerange = table.getcol('TIME_RANGE') obs_timerange = [np.min(timerange), np.max(timerange)] manual_xlim = matplotlib.dates.date2num(mjd_seconds_to_datetime(obs_timerange)) do_manual_xlim = True else: manual_xlim = None do_manual_xlim = False mysize = 'small' plt.clf() adesc = plt.subplot(321) myhspace = 0.25 mywspace = 0.25 markersize = 3 plt.subplots_adjust(hspace=myhspace, wspace=mywspace) plt.title(vis) list_of_date_times = mjd_seconds_to_datetime(mjdsec1) timeplot = matplotlib.dates.date2num(list_of_date_times) plt.plot_date(timeplot, pressure, markersize=markersize) if len(unique_stations) > 1: list_of_date_times = mjd_seconds_to_datetime(mjdsec2) timeplot2 = matplotlib.dates.date2num(list_of_date_times) plt.plot_date(timeplot2, pressure2, markersize=markersize, color='r') if do_manual_xlim is True: plt.xlim(manual_xlim) resizeFonts(adesc, myfontsize) plt.ylabel('Pressure (mb)', size=mysize) adesc.xaxis.set_major_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 30)))) adesc.xaxis.set_minor_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 10)))) adesc.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%H:%M')) adesc.fmt_xdata = matplotlib.dates.DateFormatter('%H:%M') RescaleXAxisTimeTicks(plt.xlim(), adesc) adesc.xaxis.grid(True, which='major') adesc.yaxis.grid(True, which='major') adesc = plt.subplot(322) plt.plot_date(timeplot, temperature, markersize=markersize) if len(unique_stations) > 1: list_of_date_times = mjd_seconds_to_datetime(mjdsec2) timeplot2 = matplotlib.dates.date2num(list_of_date_times) plt.plot_date(timeplot2, temperature2, markersize=markersize, color='r') if do_manual_xlim is True: plt.xlim(manual_xlim) resizeFonts(adesc, myfontsize) plt.ylabel('Temperature (C)', size=mysize) adesc.xaxis.set_major_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 30)))) adesc.xaxis.set_minor_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 10)))) adesc.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%H:%M')) adesc.fmt_xdata = matplotlib.dates.DateFormatter('%H:%M') RescaleXAxisTimeTicks(plt.xlim(), adesc) adesc.xaxis.grid(True, which='major') adesc.yaxis.grid(True, which='major') colors = ['b', 'r'] xinc = 0.4 labxstart = 0.5-(min(len(unique_station_names), 2)-1.0)*xinc/2.0 for idx, station_name in enumerate(unique_station_names): if idx > 1: continue y0 = 1.05 plt.text(labxstart+idx*xinc, y0, station_name, color=colors[idx], transform=adesc.transAxes, ha='center') adesc = plt.subplot(323) plt.plot_date(timeplot, relative_humidity, markersize=markersize) if len(unique_stations) > 1: list_of_date_times = mjd_seconds_to_datetime(mjdsec2) timeplot2 = matplotlib.dates.date2num(list_of_date_times) plt.plot_date(timeplot2, relative_humidity2, markersize=markersize, color='r') if do_manual_xlim is True: plt.xlim(manual_xlim) resizeFonts(adesc, myfontsize) plt.ylabel('Relative Humidity (%)', size=mysize) adesc.xaxis.set_major_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 30)))) adesc.xaxis.set_minor_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 10)))) adesc.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%H:%M')) adesc.fmt_xdata = matplotlib.dates.DateFormatter('%H:%M') RescaleXAxisTimeTicks(plt.xlim(), adesc) adesc.xaxis.grid(True, which='major') adesc.yaxis.grid(True, which='major') pid = 4 if dew_point is not None: adesc = plt.subplot(3, 2, pid) plt.plot_date(timeplot, dew_point, markersize=markersize) if len(unique_stations) > 1: list_of_date_times = mjd_seconds_to_datetime(mjdsec2) timeplot2 = matplotlib.dates.date2num(list_of_date_times) plt.plot_date(timeplot2, dew_point2, markersize=markersize, color='r') if do_manual_xlim is True: plt.xlim(manual_xlim) resizeFonts(adesc, myfontsize) plt.ylabel('Dew point (C)', size=mysize) adesc.xaxis.set_major_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 30)))) adesc.xaxis.set_minor_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 10)))) adesc.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%H:%M')) adesc.fmt_xdata = matplotlib.dates.DateFormatter('%H:%M') RescaleXAxisTimeTicks(plt.xlim(), adesc) adesc.xaxis.grid(True, which='major') adesc.yaxis.grid(True, which='major') pid += 1 adesc = plt.subplot(3, 2, pid) plt.plot_date(timeplot, wind_speed, markersize=markersize) if len(unique_stations) > 1: list_of_date_times = mjd_seconds_to_datetime(mjdsec2) timeplot2 = matplotlib.dates.date2num(list_of_date_times) plt.plot_date(timeplot2, wind_speed2, markersize=markersize, color='r') if do_manual_xlim is True: plt.xlim(manual_xlim) resizeFonts(adesc, myfontsize) # Get date of observation. date_string = mjd_seconds_to_datetime(mjdsec[0:])[0].strftime('%Y-%m-%d') xlabel = 'Universal Time (%s)' % date_string plt.xlabel(xlabel, size=mysize) plt.ylabel('Wind speed (m/s)', size=mysize) adesc.xaxis.set_major_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 30)))) adesc.xaxis.set_minor_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 10)))) adesc.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%H:%M')) adesc.fmt_xdata = matplotlib.dates.DateFormatter('%H:%M') RescaleXAxisTimeTicks(plt.xlim(), adesc) adesc.xaxis.grid(True, which='major') adesc.yaxis.grid(True, which='major') pid += 1 adesc = plt.subplot(3, 2, pid) plt.xlabel(xlabel, size=mysize) plt.ylabel('Wind direction (deg)', size=mysize) plt.plot_date(timeplot, wind_direction, markersize=markersize) if len(unique_stations) > 1: list_of_date_times = mjd_seconds_to_datetime(mjdsec2) timeplot2 = matplotlib.dates.date2num(list_of_date_times) plt.plot_date(timeplot2, wind_direction2, markersize=markersize, color='r') if do_manual_xlim is True: plt.xlim(manual_xlim) resizeFonts(adesc, myfontsize) adesc.xaxis.set_major_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 30)))) adesc.xaxis.set_minor_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 10)))) adesc.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%H:%M')) adesc.fmt_xdata = matplotlib.dates.DateFormatter('%H:%M') RescaleXAxisTimeTicks(plt.xlim(), adesc) adesc.xaxis.grid(True, which='major') adesc.yaxis.grid(True, which='major') if len(figfile) < 1: weather_file = vis+'.weather.png' else: weather_file = figfile plt.savefig(weather_file) plt.draw() LOG.info("Wrote file = %s" % weather_file)
[docs]def ComputeDewPointCFromRHAndTempC(relativeHumidity, temperature): """ inputs: relativeHumidity in percentage, temperature in C output: in degrees C Uses formula from http://en.wikipedia.org/wiki/Dew_point#Calculating_the_dew_point """ es = 6.112*np.exp(17.67*temperature/(temperature+243.5)) e = relativeHumidity*0.01*es dewPoint = 243.5*np.log(e/6.112)/(17.67-np.log(e/6.112)) return dewPoint
[docs]def computeWVP(d): """ This simply converts the specified temperature (in Celsius) to water vapor pressure, which can be used to estimate the relative humidity from the measured dew point. """ # d is in Celsius t = d+273.15 w = np.exp(-6096.9385/t + 21.2409642 - 2.711193e-2*t + 1.673952e-5*t**2 + 2.433502*np.log(t)) return w
[docs]def resizeFonts(adesc, fontsize): """ Plotting utility routine """ y_format = matplotlib.ticker.ScalarFormatter(useOffset=False) adesc.yaxis.set_major_formatter(y_format) adesc.xaxis.set_major_formatter(y_format) plt.setp(adesc.get_xticklabels(), fontsize=fontsize) plt.setp(adesc.get_yticklabels(), fontsize=fontsize)
[docs]def RescaleXAxisTimeTicks(xlim, adesc): """ Plotting utility routine """ if xlim[1] - xlim[0] < 10/1440.: adesc.xaxis.set_major_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 1)))) adesc.xaxis.set_minor_locator(matplotlib.dates.SecondLocator(bysecond=list(range(0, 60, 30)))) elif xlim[1] - xlim[0] < 0.5/24.: adesc.xaxis.set_major_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 5)))) adesc.xaxis.set_minor_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 1)))) elif xlim[1] - xlim[0] < 1/24.: adesc.xaxis.set_major_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 10)))) adesc.xaxis.set_minor_locator(matplotlib.dates.MinuteLocator(byminute=list(range(0, 60, 2))))