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
 21import csv
 22import os
 23from datetime import datetime, time, timedelta
 24from pathlib import Path
 25
 26import numpy as np
 27import pandas as pd
 28
 29# debug
 30# try:
 31#     aqui = Path(__file__).parent
 32# except:
 33#     aqui = Path()
 34# Ruta a los datos de estaciones
 35aqui = Path(__file__).parent
 36"""@private"""
 37ruta_data = aqui / "DB_DMC"
 38"""@private"""
 39
 40# TODO test this:
 41# assert len(list(ruta_data.glob("*.csv"))) > 1
 42# assert (ruta_data / "Estaciones.csv").is_file()
 43
 44
 45def transform_weather_file(filename):
 46    """Flip the direction of the wind in a weather file. Column containing wind direction must be called "WD".\
 47        this changes the contents of the given file."""
 48    temp_name = f"{filename}_tmp"
 49    with open(filename, newline='') as csvfile, open(temp_name, mode='w') as outfile:
 50        reader = csv.DictReader(csvfile)
 51        writer = csv.DictWriter(outfile, reader.fieldnames)
 52        writer.writeheader()
 53        for i in reader:
 54            i['WD'] = flip_wind(float(i["WD"]))
 55            writer.writerow(i)
 56    if os.path.isfile(filename):
 57        os.remove(filename)
 58        os.rename(temp_name, filename)
 59
 60
 61def is_qgis_running():
 62    qgis = False
 63    try:
 64        from qgis.core import QgsApplication  # type: ignore[import]
 65
 66        if QgsApplication.instance():  # qgis is running
 67            qgis = True
 68    except ModuleNotFoundError:
 69        qgis = False
 70    return qgis
 71
 72
 73qgis = is_qgis_running()
 74"""@private"""
 75
 76
 77def qgis_distance(fila, lat, lon):
 78    """Qgis distance calculation, used when QGIS is already running..."""
 79    if qgis:
 80        from qgis.core import QgsDistanceArea, QgsPointXY, QgsUnitTypes  # type: ignore[import] # , QgsCoordinateReferenceSystem, QgsProject
 81        print("Import QgsDistance")
 82    dist = QgsDistanceArea()
 83    dist.setEllipsoid("WGS84")
 84    # dist.setSourceCrs(QgsCoordinateReferenceSystem("EPSG:4326"), QgsProject.instance().transformContext())
 85    assert dist.lengthUnits() == QgsUnitTypes.DistanceMeters
 86    p1 = QgsPointXY(fila["lon"], fila["lat"])
 87    p2 = QgsPointXY(lon, lat)
 88    d = dist.measureLine(p1, p2)
 89    # return d
 90    return dist.convertLengthMeasurement(d, QgsUnitTypes.DistanceDegrees)
 91
 92
 93def eucl_distance(fila, lat, lon):
 94    """Euclidean distance calculation, used when QGIS is not running"""
 95    if lat == fila["lat"] and lon == fila["lon"]:
 96        return 0
 97    return np.sqrt((fila["lat"] - lat) ** 2 + (fila["lon"] - lon) ** 2)
 98
 99
