Commit 1549df81 authored by Renzo Marini's avatar Renzo Marini
Browse files

modo alt, llamadas por punto y area, fixes menores

parent 0387c1fd
# -*- coding: utf-8 -*-
import argparse
import copy
import pygrib
import numpy as np
import pickle
......@@ -12,104 +13,163 @@ def daterange(start_date, end_date):
yield start_date + timedelta(n)
# El grib es cargado en diccionarios anidados tal como se muestra a continuación:
# self.param_data[modelo][parámetro][día][tiempo inicial][tiempo de pronóstico].
# self.model["data"][día]["normal" | "alt"][tiempo de pronóstico].
# Cada elemento del diccionario final tiene un array con los datos de todas las perturbaciones.
# Acerca de los mensajes grib, se asume que contienen información de un mismo centro, parámetro,
# área, tiempo inicial y steps.
class GribParser:
def __init__(self, grib_path):
def __init__(self, grib_path, alt):
if not os.path.exists('Pickles'):
os.makedirs('Pickles')
pickle_path = "Pickles/" + grib_path.split("/")[-1][:-5] + '.pickle'
if os.path.isfile(pickle_path):
with open(pickle_path, 'rb') as handle:
self.param_data = pickle.load(handle)
return
self.pickle_path = "Pickles/" + grib_path.split("/")[-1][:-5] + '.pickle'
if os.path.isfile(self.pickle_path):
with open(self.pickle_path, 'rb') as handle:
self.model = pickle.load(handle)
if alt and not self.model["metadata"]["alt"]:
self.create_non_cumulative_values()
self.dump()
else:
self.model = {}
self.model["data"] = {}
grbs = pygrib.open(grib_path)
for grb in grbs:
data, lats, lons = grb.data()
date = str(grb["dataDate"])
if date not in self.model["data"]:
self.model["data"][date] = {}
self.model["data"][date]["normal"] = {}
forecast_time = str(grb["stepRange"])
param_by_date = self.model["data"][date]["normal"]
if forecast_time not in param_by_date:
param_by_date[forecast_time] = data[np.newaxis, :, :]
else:
param_by_date[forecast_time] = \
np.concatenate((param_by_date[forecast_time], data[np.newaxis, :, :]))
self.load_metadata(grb, alt)
if alt:
self.create_non_cumulative_values()
self.dump()
self.define_grid()
# Define el conjunto de puntos que son considerados como válidos dentro del grid.
# Inicialmente, son todos los puntos que forman parte del área definida en los gribs.
# Luego puede reducirse si se quiere trabajar con un subgrid.
def define_grid(self, borders = []):
if not borders:
self.first_point = [self.model["metadata"]["firstLat"], self.model["metadata"]["firstLon"]]
self.last_point = [self.model["metadata"]["lastLat"], self.model["metadata"]["lastLon"]]
borders = [p for p in self.first_point + self.last_point]
borders = [int(x) for x in borders]
grib_info = pygrib.fromstring(self.model["metadata"]["additional_info"])
self.i_step = grib_info["iDirectionIncrementInDegrees"]
self.j_step = grib_info["jDirectionIncrementInDegrees"]
i_range = range(min(borders[0], borders[2]), max(borders[0], borders[2]) + 1, int(self.i_step))
j_range = range(min(borders[1], borders[3]), max(borders[1], borders[3]) + 1, int(self.j_step))
self.grid = [[i, j] for i in i_range for j in j_range]
def load_step_range(self, step_range, initial_time):
# Solo muestra los steps, sin el initial_time. Se usa en la primer fila del CSV.
formatted_step_range = []
if len(step_range[-1].split("-")) == 1:
# El formato del step range es -> step
new_step_range = copy.deepcopy(step_range)
new_step_range.sort(key=int)
formatted_step_range = new_step_range
else:
# El formato del step range es -> initial_time-step, salvo el primer step que es initial_time
formatted_step_range = []
for step in step_range:
step = step.split("-")
if len(step) == 1:
formatted_step_range.append(step[0])
else:
formatted_step_range.append(step[1])
formatted_step_range.sort(key=int)
new_step_range = [formatted_step_range[0]] + [initial_time + "-" + x for x in formatted_step_range[1:]]
self.model["metadata"]["step_range"] = new_step_range
self.model["metadata"]["formatted_step_range"] = formatted_step_range
self.param_data = {}
self.centers = {}
grbs = pygrib.open(grib_path)
for grb in grbs:
data, lats, lons = grb.data()
center = grb["centre"]
if center not in self.param_data:
self.param_data[center] = {}
param = grb["shortName"]
if param not in self.param_data[center]:
self.param_data[center][param] = {}
date = str(grb["dataDate"])
if date not in self.param_data[center][param]:
self.param_data[center][param][date] = {}
initial_time = str(grb["dataTime"])
if initial_time not in self.param_data[center][param][date]:
self.param_data[center][param][date][initial_time] = {}
forecast_time = str(grb["stepRange"])
if forecast_time not in self.param_data[center][param][date][initial_time]:
self.param_data[center][param][date][initial_time][forecast_time] = data[np.newaxis, :, :]
else:
self.param_data[center][param][date][initial_time][forecast_time] = \
np.concatenate((self.param_data[center][param][date][initial_time][forecast_time], data[np.newaxis, :, :]))
self.param_data["metadata"] = {}
self.param_data["metadata"]["center"] = center
self.param_data["metadata"]["parameter"] = param
self.param_data["metadata"]["step_range"] = self.param_data[center][param][date][initial_time].keys()
self.param_data["metadata"]["firstLat"] = grb["latitudeOfFirstGridPointInDegrees"]
self.param_data["metadata"]["firstLon"] = grb["longitudeOfFirstGridPointInDegrees"]
self.param_data["metadata"]["lastLat"] = grb["latitudeOfLastGridPointInDegrees"]
self.param_data["metadata"]["lastLon"] = grb["longitudeOfLastGridPointInDegrees"]
def load_metadata(self, grb, alt):
self.model["metadata"] = {}
self.model["metadata"]["center"] = grb["centre"]
self.model["metadata"]["parameter"] = grb["shortName"]
self.model["metadata"]["initial_time"] = str(grb["dataTime"])
self.load_step_range(self.model["data"][str(grb["dataDate"])]["normal"].keys(), str(grb["dataTime"]))
self.model["metadata"]["firstLat"] = grb["latitudeOfFirstGridPointInDegrees"]
self.model["metadata"]["firstLon"] = grb["longitudeOfFirstGridPointInDegrees"]
self.model["metadata"]["lastLat"] = grb["latitudeOfLastGridPointInDegrees"]
self.model["metadata"]["lastLon"] = grb["longitudeOfLastGridPointInDegrees"]
self.model["metadata"]["alt"] = alt
self.model["metadata"]["additional_info"] = grb.tostring()
with open(pickle_path, 'wb') as handle:
pickle.dump(dict(self.param_data), handle, protocol=pickle.HIGHEST_PROTOCOL)
def dump(self):
with open(self.pickle_path, 'wb') as handle:
pickle.dump(dict(self.model), handle, protocol=pickle.HIGHEST_PROTOCOL)
def valid_input(self, label, value):
if label == "latitud":
firstCoord = self.param_data["metadata"]["firstLat"]
lastCoord = self.param_data["metadata"]["lastLat"]
else:
firstCoord = self.param_data["metadata"]["firstLon"]
lastCoord = self.param_data["metadata"]["lastLon"]
if value < min(firstCoord, lastCoord) or value > max(firstCoord, lastCoord):
print("ERROR: ", label, " fuera del rango permitido\n")
print("El rango permitido es de ", min(firstCoord, lastCoord), " a ", max(firstCoord, lastCoord))
def create_non_cumulative_values(self):
self.model["metadata"]["alt"] = True
step_range = self.model["metadata"]["step_range"]
data = self.model["data"]
for day in data:
alt = {}
alt[step_range[0]] = np.copy(data[day]["normal"][step_range[0]])
for step, next_step in zip(step_range[0:], step_range[1:]):
previous_array = data[day]["normal"][step]
alt[next_step] = np.copy(data[day]["normal"][next_step])
for i in range(0, previous_array.shape[0]):
alt[next_step][i, :, :] -= previous_array[i, :, :]
data[day]["alt"] = alt
def valid_input(self, point):
if point not in self.grid:
print "\nERROR: el punto elegido está fuera del área\n"
print " - El primer punto del área es", self.first_point
print " - El último punto del área es", self.last_point
print " - El incremento es de", self.i_step, "x", self.j_step, "grados"
return False
return True
return True
def dump_percentile(self, n, lat, lon, initial_time = "0"):
if not self.valid_input("latitud", lat) or not self.valid_input("longitud", lon):
return
lat_rel = int(lat - self.param_data["metadata"]["firstLat"])
lon_rel = int(lon - self.param_data["metadata"]["firstLon"])
center = self.param_data["metadata"]["center"]
parameter = self.param_data["metadata"]["parameter"]
step_range = [initial_time] + [x.split("-")[1] for x in self.param_data["metadata"]["step_range"] if len(x) > 1]
step_range.sort(key=int)
dates = self.param_data[center][parameter].keys()
def dump_percentile(self, n, point, alt):
param_data = self.model["data"]
initial_time = self.model["metadata"]["initial_time"]
lat_rel = int(point[0] - self.model["metadata"]["firstLat"])
lon_rel = int(point[1] - self.model["metadata"]["firstLon"])
center = self.model["metadata"]["center"]
parameter = self.model["metadata"]["parameter"]
step_range = self.model["metadata"]["step_range"]
formatted_step_range = self.model["metadata"]["formatted_step_range"]
dates = param_data.keys()
dates.sort()
date_format = "%Y%m%d"
new_rows = []
not_found_field = "NaN"
start_date = datetime.strptime(dates[0], date_format).date()
end_date = datetime.strptime(dates[-1], date_format).date() + timedelta(days=1)
for curr_date_object in daterange(start_date, end_date):
curr_date = curr_date_object.strftime(date_format)
row = [curr_date]
if curr_date in dates:
data = self.param_data[center][parameter][curr_date][initial_time]
data = param_data[curr_date]
if alt:
data = data["alt"]
else:
data = data["normal"]
for step in step_range:
step_key = '0' if step == '0' else initial_time + "-" + step
if step_key in data:
row.append(np.percentile(data[step_key][:, int(lat_rel), int(lon_rel)], n))
if step in data:
row.append(np.percentile(data[step][:, int(lat_rel), int(lon_rel)], n))
else:
row.append("-")
row.append(not_found_field)
else:
row = row + ["-" for x in step_range]
row = row + [not_found_field for x in step_range]
new_rows.append(row)
if not os.path.exists('CSV'):
os.makedirs('CSV')
if alt:
parameter = parameter + "[alt]"
csv_path = "CSV/" + center + "_" + parameter + "_per" + str(n) + "_" \
+ initial_time + "hs_" + "lat" + str(lat) + "_lon" + str(lon) + ".csv"
first = [""] + step_range
+ initial_time + "hs_" + "lat" + str(point[0]) + "_lon" + str(point[1]) + ".csv"
first = [""] + formatted_step_range
if os.path.isfile(csv_path):
with open(csv_path, 'rb') as readFile:
reader = csv.reader(readFile)
......@@ -128,14 +188,31 @@ class GribParser:
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('--alt', dest='alt', action='store_true')
parser.add_argument('--path', type=str, required=True,
help='Ingrese el path del grib')
group = parser.add_mutually_exclusive_group(required=True)
group.add_argument('--point', type=int, nargs=2,
help='Ingrese la latitud y longitud del punto')
group.add_argument('--area', type=int, nargs=4,
help='Ingrese la latitud y longitud del primer y último punto del área')
group.add_argument('--all', dest='all', action='store_true')
parser.add_argument('--per', type=int, required=True,
help='Ingrese el percentil')
parser.add_argument('--lat', type=int, required=True,
help='Ingrese la latitud')
parser.add_argument('--lon', type=int, required=True,
help='Ingrese la longitud')
args = parser.parse_args()
grib_parser = GribParser(args.path)
grib_parser.dump_percentile(args.per, args.lat, args.lon)
grib_parser = GribParser(args.path, args.alt)
if args.point:
if grib_parser.valid_input(args.point):
grib_parser.dump_percentile(args.per, args.point, args.alt)
else:
if args.area:
area = args.area
elif args.all:
first_point = [grib_parser.model["metadata"]["firstLat"], grib_parser.model["metadata"]["firstLon"]]
last_point = [grib_parser.model["metadata"]["lastLat"], grib_parser.model["metadata"]["lastLon"]]
area = first_point + last_point
if grib_parser.valid_input(area[:2]) and grib_parser.valid_input(area[2:]):
grib_parser.define_grid(area)
for point in grib_parser.grid:
grib_parser.dump_percentile(args.per, point, args.alt)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment