from station import Station
from datetime import date, timedelta, datetime
import matplotlib.pyplot as plt
import pathos.pools as pools
import numpy as np
import warnings
from os import path
from google.cloud import bigquery
# Gathers weather data for a given location and time range
class Record:
# IMPORTANT: multithreading on windows requires code to be inside of
# if __name__ == '__main__':
# ...
# otherwise it will recursively run the entire code
multithreaded = True
date_format = "%Y%m%d" # YYYYMMDD
attributes = ['min', # Min temperature F°
'temp', # Mean temperature F°
'max', # Max temperature F°
'wdsp', # Mean wind speed
'mxspd', # Max wind speed (sustained)
#'gust', # Max wind speed (instant) - IGNORED
'prcp', # Total precipitation
'dewp', # Mean dew point F°
'sndp', # Snow depth
'stp', # Mean station pressure
'slp', # Mean sea level pressure
'visib', # Mean visibility
'fog', # Presence of fog
'rain', # Presence of rain
'snow', # Presence of snow
'hail', # Presence of hail
'thunder', # Presence of thunder
'tornado'] # Presence of a tornado
max_station_distance = 50.0 # Distance in km after which a station has no weight
def __init__(self, name, shape, start_date, end_date=None, contour_dist=25, max_stations=12,
interactive=False, save_plot=False):
"""
:param name: The name of the location
:param shape: The (multi)polygon describing the shape of the location
:param start_date: Beginning of the record period (YYYYMMDD)
:param end_date: End of the record period. If None, the day before the current day is used.
:param contour_dist: Range within which weather station data is to be trusted.
:param max_stations: Maximum number of stations used to gather data
:param interactive: Whether or not to display a plot of the stations
:param save_plot: Whether or not to save this plot to disk
"""
self.name = name
self.shape = shape
self.contour_dist = contour_dist
if contour_dist > Record.max_station_distance:
warnings.warn('contour distance of %f is greater than maximum allowed distance of %f' %
(contour_dist, Record.max_station_distance))
self.max_stations = max_stations
self.start_date = datetime.strptime(start_date, Record.date_format).date()
if end_date is None:
self.end_date = date.today() - timedelta(days=1)
else:
self.end_date = datetime.strptime(end_date, Record.date_format).date()
self.interactive = interactive
self.save_plot = save_plot
# Gather the data and build the record
self.collection_info = ""
self.__build()
def __log(self, message):
self.collection_info += "# " + message.replace("\n", "\n# ") + "\n"
def __plot_stations(self, closest_stations):
shapes = self.shape if self.shape.geom_type == 'MultiPolygon' else [self.shape]
for sub_shape in shapes:
x, y = sub_shape.exterior.xy
plt.plot(x, y, c='blue')
for station, distance in closest_stations:
plt.plot(station.longitude, station.latitude, 'o', c='red')
plt.text(station.longitude, station.latitude, station.name + ("\n%.1fkm" % distance))
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Weather stations in ' + self.name)
if self.save_plot:
plt.savefig(self.name + '.png')
if self.interactive:
plt.show()
def __build(self):
self.__log(self.name.upper())
print("Gathering data from " + self.name +
" from " + self.start_date.strftime("%d/%m/%Y") +
" to " + self.end_date.strftime("%d/%m/%Y"))
self.__log("From " + self.start_date.strftime("%d/%m/%Y") +
" to " + self.end_date.strftime("%d/%m/%Y") + "\n")
# Retrieve stations within <contour_dist> that are overlapping with the time range
closest_stations = Station.find_stations_in_geometry(
self.shape,
self.contour_dist,
self.start_date,
self.end_date
)
if len(closest_stations) == 0:
print("Found no stations close enough! Try increasing the search radius.")
return
self.max_stations = min(self.max_stations, len(closest_stations))
# Prune unnecessary stations
print("Using the first %d of the %d stations found within %dkm:" %
(self.max_stations, len(closest_stations), self.contour_dist))
closest_stations = closest_stations[:self.max_stations]
stations, distances = zip(*closest_stations)
# Display them
self.__log("STATIONS:")
base_weights = [Record.distance_weight(distance) for distance in distances]
base_weights = np.array(base_weights)
for (station, distance), weight in zip(closest_stations, base_weights):
message = station.name[:20].ljust(20) + (" at %.2fkm" % distance).ljust(14) + \
"(trust: %.2f%%)" % (weight * 100)
print(message)
self.__log(message)
if self.interactive or self.save_plot:
self.__plot_stations(closest_stations)
# Gather data
self.data = []
if Record.multithreaded:
thread_pool = pools.ProcessPool(4)
date = self.start_date
year = -1
while date <= self.end_date:
if year != date.year:
# Retrieve data from the year <year> for each station
year = date.year
print("\nCollecting data for the year %d" % year)
if Record.multithreaded:
job = lambda station: station.retrieve_obs(year)
all_yearly_data = thread_pool.map(job, stations)
else:
all_yearly_data = [station.retrieve_obs(year) for station in stations]
# Gather daily data
date_key = date.strftime(Record.date_format)
daily_data, daily_trusts = [], []
for yearly_data, trust in zip(all_yearly_data, base_weights):
if yearly_data is None:
continue
if date_key in yearly_data:
daily_data.append(yearly_data[date_key])
daily_trusts.append(trust)
# Case of missing data for the entire day
if len(daily_data) == 0:
print("Got no data for " + str(date))
# Populate attributes
datum = {}
for attribute in Record.attributes:
# Perform a weighted average
value = 0
total_weight = 0
for daily_datum, daily_trust in zip(daily_data, daily_trusts):
if daily_datum[attribute] is None:
continue
value += daily_datum[attribute] * daily_trust
total_weight += daily_trust
if total_weight == 0:
datum[attribute] = None
continue
value = value / total_weight
datum[attribute] = value
# Add datum to record
self.data.append((date, datum))
date += timedelta(days=1)
@staticmethod
def distance_weight(distance):
"""
Computes a coefficient of trust based on the distance between the observations and the
target location.
:param distance: [float] The distance in kilometers
:return: [float] The coefficient of trust (within [0 - 1])
"""
# A plot of the function can be found here: https://i.imgur.com/QjzAT18.png
if distance >= Record.max_station_distance:
return 0
return 1 - (distance / Record.max_station_distance) ** 2
def export_as_csv(self, filepath):
just_size = 11
csvfile = open(filepath, 'w')
# Write the header
csvfile.write(self.collection_info + "\n")
csvfile.write("DATE,".ljust(just_size))
for attribute in Record.attributes:
attribute = attribute.upper() + ("," if attribute != Record.attributes[-1] else "")
csvfile.write(attribute.ljust(just_size))
csvfile.write("\n")
# Write each row
for date, datum in self.data:
csvfile.write((date.strftime(Record.date_format)+ ", ").ljust(just_size))
for attribute in Record.attributes:
if not attribute in datum:
raise Exception('Corrupted record data')
if datum[attribute] is None:
text = "NA"
else:
text = "%.3f" % datum[attribute]
if attribute != Record.attributes[-1]:
text += ", "
csvfile.write(text.ljust(just_size))
csvfile.write("\n")
print("\nSuccesfully written " + filepath + "!")
csvfile.close()
def export_in_bigquerry(self):
access_key = path.realpath("../../access_key.json")
bigquery_client = bigquery.Client.from_service_account_json(access_key)
dataset_ref = bigquery_client.dataset("weather")
table_ref = bigquery_client.get_table(dataset_ref.table("noaa"))
data = []
for date, datum in self.data:
entry = {
"city": self.name,
"date": date.strftime("%Y-%m-%d")
}
for attribute in Record.attributes:
if not attribute in datum:
raise Exception("Corrupted record data")
if datum[attribute] is None:
entry[attribute] = None
else:
entry[attribute] = "%.3f" % datum[attribute]
data.append(entry)
result = bigquery_client.create_rows(table_ref, data)
if result != []:
print("Error storing weather")
print(result)
else:
print("Storage successful")
@staticmethod
def read_from_csv(filepath):
csvfile = open(filepath, 'r')
# Skip header
for line in csvfile:
line = line.rstrip()
if line != "" and line[0] != '#':
break
# Match row names
csv_rownames = line.lower().replace(" ", "").split(',')
true_rownames = ['date'] + Record.attributes
if csv_rownames != true_rownames:
raise Exception('\nHeader mismatch:\n' + str(csv_rownames) +
'\n\nAgainst:\n' + str(true_rownames))
# Read contents
data = []
for line in csvfile:
text_datum = line.rstrip().replace(" ", "").split(',')
date = datetime.strptime(text_datum[0], Record.date_format).date()
values = [None if value == 'NA' else float(value) for value in text_datum[1:]]
datum = {attribute: value for attribute, value in zip(Record.attributes, values)}
data.append((date, datum))
csvfile.close()
return data