100def distance(fila, lat, lon):
101    """Select distance calculation method whether QGIS is running or not"""
102    if qgis:
103        return qgis_distance(fila, lat, lon)
104    else:
105        return eucl_distance(fila, lat, lon)
106
107
108def flip_wind(a):
109    """Leeward to Windward wind angle direction flip. Barlovento a sotavento. Downwind to upwind."""
110    return round((a + 180) % 360, 2)
111
112
113def generate(
114    lat=-36.0,
115    lon=-73.2,
116    start_datetime=datetime(1989, 1, 12, 12, 0, 0),
117    rowres=60,
118    numrows=12,
119    numsims=100,
120    percentile=0.5,
121    outdir=Path("weather"),
122    dn=3,
123):
124    """Carolina Lorem Ipsum Dolor Sit Amet Consectetur Adipiscing Elit
125
126    Args:
127
128        lat (float): latitude-coordinate of the ignition point, EPSG 4326
129        lon (float): longitude-coordinate of the ignition point, EPSG 4326
130        start_datetime (datetime): starting time of the weather scenario
131        rowres (int): time resolution in minutes (not implemented yet)
132        numrows (int): number of hours in the weather scenario
133        numsims (int): number of weather scenarios
134        percentile (float): daily maximum temperature quantil
135        outdir (Path): output directory
136        dn (int): number of closest stations to base the scenario on
137
138    Return:
139
140        retval (int): 0 if successful, 1 otherwise, 2...
141        outdict (dict): output dictionary at least 'filelist': list of filenames created
142    """
143
144    filelist = []
145    try:
146        if not outdir.is_dir():
147            outdir.mkdir()
148            # print(f"Creating directory {outdir.absolute()}")
149
150        # leer estaciones
151        list_stn = pd.read_csv(ruta_data / "Estaciones.csv")
152        # calcular distancia a input point
153        list_stn["Distancia"] = list_stn.apply(distance, args=(lat, lon), axis=1)
154
155        #
156        # Vincenty’s formula >> euclidean distance
157        #
158        # list_stn["d1"] = list_stn.apply(distancia, args=(lat, lon), axis=1)
159        # list_stn["d2"] = list_stn.apply(distancia2, args=(lat, lon), axis=1)
160        # list_stn["dd"] = list_stn["d1"] - list_stn["d2"]
161        # list_stn[["dd","d1","d2"]]
162        # dd = list_stn[["dd","d1","d2"]].sort_values(by=["d1"], ascending=True)
163        # assert dd['d2'].is_monotonic_increasing
164        #           dd        d1        d2
165        # 0   0.006414  0.791903  0.785489
166        # 2   0.196967  1.302510  1.105544
167        # 1   0.094724  1.725775  1.631051
168        # 4   0.031122  1.868370  1.837248 <- !
169        # 3   0.305407  2.005146  1.699739
170        # 5   0.062915  2.387851  2.324937
171        # 6   0.029269  2.825473  2.796204 <- !
172        # 12  0.155427  2.830431  2.675004
173
174        # get 3 closest stations
175        stn = list_stn.sort_values(by=["Distancia"]).head(dn)["nombre"].tolist()
176
177        meteos = pd.DataFrame()
178        for st in stn:
179            # df = pd.read_csv(ruta_data / f"{st}.csv", sep=",", index_col=0, parse_dates=True)
180            # https://stackoverflow.com/questions/29206612/difference-between-data-type-datetime64ns-and-m8ns
181            df = pd.read_csv(ruta_data / f"{st}.csv", sep=",", index_col=0)
182            df.index = pd.to_datetime(df.index, errors="coerce")
183            df["station"] = st
184            # serie temperatura diaria
185            ser_tmp = df["TMP"].resample("D").max()
186            qn_date = ser_tmp[ser_tmp >= ser_tmp.quantile(percentile)].index
187            meteos = pd.concat([meteos, df[df.index.floor("D").isin(qn_date)].reset_index()], ignore_index=True)  # ?
188        # meteos["datetime"] = pd.to_datetime(meteos["datetime"], errors="coerce") # ?
189        assert meteos["datetime"].dtype == "datetime64[ns]"
190        # available days by stations
191        days = meteos.groupby(meteos.datetime.dt.date).first()["station"]
192
193        for i in range(numsims):
194            # draw station and day
195            cd = 0
196            ch = 0
197            while True:
198                station = np.random.choice(stn)
199                # TODO mejora ocupar estaciones sorteadas, en vez de posiblemente repetir
200                # station = np.random.choice(stn, size=len(stn), replace=False)
201                chosen_days = days[days == station]
202                if chosen_days.empty:
203                    if cd > dn * 10:  # 10 veces el numero de estaciones asegura que todas fueron sorteadas
204                        # print("Not enough data days", cd, ch)
205                        # break
206                        return 1, {"filelist": [], "exception": "No [enough] data in closest stations"}
207                    cd += 1
208                    continue
209                day = np.random.choice(chosen_days.index)
210                # retroceder un dia cada vez que falta
211                start = datetime.combine(day - timedelta(days=ch), start_datetime.time())
212                chosen_meteo = meteos[(meteos["datetime"] >= start) & (meteos["station"] == station)]
213                if len(chosen_meteo) < numrows:
214                    if ch > len(meteos):
215                        # print("Not enough data hours", cd, ch)
216                        # break
217                        return 1, {"filelist": [], "exception": "Not enough data"}
218                    ch += 1
219                    continue
220                break
221            # take rows
222            chosen_meteo = chosen_meteo.head(numrows)
223            # print(f"Selected {chosen_meteo.shape} from {station} on {day}")
224            # drop station
225            # TODO no drop ?
226            # chosen_meteo = chosen_meteo.drop(columns=["station"])
227            # scenario name
228            chosen_meteo.loc[:, "Scenario"] = "DMC" if numsims == 1 else f"DMC_{i+1}"
229            # TODO sobra: datetime format
230            # chosen_meteo.loc[:, "datetime"] = chosen_meteo["datetime"].dt.strftime("%Y-%m-%dT%H:%M:%S")
231            # TODO no aplanar al mismo dia
232            # chosen_meteo.loc[:, "datetime"] = [
233            #     chosen_meteo["datetime"].iloc[0] + timedelta(hours=i) for i in range(numrows)
234            # ]
235            # reorder
236            # chosen_meteo = chosen_meteo[["Scenario", "datetime", "WS", "WD", "TMP", "RH"]]
237            all_cols = chosen_meteo.columns.tolist()
238            first_cols = ["Scenario", "datetime", "WS", "WD", "TMP", "RH"]
239            for col in first_cols:
240                all_cols.remove(col)
241            chosen_meteo = chosen_meteo[first_cols + all_cols]
242            # write
243            tmpfile = outdir / ("Weather.csv" if numsims == 1 else f"Weather{i + 1}.csv")
244            filelist += [tmpfile.name]
245            chosen_meteo.to_csv(tmpfile, header=True, index=False)
246            # print(f"Writing {tmpfile.name} with {len(chosen_meteo)} rows")
247
248        return 0, {"filelist": filelist}
249    except Exception as e:
250        return 1, {"filelist": filelist, "exception": e}
251
252
253if __name__ == "__main__":
254    return_code, return_dict = generate()
255    if return_code == 0:
256        filelist = return_dict["filelist"]
257        print(f"Generated {len(filelist)} files: {filelist[0]}..{filelist[-1]}")
258    else:
259        print(f"Error generating files: {return_dict['exception']}")
260
261    sys.exit(return_code)