versión inicial

.gitignore

+*.csv
+CSV
+*.pyc
+*.grib
+*.request
+*.pickle

constants.py

+# -*- coding: utf-8 -*-
+
+origin = [str(x) for x in list(range(0, 11))]
+o_map = {}
+o_map[0] = ("BoM", "ammc")
+o_map[1] = ("CMA", "babj")
+o_map[2] = ("CPTEC", "sbsj")
+o_map[3] = ("ECCC", "cwao")
+o_map[4] = ("ECMWF", "ecmf")
+o_map[5] = ("JMA", "rjtd")
+o_map[6] = ("KMA", "rksl")
+o_map[7] = ("Meteo_France", "lfpw")
+o_map[8] = ("NCEP", "kwbc")
+o_map[9] = ("NCMRWF", "dems")
+o_map[10] = ("UKMO", "egrr")
+
+time = ["00", "06", "12", "18"]
+
+type = ["pf", "cf", "fc"]
+type_map = {
+	"pf": "Perturbed Forecast", 
+	"cf": "Control Forecast",
+	"fc": "Forecast"
+	}
+
+levelist = {
+	"pt": ["320"], 
+	"pv": ["2"],
+	"pl": ["1000", "925", "850", "700", "500", "300", "250", "200", "50"]
+	}
+
+step = ["0", "6", "12", "18", "24", "30", "36", "42", "48", "54", "60", "66", "72", "78", "84", "90", "96", "102", "108", "114", "120", "126", "132", "138", "144", "150", "156", "162", "168", "174", "180", "186", "192", "198", "204", "210", "216", "222", "228", "234", "240", "246", "252", "258", "264", "270", "276", "282", "288", "294", "300", "306", "312", "318", "324", "330", "336", "342", "348", "354", "360", "366", "372", "378", "384"]
+
+levtype = ["sfc", "pl", "pv", "pt"]
+lt_map = {
+	"sfc": "Surface",
+	"pl": "Pressure levels",
+	"pv": "Potential vorticity",
+	"pt": "Potential temperature"
+	}
+
+ids = {}
+
+ids["sfc"] = {}
+ids["sfc"]["2 metre dewpoint temperature"] = "168"
+ids["sfc"]["2 metre temperature"] = "167"
+ids["sfc"]["10 metre U wind component"] = "165"
+ids["sfc"]["10 metre V wind component"] = "166"
+ids["sfc"]["Convective available potential energy"] = "59"
+ids["sfc"]["Convective inhibition"] = "228001"
+ids["sfc"]["Field capacity"] = "190170"
+ids["sfc"]["Land-sea mask"] = "172"
+ids["sfc"]["Maximum temperature at 2 metres in the last 6 hours"] = "121"
+ids["sfc"]["Mean sea level pressure"] = "151"
+ids["sfc"]["Minimum temperature at 2 metres in the last 6 hours"] = "122"
+ids["sfc"]["Orography"] = "228002"
+ids["sfc"]["Skin temperature"] = "235"
+ids["sfc"]["Snow Fall water equivalent"] = "228144"
+ids["sfc"]["Snow depth water equivalent"] = "228141"
+ids["sfc"]["Soil Moisture"] = "228039"
+ids["sfc"]["Soil Temperature"] = "228139"
+ids["sfc"]["Soil moisture top 20 cm"] = "228086"
+ids["sfc"]["Soil temperature top 20 cm"] = "228095"
+ids["sfc"]["Sunshine duration"] = "189"
+ids["sfc"]["Surface latent heat flux"] = "147"
+ids["sfc"]["Surface net solar radiation"] = "176"
+ids["sfc"]["Surface net thermal radiation"] = "177"
+ids["sfc"]["Surface pressure"] = "134"
+ids["sfc"]["Surface sensible heat flux"] = "146"
+ids["sfc"]["Top net thermal radiation"] = "179"
+ids["sfc"]["Total Cloud Cover"] = "164"
+ids["sfc"]["Total Precipitation"] = "228"
+ids["sfc"]["Total column water"] = "136"
+ids["sfc"]["Wilting point"] = "228171"
+
+ids["pl"] = {}
+ids["pl"]["Geopotential Height"] = "156"
+ids["pl"]["Specific humidity"] = "133"
+ids["pl"]["Temperature"] = "130"
+ids["pl"]["U component of wind"] = "131"
+ids["pl"]["V component of wind"] = "132"
+
+ids["pv"] = {}
+ids["pv"]["Potential temperature"] = "3"
+ids["pv"]["U component of wind"] = "131"
+ids["pv"]["V component of wind"] = "132"
+
+ids["pt"] = {}
+ids["pt"]["Potencial vorticity"] = "60"

