fire2a.meteo

  1#!python3
  2"""<b>Meteo Kitral</b> generates weather scenarios files for Cell2Fire-W Simulator using the Kitral fuel model standard.<br>
  3Using real Chilean weather station data from the Valparaíso to the Araucanía region in summer; Defining the target area (32S to 40S)<br>
  4<br>
  5Usage:<br>
  6- Selecting a point layer as <b>location</b> will pick the three nearest weather stations for sampling.<br>
  7- <b>Start hour</b>: Time of day from where to start picking station data.<br>
  8- <b>Temperature quantile</b>: A number greater than or equal to 0 and less than 1. It takes the daily maximum temperature values that are above the desired proportion. Example: quantile 0.75 takes the days when the daily maximum temperature is above the 75 &#37; <br>
  9- <b>Length of each scenario </b>: Indicates the duration, in hours, of each scenario<br>
 10- <b>Number_of_simulations</b>: files to generate<br>
 11- <b>output_directory</b>: where the files are written containing Weather(+digit).csv numbered files with each weather scenario<br>
 12<br>
 13Future Roadmap:<br>
 14- <b>step resolution</b>: Do other than hourly weather scenarios, to be used with the --Weather-Period-Length option (that defaults to 60)<br>
 15- Draw an animated vector layer representing the weather scenarios as arrows<br>
 16"""
 17__author__ = "Caro"
 18__revision__ = "$Format:%H$"
 19
 20import sys
 21from datetime import datetime, time, timedelta
 22from pathlib import Path
 23
 24import numpy as np
 25import pandas as pd
 26
 27# debug
 28try:
 29    aqui = Path(__file__).parent
 30except:
 31    aqui = Path()
 32# Ruta a los datos de estaciones
 33ruta_data = aqui / "DB_DMC"
 34assert len(list(ruta_data.glob("*.csv"))) > 1
 35assert (ruta_data / "Estaciones.csv").is_file()
 36
 37
 38def file_name(i, numsims):
 39    if numsims > 1:
 40        return f"Weather{i+1}.csv"
 41    return "Weather.csv"
 42
 43
 44def scenario_name(i, numsims):
 45    if numsims > 1:
 46        return f"DMC_{i+1}"
 47    return "DMC"
 48
 49
 50ok = False
 51try:
 52    from qgis.core import QgsApplication  # type: ignore[import]
 53    from qgis.core import QgsDistanceArea, QgsPointXY, QgsUnitTypes  # , QgsCoordinateReferenceSystem, QgsProject
 54
 55    if QgsApplication.instance():  # qgis is running
 56
 57        def distancia(fila, lat, lon):
 58            dist = QgsDistanceArea()
 59            dist.setEllipsoid("WGS84")
 60            # dist.setSourceCrs(QgsCoordinateReferenceSystem("EPSG:4326"), QgsProject.instance().transformContext())
 61            assert dist.lengthUnits() == QgsUnitTypes.DistanceMeters
 62            p1 = QgsPointXY(fila["lon"], fila["lat"])
 63            p2 = QgsPointXY(lon, lat)
 64            d = dist.measureLine(p1, p2)
 65            # return d
 66            return dist.convertLengthMeasurement(d, QgsUnitTypes.DistanceDegrees)
 67
 68        ok = True
 69except ModuleNotFoundError:
 70    ok = False
 71
 72if not ok:
 73
 74    def distancia(fila, lat, lon):
 75        if lat == fila["lat"] and lon == fila["lon"]:
 76            return 0
 77        return np.sqrt((fila["lat"] - lat) ** 2 + (fila["lon"] - lon) ** 2)
 78
 79
 80def meteo_to_c2f(alfa):
 81    if alfa >= 0 and alfa < 180:
 82        return round(alfa + 180, 2)
 83    elif alfa >= 180 and alfa <= 360:
 84        return round(alfa - 180, 2)
 85    return np.nan
 86
 87
 88def barlo_sota(a):
 89    return round((a + 180) % 360, 2)
 90
 91
 92def generate(
 93    lat=-36.0,
 94    lon=-73.2,
 95    start_datetime=datetime(1989, 1, 12, 12, 0, 0),
 96    rowres=60,
 97    numrows=12,
 98    numsims=100,
 99    percentile=0.5,
100    outdir=Path("weather"),
101    dn=3,
102):
103    """Carolina Lorem Ipsum Dolor Sit Amet Consectetur Adipiscing Elit
104    Args:
105        lat (float): latitude-coordinate of the ignition point, EPSG 4326
106        lon (float): longitude-coordinate of the ignition point, EPSG 4326
107        start_datetime (datetime): starting time of the weather scenario
108        rowres (int): time resolution in minutes (not implemented yet)
109        numrows (int): number of hours in the weather scenario
110        numsims (int): number of weather scenarios
111        percentile (float): daily maximum temperature quantil
112        outdir (Path): output directory
113        dn (int): number of closest stations to base the scenario on
114    Return:
115        retval (int): 0 if successful, 1 otherwise, 2...
116        outdict (dict): output dictionary at least 'filelist': list of filenames created
117    """
118
119    filelist = []
120    try:
121        if not outdir.is_dir():
122            outdir.mkdir()
123            # print(f"Creating directory {outdir.absolute()}")
124
125        # leer estaciones
126        list_stn = pd.read_csv(ruta_data / "Estaciones.csv")
127        # calcular distancia a input point
128        list_stn["Distancia"] = list_stn.apply(distancia, args=(lat, lon), axis=1)
129
130        #
131        # Vincenty’s formula >> euclidean distance
132        #
133        # list_stn["d1"] = list_stn.apply(distancia, args=(lat, lon), axis=1)
134        # list_stn["d2"] = list_stn.apply(distancia2, args=(lat, lon), axis=1)
135        # list_stn["dd"] = list_stn["d1"] - list_stn["d2"]
136        # list_stn[["dd","d1","d2"]]
137        # dd = list_stn[["dd","d1","d2"]].sort_values(by=["d1"], ascending=True)
138        # assert dd['d2'].is_monotonic_increasing
139        #           dd        d1        d2
140        # 0   0.006414  0.791903  0.785489
141        # 2   0.196967  1.302510  1.105544
142        # 1   0.094724  1.725775  1.631051
143        # 4   0.031122  1.868370  1.837248 <- !
144        # 3   0.305407  2.005146  1.699739
145        # 5   0.062915  2.387851  2.324937
146        # 6   0.029269  2.825473  2.796204 <- !
147        # 12  0.155427  2.830431  2.675004
148
149        # get 3 closest stations
150        stn = list_stn.sort_values(by=["Distancia"]).head(dn)["nombre"].tolist()
151
152        meteos = pd.DataFrame()
153        for st in stn:
154            # df = pd.read_csv(ruta_data / f"{st}.csv", sep=",", index_col=0, parse_dates=True)
155            # https://stackoverflow.com/questions/29206612/difference-between-data-type-datetime64ns-and-m8ns
156            df = pd.read_csv(ruta_data / f"{st}.csv", sep=",", index_col=0)
157            df.index = pd.to_datetime(df.index, errors="coerce")
158            df["station"] = st
159            # serie temperatura diaria
160            ser_tmp = df["TMP"].resample("D").max()
161            qn_date = ser_tmp[ser_tmp >= ser_tmp.quantile(percentile)].index
162            meteos = pd.concat([meteos, df[df.index.floor("D").isin(qn_date)].reset_index()], ignore_index=True)  # ?
163        # meteos["datetime"] = pd.to_datetime(meteos["datetime"], errors="coerce") # ?
164        assert meteos["datetime"].dtype == "datetime64[ns]"
165        # available days by stations
166        days = meteos.groupby(meteos.datetime.dt.date).first()["station"]
167
168        for i in range(numsims):
169            # draw station and day
170            cd = 0
171            ch = 0
172            while True:
173                station = np.random.choice(stn)
174                # TODO mejora ocupar estaciones sorteadas, en vez de posiblemente repetir
175                # station = np.random.choice(stn, size=len(stn), replace=False)
176                chosen_days = days[days == station]
177                if chosen_days.empty:
178                    if cd > dn * 10:  # 10 veces el numero de estaciones asegura que todas fueron sorteadas
179                        # print("Not enough data days", cd, ch)
180                        # break
181                        return 1, {"filelist": [], "exception": "No [enough] data in closest stations"}
182                    cd += 1
183                    continue
184                day = np.random.choice(chosen_days.index)
185                # retroceder un dia cada vez que falta
186                start = datetime.combine(day - timedelta(days=ch), start_datetime.time())  
187                chosen_meteo = meteos[(meteos["datetime"] >= start) & (meteos["station"] == station)]
188                if len(chosen_meteo) < numrows:
189                    if ch > len(meteos):
190                        # print("Not enough data hours", cd, ch)
191                        # break
192                        return 1, {"filelist": [], "exception": "Not enough data"}
193                    ch += 1
194                    continue
195                break
196            # take rows
197            chosen_meteo = chosen_meteo.head(numrows)
198            # print(f"Selected {chosen_meteo.shape} from {station} on {day}")
199            # drop station
200            # TODO no drop ?
201            # chosen_meteo = chosen_meteo.drop(columns=["station"])
202            # wind direction
203            chosen_meteo.loc[:, "WD"] = chosen_meteo["WD"].apply(barlo_sota)
204            # scenario name
205            chosen_meteo.loc[:, "Scenario"] = scenario_name(i, numsims)
206            # TODO sobra: datetime format
207            # chosen_meteo.loc[:, "datetime"] = chosen_meteo["datetime"].dt.strftime("%Y-%m-%dT%H:%M:%S")
208            # TODO no aplanar al mismo dia
209            # chosen_meteo.loc[:, "datetime"] = [
210            #     chosen_meteo["datetime"].iloc[0] + timedelta(hours=i) for i in range(numrows)
211            # ]
212            # reorder
213            # chosen_meteo = chosen_meteo[["Scenario", "datetime", "WS", "WD", "TMP", "RH"]]
214            all_cols = chosen_meteo.columns.tolist()
215            first_cols = ["Scenario", "datetime", "WS", "WD", "TMP", "RH"]
216            for col in first_cols:
217                all_cols.remove(col)
218            chosen_meteo = chosen_meteo[first_cols + all_cols]
219            # write
220            tmpfile = outdir / file_name(i, numsims)
221            filelist += [tmpfile.name]
222            chosen_meteo.to_csv(tmpfile, header=True, index=False)
223            # print(f"Writing {tmpfile.name} with {len(chosen_meteo)} rows")
224
225        return 0, {"filelist": filelist}
226    except Exception as e:
227        return 1, {"filelist": filelist, "exception": e}
228
229
230if __name__ == "__main__":
231    return_code, return_dict = generate()
232    if return_code == 0:
233        filelist = return_dict["filelist"]
234        print(f"Generated {len(filelist)} files: {filelist[0]}..{filelist[-1]}")
235    else:
236        print(f"Error generating files: {return_dict['exception']}")
237
238    sys.exit(return_code)