grib_analysis.py

+# -*- coding: utf-8 -*-
+import argparse
+import pygrib
+import numpy as np
+import pickle
+import os, sys
+import csv
+from datetime import timedelta, date, datetime
+
+def daterange(start_date, end_date):
+    for n in range(int ((end_date - start_date).days)):
+        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].
+# 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):
+        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.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"]
+        
+        with open(pickle_path, 'wb') as handle:
+            pickle.dump(dict(self.param_data), 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))
+            return False
+        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()
+        dates.sort()
+        date_format = "%Y%m%d"
+        new_rows = []
+        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]
+                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))
+                    else:
+                        row.append("-")
+            else:
+                row = row + ["-" for x in step_range]
+            new_rows.append(row)
+            
+        if not os.path.exists('CSV'):
+            os.makedirs('CSV')
+        csv_path = "CSV/" + center + "_" + parameter + "_per" + str(n) + "_" \
+                + initial_time + "hs_" + "lat" + str(lat) + "_lon" + str(lon) + ".csv"
+        first = [""] + step_range
+        if os.path.isfile(csv_path):
+            with open(csv_path, 'rb') as readFile:
+                reader = csv.reader(readFile)
+                rows = list(reader)[1:]
+                if int(rows[0][0]) > int(new_rows[-1][0]):
+                    rows = [first] + new_rows + rows
+                else:
+                    rows = [first] + rows + new_rows
+            os.remove(csv_path)
+        else:
+            rows = [first] + new_rows
+        with open(csv_path, 'wb') as csvfile:
+            writer = csv.writer(csvfile)
+            writer.writerows(rows)
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--path', type=str, required=True,
+                    help='Ingrese el path del grib')
+    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)

script.py

+#!/usr/bin/env python
+    # -*- coding: utf-8 -*-
+from ecmwfapi import ECMWFDataServer
+import constants
+import sys
+import datetime
+
+def download(parameters):
+    server = ECMWFDataServer()
+    server.retrieve(parameters)
+
+def print_help():
+    print("[Esto explicaría las distintas opciones para correr el script]")
+    sys.exit()
+    
+def handle_error(reason, parameter = None, levtype = None):
+    print("\nERROR: " + reason + "\n")
+    if parameter is not None:
+        if levtype is not None:
+            if parameter == "param":
+                print("Las variables posibles para levtype = " + levtype + " son:")
+                for variable in constants.ids[levtype]:
+                    print(">> " + variable + " | Id: " + str(constants.ids[levtype][variable]))
+            elif parameter == "levelist":
+                print("Las niveles posibles para levtype = " + levtype + " son:")
+                for nivel in constants.levelist[levtype]:
+                    print(">> " + str(nivel))
+        elif parameter == "origin":
+            for number in constants.o_map:
+                print(str(number) + " -> " + constants.o_map[number][0])
+        elif parameter == "type":
+            print("Los valores posibles son " + ", ".join([str(x) for x in getattr(constants, parameter)][:-1]) + " y " + str(getattr(constants, parameter)[-1]))
+            for abbr in constants.type_map:
+                print(str(abbr) + " -> " + constants.type_map[abbr])
+        elif parameter == "levtype":
+            print("Los valores posibles son " + ", ".join([str(x) for x in getattr(constants, parameter)][:-1]) + " y " + str(getattr(constants, parameter)[-1]))
+            for abbr in constants.lt_map:
+                print(str(abbr) + " -> " + constants.lt_map[abbr])
+        elif parameter == "number":
+            print("Los valores posibles son números del 1 al 50. Opcionalmente, puede utilizarse el valor 'all' como atajo para designar todos.")
+        elif parameter == "origin_time":
+            print("Los valores que toma origin_time son listas con el formato [origin time time ...]. Por ejemplo: 1 00 06 12.")
+            print("Para origin, los identificadores están dados por las siguientes asociaciones:")
+            for number in constants.o_map:
+                print(str(number) + " -> " + constants.o_map[number][0])
+            print("Los valores posibles para time son 00, 06, 12 y 18")
+        else:
+            print("Los valores posibles son " + ", ".join([str(x) for x in getattr(constants, parameter)][:-1]) + " y " + str(getattr(constants, parameter)[-1]))
+        print("\n")
+            
+    print("Para obtener ayuda, correr el script con la opción --help")
+    sys.exit()
+
+if __name__ == "__main__":
+
+    # Parsear parámetros
+    required_parameters = ["origin_time", "levtype", "start_date", "end_date", "param", "step", "type"]
+    optional_parameters = ["grid", "area"]
+    conditional_parameters = ["number", "levelist"]
+    parameters = {}
+    for par in sys.argv[1:]:
+        if par[:2] == "--":
+            curr_par = par[2:].lower()
+            if curr_par == "origin_time":
+                if curr_par in parameters:
+                    parameters[curr_par].append([])
+                else:
+                    parameters[curr_par] = [[]]
+            elif curr_par == "help":
+                print_help()
+            else:
+                parameters[curr_par] = []
+        elif curr_par == "origin_time":
+            parameters[curr_par][-1].append(par.lower())
+        else:
+            parameters[curr_par].append(par.lower())
+
+    # Validar parametros en el caso general
+    all_parameters = required_parameters + optional_parameters + conditional_parameters
+    for par in parameters:
+        if par not in all_parameters:
+            handle_error("El parámetro --" + par + " es inválido")     
+        elif par in ["start_date", "end_date"]:
+            parameters[par][0] = parameters[par][0].replace("/", '-')
+            partition = parameters[par][0].split("-")
+            if len(partition) != 3 or len(partition[0]) != 4:
+                handle_error("El formato de " + par + " es AAAA-MM-DD")  
+        elif par in ["step", "type", "levtype"]:
+            for curr_par in parameters[par]:
+                if curr_par not in [str(x) for x in getattr(constants, par)]:
+                    handle_error("El valor ingresado para el parámetro " + par + " es inválido.", par)
+        elif par == "origin_time":
+            for curr_origin_time in parameters[par]:
+                if curr_origin_time[0] in constants.origin:
+                    for time in curr_origin_time[1:]:
+                        if time not in constants.time:
+                             handle_error("El valor ingresado para el parámetro " + par + " es inválido.", par)
+                else:
+                    handle_error("El valor ingresado para el parámetro " + par + " es inválido.", par)
+        elif par == "param":
+            for curr_par in parameters["param"]:
+            #                if curr_par not in [str(id) for levtype_variables in constants.ids.values() for id in levtype_variables.values()]:
+                if curr_par not in [id for id in constants.ids[parameters["levtype"][0]].values()]:
+                    handle_error(curr_par + " no es un parámetro posible en type of level " + parameters["levtype"][0], par, parameters["levtype"][0])
+                                        
+    # Chequear que todos los parámetros requiridos fueron recibidos
+    for par in required_parameters:
+        if par not in parameters:
+            if par in ["start_date", "end_date"]:
+                handle_error("El parámetro " + par + " es obligatorio. El formato de " + par + " es AAAA-MM-DD")
+            elif par == "param":
+                handle_error("El parámetro " + par + " es obligatorio.", par, parameters["levtype"][0])
+            else:
+                handle_error("El parametro " + par + " es obligatorio.", par)
+        elif len(parameters[par]) == 0:
+            handle_error("No se proporcionó ningún valor para " + par)
+    
+    # Validar levelist en los casos en que type of level es distinto de Surface
+    if parameters["levtype"][0] != "sfc":
+        if "levelist" not in parameters:
+            handle_error("Es necesario indicar parámetro levelist cuando levtype es " + parameters["levtype"][0], "levelist", parameters["levtype"][0])
+        for level in parameters["levelist"]:
+            if level not in constants.levelist[parameters["levtype"][0]]:
+                handle_error("Es valor " + level + " no es válido para el parámetro levelist, dado levtype = " + parameters["levtype"][0], "levelist", parameters["levtype"][0])
+                
+    # Parsear y validar el parametro number en el caso perturbed forecast
+    if parameters["type"][0] == "pf":
+        if "number" not in parameters:
+            handle_error("El parámetro number es obligatorio en el caso en que type es pf (Perturbed forecast).", "number")
+        else:
+            if parameters["number"][0] == "all":
+                parameters["number"] = [str(x) for x in range(1, 51)]
+            else:
+                try:
+                    numbers = [int(x) for x in parameters["number"]]
+                    for number in numbers:
+                        if number < 1 or number > 50:
+                            handle_error("El valor ingresado para el parámetro number es inválido", "number")
+                except ValueError:
+                    handle_error("El valor ingresado para el parámetro number es inválido", "number")
+    
+    # Se crea un request por cada conjunto de valores origin_time compatibles en el tiempo
+    # (No sé si es necesario, pero es lo que hace la interfaz web de TIGGE)
+    requests = {}
+    for curr_origin_time in parameters["origin_time"]:
+        curr_path = requests
+        for time in curr_origin_time[1:]:
+            if time not in curr_path:
+                curr_path[time] = {}
+            curr_path = curr_path[time]
+        if "origin" not in curr_path:
+            curr_path["origin"] = []
+        curr_path["origin"].append(curr_origin_time[0])
+    
+    # Esqueleto del request
+    api_parameters = {
+        "class": "ti",
+        "dataset": "tigge",
+        "expect": "any",
+        "expver": "prod",
+        "grid": "1.0/1.0",
+        "area":"-28/300/-36/309"
+    }
+    
+    # Cargar parte de los parámetros del request
+    api_parameters["date"] = parameters["start_date"][0] + "/to/" + parameters["end_date"][0]
+    api_parameters["levtype"] = parameters["levtype"][0]
+    insertion_set = ["levtype", "param", "step", "type"]
+    if parameters["type"][0] == "pf":
+        insertion_set.append("number")
+    if parameters["levtype"][0] != "sfc":
+        insertion_set.append("levelist")
+    for par in optional_parameters:
+        if par in parameters:
+            insertion_set.append(par)
+    for par in insertion_set:
+        value = "/".join(parameters[par])
+        api_parameters[par] = value
+    
+    remaining_requests = [requests]
+    remaining_times = [""]
+    while len(remaining_requests) > 0:
+        curr_request = remaining_requests.pop(0)
+        curr_time = remaining_times.pop(0)
+        for time in constants.time:
+            if time in curr_request:
+                remaining_requests.append(curr_request[time])
+                remaining_times.append(curr_time + "/" + time + ":00:00")
+        if curr_time and "origin" in curr_request:
+            origin = ""
+            origin_target = ""
+            for origin_id in curr_request["origin"]:
+                origin = origin + "/" + constants.o_map[int(origin_id)][1]
+                origin_target = origin_target + "-" + constants.o_map[int(origin_id)][0]
+            # Terminar de cargar los parámetros del request
+            api_parameters["origin"] = origin[1:]
+            api_parameters["time"] = curr_time[1:]
+            api_parameters["target"] = origin_target[1:] + "_" + api_parameters["type"] + "_" + api_parameters["date"].replace("/", "-") + ".grib"
+            # Hacer request
+            print(api_parameters)
+            download(api_parameters)