diff --git a/.gitignore b/.gitignore index 77c12fc8..bc7e47f8 100644 --- a/.gitignore +++ b/.gitignore @@ -1,25 +1,70 @@ -#Ignoring files and folders by github +slurm* +output/ +#dir +src/__pycache__/ +output/* +old/ +*.png +*.npy +#gvim +*.bak +*.swp +*.swo +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class -#build directories -bin -build -lib -sconsbuild -forefire2 +# C extensions +*.so -#python -__pycache__ +# Distribution / packaging +.Python +env/ +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +*.egg-info/ +.installed.cfg +*.egg -# compilation files -*.dblite -*.o +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec -# editors -.vscode +# Installer logs +pip-log.txt +pip-delete-this-directory.txt -#mac files -**/.DS_Store +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*,cover -# jupyter notebooks checkpoints -.ipynb_checkpoints -CMakeLists.txt +# Translations +*.mo +*.pot + +# Django stuff: +*.log + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ diff --git a/cmake-build.sh b/cmake-build.sh old mode 100644 new mode 100755 diff --git a/src/Andrews.cpp b/src/Andrews.cpp new file mode 100644 index 00000000..dbe456fd --- /dev/null +++ b/src/Andrews.cpp @@ -0,0 +1,214 @@ +/* + +Copyright (C) 2012 ForeFire Team, SPE, Universit� de Corse. + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU Lesser General Public +License as published by the Free Software Foundation; either +version 2.1 of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public +License along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 US + + */ + + + +#include "PropagationModel.h" +#include "FireDomain.h" +#include +using namespace std; +namespace libforefire { + +class Andrews: public PropagationModel { + + /*! name the model */ + static const string name; + + /*! boolean for initialization */ + static int isInitialized; + double windReductionFactor; + /*! properties needed by the model */ + + size_t slope; + size_t Md; + size_t normalWind; + size_t Sigmad; + size_t sd; + size_t Rhod; + size_t e; + size_t Mchi; + size_t DeltaH; + + /*! coefficients needed by the model */ + + /*! local variables */ + + /*! result of the model */ + double getSpeed(double*); + +public: + Andrews(const int& = 0, DataBroker* db=0); + virtual ~Andrews(); + + string getName(); + +}; + +PropagationModel* getAndrewsModel(const int& = 0, DataBroker* db=0); + + +/* name of the model */ +const string Andrews::name = "Andrews"; + +/* instantiation */ +PropagationModel* getAndrewsModel(const int & mindex, DataBroker* db) { + return new Andrews(mindex, db); +} + +/* registration */ +int Andrews::isInitialized = + FireDomain::registerPropagationModelInstantiator(name, getAndrewsModel ); + +/* constructor */ +Andrews::Andrews(const int & mindex, DataBroker* db) +: PropagationModel(mindex, db) { + /* defining the properties needed for the model */ + windReductionFactor = params->getDouble("windReductionFactor"); + + slope = registerProperty("slope"); + Md = registerProperty("moisture"); + + normalWind = registerProperty("normalWind"); + Sigmad = registerProperty("fuel.Sigmad"); + sd = registerProperty("fuel.sd"); + Rhod = registerProperty("fuel.Rhod"); + e = registerProperty("fuel.e"); + Mchi = registerProperty("fuel.Mchi"); + DeltaH = registerProperty("fuel.DeltaH"); + + /* allocating the vector for the values of these properties */ + if ( numProperties > 0 ) properties = new double[numProperties]; + + /* registering the model in the data broker */ + dataBroker->registerPropagationModel(this); + + /* Definition of the coefficients */ +} + +/* destructor (shoudn't be modified) */ +Andrews::~Andrews() { +} + +/* accessor to the name of the model */ +string Andrews::getName(){ + return name; +} + +/* *********************************************** */ +/* Model for the propagation velovity of the front */ +/* *********************************************** */ + +double Andrews::getSpeed(double* valueOf){ + + double lRhod = valueOf[Rhod] * 0.06; // conversion kg/m^3 -> lb/ft^3 + double lMd = valueOf[Md]; + double lMchi = valueOf[Mchi]; + double lsd = valueOf[sd] / 3.2808399; // conversion 1/m -> 1/ft + double le = valueOf[e] * 3.2808399; // conversion m -> ft + if (le==0) return 0; + double lSigmad = valueOf[Sigmad] * 0.2048; // conversion kg/m^2 -> lb/ft^2 + double lDeltaH = valueOf[DeltaH] / 2326.0;// conversion J/kg -> BTU/lb + double normal_wind = valueOf[normalWind] * 196.850394 ; //conversion m/s -> ft/min + double localngle = valueOf[slope]; + // if (normal_wind > 0) cout<<"wind is"<Uf) { + normal_wind = Uf; + } + + double phiV = C * pow((Beta/Betaop), -E) * pow(normal_wind,B) ; + + double phiP = 5.275 * pow(Beta, -0.3) * pow(tanangle, 2); + + double R0 = (Ir * chi) / (lRhobulk * epsilon * Qig); + + double R = R0 * (1 + phiV + phiP); + + if(R < R0) R = R0; + + if(R > 0.0) { + //cout << "Andrews =" << R * 0.00508 << endl; + return R * 0.005588 ; // ft/min -> m/s + }else{ + //cout << " Rhod "<< lRhod << " lMd "<< lMd << " lsd "<< lsd << " le "<< le << " lSigmad "<< lSigmad <fd->getDomainID()<<" is putting "<fd->getDomainID()<<" is putting "<* myLayer = session->fd->getDataLayer(tmpname); - if ( myLayer ){ FFArray* myMatrix; // getting the pointer @@ -247,7 +246,6 @@ void FFPutDoubleArray(const char* mname, double* x, myMatrix->copyDataToFortran(x); } else { - cout<<"Error trying to put data from unknown layer "<* layer) { /* inserting the layer into the map of layers */ - - ilayer = layersMap.find(name); if (ilayer != layersMap.end()) { @@ -292,6 +290,7 @@ void DataBroker::registerLayer(string name, DataLayer* layer) { } } + layersMap.insert(make_pair(name, layer)); layers.push_back(layer); @@ -318,7 +317,6 @@ void DataBroker::registerLayer(string name, DataLayer* layer) { temperatureLayer = layer; } if (name.find("windU") != string::npos){ - windULayer = layer; } if (name.find("windV") != string::npos){ @@ -912,7 +910,6 @@ void DataBroker::loadFromNCFile(string filename) { NcVarAtt layerTypeNC; string layerType; for (auto it = allVariables.begin(); it != allVariables.end(); ++it){ - try { layerTypeNC = it->second.getAtt("type"); if (layerTypeNC.isNull()) { @@ -936,8 +933,6 @@ void DataBroker::loadFromNCFile(string filename) { string varName = it->first; - - if (varName == "wind") { XYZTDataLayer * wul = constructXYZTLayerFromFile(&dataFile, "wind",0); registerLayer("windU", wul); @@ -1073,7 +1068,7 @@ FFPoint DataBroker::getNetCDFSWCorner(NcVar* var) { params->getDouble("SHIFT_ALL_DATA_ORDINATES_BY")); return shiftPos + relSWC; } -double DataBroker::getNetCDFTimeOrigin(NcVar* var) { +double DataBroker::getNetCDFTimeSpan(NcVar* var) { double Lt =0; NcVarAtt att = var->getAtt("Lt"); if(!att.isNull()) att.getValues(&Lt); @@ -1091,7 +1086,7 @@ FFPoint DataBroker::getNetCDFSpatialSpan(NcVar* var) { if(!att.isNull()) att.getValues(&Lz); return FFPoint(Lx, Ly, Lz); } -double DataBroker::getNetCDFTimeSpan(NcVar* var) { +double DataBroker::getNetCDFTimeOrigin(NcVar* var) { double timeOrigin =0; NcVarAtt att = var->getAtt("t0"); if(!att.isNull()) att.getValues(&timeOrigin); @@ -1143,9 +1138,9 @@ XYZTDataLayer* DataBroker::constructXYZTLayerFromFile( NcFile* NcdataFil if (isRelevantData(SWCorner, spatialExtent)) { // - //cout << "Variable " << property <<" " <* newlayer = new XYZTDataLayer(property, SWCorner, timeOrigin, spatialExtent, Lt, nx, ny, nz, nt, data); + XYZTDataLayer* newlayer = new XYZTDataLayer(property, SWCorner, timeOrigin, spatialExtent, Lt, nx, ny, nz, nt, data); delete[] data; return newlayer; diff --git a/src/FireDomain.cpp b/src/FireDomain.cpp index dbaa4e11..5f599199 100644 --- a/src/FireDomain.cpp +++ b/src/FireDomain.cpp @@ -805,7 +805,7 @@ void FireDomain::readMultiDomainMetadata(){ XYZTDataLayer* newLayer = new XYZTDataLayer(name, origin,t0, span, timespan, nnx, nny, nnz, nnk, values); dataBroker->registerLayer(name, newLayer); - //cout<<"adding scalar Layer "<= 0) -cout << "formObs " << et << ' ' << bt << ' ' << at << ' ' << " - " - << burningTime/500 << '|' << residenceTime + + if (at >= 0) + cout << "formObs " << et << ' ' << bt << ' ' << at << ' ' << " - " + << burningTime/500 << '|' << residenceTime << ' ' << nominalHeatFlux_f/1.8e6 << '|' < 0. ) { + if ( overspeed > 0. ) { return R0 + overspeed; } else { return R0; diff --git a/tools/postprocessing/pMNHFF2VTK.py b/tools/postprocessing/pMNHFF2VTK.py index f6d06c80..57b4e01a 100755 --- a/tools/postprocessing/pMNHFF2VTK.py +++ b/tools/postprocessing/pMNHFF2VTK.py @@ -14,6 +14,9 @@ from copy import copy import glob import netCDF4 as nc4 +import pdb + + class LaserShot(): valueName = "PR2" @@ -809,7 +812,7 @@ def ffFrontsToVtk(inFFpattern = "", outPath = ""): def ffmnhFileToVtk(inpattern="", pgdFile="", outPath="", cleanFile=False, lidarIn=None, lidarOut=None, - startStep=-1, endStep=-1, genDomainOrigin=None, genDomainExtent=None, norecompute=False, + startStep=-1, endStep=-1, dtStep = 1, genDomainOrigin=None, genDomainExtent=None, norecompute=False, quitAfterCompute=False, xcfdName=None, inputcdfvars=(), vect_vars=("U", "V", "W")): # Obtenir le préfixe du fichier @@ -892,6 +895,7 @@ def ffmnhFileToVtk(inpattern="", pgdFile="", outPath="", cleanFile=False, lidarI outname = "%s/%s.full.%d.vts"%(outPath,fprefix,stepV) if (cleanFile and os.path.isfile(outname) ): print("Step %d already post-processed, cleaning up"%stepV) + pdb.set_trace() delCmd = "rm %s.*.%d"%(inFFpattern,stepV) print(delCmd) #os.system(delCmd) @@ -976,7 +980,7 @@ def ffmnhFileToVtk(inpattern="", pgdFile="", outPath="", cleanFile=False, lidarI if startStep > -1 and endStep > -1: - selectedSteps=selectedSteps[startStep:endStep] + selectedSteps=selectedSteps[startStep:endStep:dtStep] print(len(selectedSteps), " time steps to be computed on ", len(Allsteps) , " list :", selectedSteps, " ") @@ -1190,7 +1194,7 @@ def main(): parser.add_argument('-fields', nargs=3, help='Input : mnhdump file pattern. PGD nc file, VTK Output directory ') parser.add_argument('-clean', action='store_true', help='Enable this flag to clean the output files before processing.') parser.add_argument('-lidar', nargs=2, help='Lidar emulator input and output files. Provide two file paths, in and out csv.') - parser.add_argument('-steps', nargs=2, type=int, help='Start and end generation steps for processing. Provide two integers, start and end') + parser.add_argument('-steps', nargs=3, type=int, help='Start and end generation steps for processing. Provide two integers, start and end') parser.add_argument('-norecompute', action='store_true', help='Enable this flag to skip recompute phase but make the index vtk file') parser.add_argument('-quit', action='store_true', help='Enable this flag to quit after computation without further processing but display number of steps in data') parser.add_argument('-cdf', help='Provide a prefix for CDF output (one file per step). Example: "CDFOUT/step"') @@ -1212,6 +1216,8 @@ def main(): lidarOut=args.lidar[1] if args.lidar else None, startStep=args.steps[0] if args.steps else -1, endStep=args.steps[1] if args.steps else -1, + dtStep=args.steps[2] if args.steps else 1, + norecompute=args.norecompute, quitAfterCompute=args.quit, xcfdName=args.cdf diff --git a/tools/postprocessing/vtk2gpkg.py b/tools/postprocessing/vtk2gpkg.py new file mode 100644 index 00000000..e4b0e188 --- /dev/null +++ b/tools/postprocessing/vtk2gpkg.py @@ -0,0 +1,365 @@ +import numpy as np +import pyvista as pv +import geopandas as gpd +from shapely.geometry import Point, Polygon +import pandas as pd +import matplotlib.pyplot as plt +import sys +import cv2 +import os +import xarray as xr +from scipy import ndimage +import numpy as np +import xarray as xr +from scipy.ndimage import distance_transform_edt +import dask.array as da +from scipy.ndimage import gaussian_filter +import glob +from affine import Affine +from rasterio.features import rasterize +from scipy import interpolate +import rasterio +import xarray as xr +from rasterio.crs import CRS +import pyproj + +################################ +def fill_missing_nearest_2d(da): + """ + Fill NaN values in a 2D xarray DataArray with the value of the nearest valid cell (in Euclidean distance). + + Parameters + ---------- + da : xr.DataArray + A 2D DataArray (dims might be something like ("y", "x")) that contains NaNs. + + Returns + ------- + xr.DataArray + A new DataArray of the same shape and coordinates, but with NaN values filled + by the nearest non-NaN neighbor in Euclidean distance. + + Notes + ----- + 1. This approach uses scipy.ndimage.distance_transform_edt to compute, for each NaN cell, + the indices of the closest non-NaN cell. + 2. If your grid is very large, this can be memory-intensive. Consider chunking or a + KD-tree approach if needed. + 3. This method does not differentiate between “inside” or “outside” convex regions + of valid data; it simply finds the closest valid cell for every NaN location. + """ + if da.ndim != 2: + raise ValueError("This function is designed for 2D DataArrays only.") + + # Get raw numpy array and create a mask of where values are NaN + arr = da.values + mask_nan = np.isnan(arr) + + # distance_transform_edt returns: + # dist : array of distances + # idx : array of indices of nearest non-NaN cell + dist, idx = distance_transform_edt(mask_nan, return_indices=True) + + # Copy original data so we can fill in + arr_filled = arr.copy() + + # Fill in NaN values by indexing the nearest non-NaN cell + arr_filled[mask_nan] = arr_filled[idx[0, mask_nan], idx[1, mask_nan]] + + # Return a new xarray DataArray with the same dims/coords + return xr.DataArray(arr_filled, coords=da.coords, dims=da.dims) + + +simuName = os.getcwd().split('/')[-1] + +#load fartsite tif ros +##################### +# Open the .tif file with rasterio +dirfarsite = '/data/paugam/Data/ElPontDeVilamora/FARSITE-UPC/' +with rasterio.open(dirfarsite+'ros_farsite.tif') as src: + # Read the image data into a numpy array + data = src.read(1) # Assuming single-band data + + # Create an xarray DataArray + farsite = xr.DataArray(data, dims=["y", "x"], coords={"y": src.bounds.top - src.res[1] * np.arange(src.height), + "x": src.bounds.left + src.res[0] * np.arange(src.width)}) + +# Assign CRS using the EPSG code 28831 +farsite.attrs['crs'] = CRS.from_epsg(25831).to_string() +farsite = farsite.where(farsite != -9999, np.nan) + +non_nan_mask = farsite.notnull() + +# Get the indices of the non-NaN values (bounding box) +y_minP, y_maxP = farsite.coords["y"][non_nan_mask.any(dim="x")].min(), farsite.coords["y"][non_nan_mask.any(dim="x")].max() +x_minP, x_maxP = farsite.coords["x"][non_nan_mask.any(dim="y")].min(), farsite.coords["x"][non_nan_mask.any(dim="y")].max() +buffer = 1000 + + +#laod mnh output in lat lon +########################### +mnh = xr.open_dataset('../006_mnhSolo/Postproc/outputFile/Netcdf/FCAST_model3.nc', decode_times=False) +refdate = pd.Timestamp('2022-07-17 00:00:00') + + +# Load the vtu file and rasterize velocity +################################### +filesPVD = glob.glob("./vtkFront/*.vtu") +rosmap = farsite.copy() +rosmap.data[:,:] = -999 + +################################## +def get_interp(x,y,z,mask=None): + coord_pts = np.vstack((x.flatten(), y.flatten())).T + data = z.flatten() + fill_val = -999 + return interpolate.LinearNDInterpolator(coord_pts, data, fill_value=fill_val) + +y2d,x2d = np.meshgrid(mnh.y, mnh.x) +lon_interpolator = get_interp(x2d,y2d, mnh.lon.data) +lat_interpolator = get_interp(x2d,y2d, mnh.lat.data) +# Définir les systèmes de coordonnées : WGS84 (lat, lon) et UTM Zone 31N (EPSG:25831) +wgs84 = pyproj.CRS('EPSG:4326') # Système de référence pour latitude/longitude +utm31n = pyproj.CRS('EPSG:25831') # Système de projection UTM Zone 31N + +# Créer un transformateur pour effectuer la conversion entre les deux CRS +transformer = pyproj.Transformer.from_crs(wgs84, utm31n, always_xy=True) + +for iifile, file_ in enumerate(filesPVD): + print(file_) + front = pv.read(file_) + + points = front.points + vmod = front.point_data['vmod'] + velx = front.point_data['velocity'][:,0] + vely = front.point_data['velocity'][:,1] + velz = front.point_data['velocity'][:,2] + + data = { + 'geometry': [Point(p) for p in points], # Create Point geometries + 'vmod': vmod, # vmod data, + 'u': velx, # velocity data + 'v': vely, # velocity data + 'w': velz, # velocity data + } + + gdf = gpd.GeoDataFrame(data) # You can adjust the CRS if needed + + # Define the grid and extraction bounds from the MNH dataset + x_min, x_max = farsite.x.min().item(), farsite.x.max().item() + y_min, y_max = farsite.y.min().item(), farsite.y.max().item() + resolution_x = abs(farsite.x[1] - farsite.x[0]).item() + resolution_y = abs(farsite.y[1] - farsite.y[0]).item() + + # Create the transformation matrix for mapping coordinates to pixel locations + transform = Affine(resolution_x, 0.0, x_min, 0.0, -resolution_y, y_max) + + + # Create a list of shapes (geometries and associated value) + + shapes = [] + for geom,val in zip(gdf.geometry,gdf.vmod): + lon_interpolator(geom.x, geom.y) + x, y = transformer.transform(lon_interpolator(geom.x, geom.y), lat_interpolator(geom.x, geom.y)) + shapes.append( (Point(x,y), val ) ) + + + # Rasterize the shapes onto the grid defined by bmap_po + rasterized = rasterize( + shapes=shapes, + out_shape=farsite.shape[::-1], + transform=transform, + fill=-999, # Background value (no data value) + dtype=np.float32 + ) + + # Adjust the raster for any necessary flipping (this might depend on how your grid is oriented) + tmp = rasterized.T[::-1,::-1] # Flip vertically if needed + + # Update the bmap_po data where the raster has non-background values + rosmap.data = np.where( tmp > rosmap.data, tmp, rosmap.data ) + + +rosmap = rosmap.where(rosmap != -999, np.nan) * 60 # m/min +#rosmap = rosmap.assign_coords(lat=mnh.lat, lon=mnh.lon) +#rosmap = rosmap.drop_vars(['x', 'y']) # Drop the old x and y coordinates +#rosmap = rosmap.rename({'x': 'lon', 'y': 'lat'}) + +#set to 1d grid +#points = np.array([rosmap.lon.data.flatten(), rosmap.lat.data.flatten()]).T # (n_points, 2) +#values = rosmap.data.flatten() + +#new_lat_1d = np.linspace(np.min(mnh.lat), np.max(mnh.lat), mnh.lat.shape[0]) +#new_lon_1d = np.linspace(np.min(mnh.lon), np.max(mnh.lon), mnh.lat.shape[1]) +#new_lat_grid, new_lon_grid = np.meshgrid(new_lat_1d, new_lon_1d) + +# Perform the interpolation using griddata +#rosmap_interp = interpolate.griddata(points, values, (new_lon_grid,new_lat_grid), method='nearest') + +#rosmap = xr.DataArray( +# rosmap_interp.T, # The interpolated data +# coords={'lat': new_lat_1d, 'lon': new_lon_1d}, # 1D coordinates +# dims=['lat', 'lon'], # Dimension names +#) +#rosmap.rio.write_crs("EPSG:4326", inplace=True) +#rosmap = rosmap.rename({'lat':'y', 'lon':'x'}) +#rosmap = rosmap.rio.reproject(25831) + + + +ax = plt.subplot(121) +rosmap.plot(ax=ax, ) #vmax=15) +ax.set_title('ff') +ax.set_xlim(x_minP-buffer, x_maxP+buffer) +ax.set_ylim(y_minP-buffer, y_maxP+buffer) + +ax = plt.subplot(122) +farsite.plot(ax=ax, vmax=15) +ax.set_title('farsite') +ax.set_xlim(x_minP-buffer, x_maxP+buffer) +ax.set_ylim(y_minP-buffer, y_maxP+buffer) + + +plt.show() +sys.exit() + + + + +# Load the vts file +fileVTU = "./vtkmap3.vts" +mesh = pv.read(fileVTU) + +cell_centers = mesh.cell_centers() +x = cell_centers.points[:, 0] +y = cell_centers.points[:, 1] +bmap_cell = mesh.cell_data['atime'] + +# Get the unique x and y values to reshape data into a 2D grid +unique_x = np.unique(x) +unique_y = np.unique(y) + +# Create a grid (2D array) to map the cell data +# We use the unique x and y values to determine the grid shape +bmap2d = np.full((len(unique_y), len(unique_x)), np.nan) + +# Map the cell data onto the grid by matching cell centers to grid coordinates +for i, (xi, yi) in enumerate(zip(x, y)): + # Find the closest grid point in x and y + x_idx = np.abs(unique_x - xi).argmin() + y_idx = np.abs(unique_y - yi).argmin() + + # Set the cell data at the corresponding grid point + bmap2d[y_idx, x_idx] = bmap_cell[i] + +factorNiceFront = 4 +''' +unique_x2 = np.linspace(unique_x.min(),unique_x.max(),unique_x.shape[0]*factorNiceFront) +unique_y2 = np.linspace(unique_y.min(),unique_y.max(),unique_y.shape[0]*factorNiceFront) +bmap2d2 =ndimage.zoom(bmap2d, factorNiceFront, order=0) +bmap2d2[np.where(bmap2d2>bmap2d.max())] = bmap2d.max() +bmap2d2[np.where(bmap2d2<0)] = -9999 +sigma = 2 # The standard deviation of the Gaussian kernel +bmap2d2 = ndimage.gaussian_filter(bmap2d2, sigma=sigma) +''' + +dabmap = xr.DataArray(bmap2d, coords=[("y", unique_y), ("x", unique_x)]) + +new_x = np.linspace(dabmap.x.min(), dabmap.x.max(), dabmap.sizes['x'] * factorNiceFront) +new_y = np.linspace(dabmap.y.min(), dabmap.y.max(), dabmap.sizes['y'] * factorNiceFront) + +mnhlon = mnhrosmap.lon.interp( x=new_x, y=new_y, method='linear' ) +mnhlat = mnh.lat.interp( x=new_x, y=new_y, method='linear' ) + +# Interpolate to the new grid +dabmap_high_res = dabmap.where(dabmap>0).interp(x=new_x, y=new_y, method='linear') + +mask = da.isnan(dabmap_high_res) +dabmap_high_res = fill_missing_nearest_2d(dabmap_high_res) + +# Apply Gaussian smoothing +sigma = 1.5 # Standard deviation for Gaussian kernel +dabmap_smoothed = xr.DataArray( + gaussian_filter(dabmap_high_res.where(dabmap_high_res>0), sigma=sigma), + coords=dabmap_high_res.coords, + dims=dabmap_high_res.dims, + ) + +dabmap_smoothed = dabmap_smoothed.where(~mask, np.nan) +dabmap_smoothed = dabmap_smoothed.fillna(-9999) + +bmap2d = dabmap_smoothed.values +unique_x = dabmap_smoothed.x.values +unique_y = dabmap_smoothed.y.values + +times = np.unique(bmap2d) +gdf = None +dt = 600. + +timerange = np.arange(times[np.where(times>0)].min(), times[np.where(times>0)].max()+dt, dt) +for it, time in enumerate(timerange): + if time < 0 : continue + contours,hierarchy = cv2.findContours(np.where((bmap2d<=time)&(bmap2d>0),1,0).astype(np.uint8),cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE) + + #convert contours to lat lon + contours_ll = [] + for ic, contour in enumerate(contours): + if len(contour)<=3: continue + x_ = unique_x[contour[:,0,0]] + y_ = unique_y[contour[:,0,1]] + + lat_ = []; lon_ = [] + for ii,jj in zip(contour[:,0,0],contour[:,0,1]): + lat_.append(mnhlat[jj,ii]) + lon_.append(mnhlon[jj,ii]) + + contours_ll_ = np.array([lon_, lat_]).T.reshape(-1,1,2) + contours_ll.append(contours_ll_) + + #contours_ll_ = np.array([x_, y_]).T.reshape(-1,1,2) + #contours_ll.append(contours_ll_) + + # Create a list of polygons + # this codes shoyld only work for main polygon with holes indise. new polygon inside hole are not tackle + polygons = [] + for i, contour in enumerate(contours_ll): + if len(contour)<=3: continue + contour_points = contour.reshape(-1, 2) + + if hierarchy[0][i][3] == -1: # Has a hole + + holes = [] + # Find the child contours (holes) + for j, child_contour in enumerate(contours_ll): + if hierarchy[0][j][3] == i: + holes.append(child_contour.reshape(-1, 2)) + + if len(holes) == 0 : + polygon = Polygon(contour_points) + else: + polygon = Polygon(contour_points, holes=holes) + + polygons.append(polygon) + + # Create a GeoDataFrame from the list of polygons + if len(polygons) == 0 : + continue + gdf_ = gpd.GeoDataFrame({'geometry': polygons}) + gdf_ = gdf_.set_crs(4326) + #gdf_.geometry = gdf_.buffer(0) + gdf_['timestamp'] = refdate + pd.Timedelta(seconds=time) + gdf_['time_seconds'] = time + + if gdf is None: + gdf = gdf_ + else: + gdf = pd.concat([gdf, gdf_]).reset_index(drop=True) + + print( '{:.1f} | {:}'.format(100.*it/len(timerange), gdf_['timestamp'].iloc[0]) , end='\r') + +os.makedirs('./GPKG/', exist_ok=True) +if os.path.isfile("./GPKG/frontModel3_{:s}.gpkg".format(simuName)): + os.remove("./GPKG/frontModel3_{:s}.gpkg".format(simuName)) +gdf.to_file("./GPKG/frontModel3_{:s}.gpkg".format(simuName), driver="GPKG") + + diff --git a/tools/preprocessing/create_data.py b/tools/preprocessing/create_data.py new file mode 100644 index 00000000..73029f41 --- /dev/null +++ b/tools/preprocessing/create_data.py @@ -0,0 +1,412 @@ +import numpy as np +import matplotlib.pyplot as plt +import xarray as xr +import rioxarray +import pyforefire as forefire +import os +import sys +import argparse +import f90nml +import importlib +import cftime +from datetime import datetime +import geopandas as gpd +import glob +import shapely +import pdb +import pandas as pd + + +#homebrewed +import create_data_tools as tools + +################################################### +if __name__ == '__main__': +################################################### + + importlib.reload(tools) + print('check env variable:') + print('ntask =', os.environ['ntask']) + + + parser = argparse.ArgumentParser(description='create data for MNH et forefire') + parser.add_argument('-i','--input', help='Input MNH dir',required=False) + parser.add_argument('-d','--domainNumber', help='domain number',required=False) + parser.add_argument('-sf','--specificFile', help='to use in case of nested domain',required=False) + args = parser.parse_args() + + ################### ... start input param + + #ff param + #input_ff_param_file = './Inputs/ff-param.nam' + + ################### ... end input param + + #create dir + tools.ensure_dir('./ForeFire/') + tools.ensure_dir('./ForeFire/Outputs/') + #tools.ensure_dir('./MODEL1/') + + dir_input = args.input + if dir_input is None: dir_input = '../006_mnh/' + domainNumber = int(args.domainNumber) if (args.domainNumber is not None) else None + + #read MNH info + MNH_namelist = f90nml.read(dir_input+'/EXSEG1.nam') + expr_name = MNH_namelist['NAM_CONF']['CEXP'] + dtoutput_mnh = MNH_namelist['NAM_BACKUP']['XBAK_TIME_FREQ'] + MNH_domain = tools.read_MNH_domain(expr_name, dtoutput_mnh, dir_mnh = dir_input, domainNumber = domainNumber, specificInputMNHFile=args.specificFile) + + mnhfile = dir_input + 'Postproc/outputFile/Netcdf/{:s}_model{:d}.nc'.format(expr_name,domainNumber) + mnhdata = xr.open_dataset(mnhfile, decode_times=False) + # Extract the time variable + time_var = mnhdata["time"] + #Get the time units, and remove the non-standard prefix + units = time_var.attrs['units'].replace('fire ignition: ', '') + # Convert the time variable using cftime + times = cftime.num2date(time_var.values, units=units, calendar='standard') + # Convert to pandas datetime if needed + times_as_datetime = [datetime(year=t.year, month=t.month, day=t.day, + hour=t.hour, minute=t.minute, second=t.second,microsecond=t.microsecond) + for t in times] + # Replace the time variable in the dataset with the converted times + mnhdata["time"] = ("time", times_as_datetime) + #set mnhdata to have the size of the original domain. add 1pt on the left and 2 on the right + # Calculate spacing for each dimension + dx = mnhdata.x[1] - mnhdata.x[0] # Spacing in the x dimension + dy = mnhdata.y[1] - mnhdata.y[0] # Spacing in the y dimension + dz = mnhdata.z[1] - mnhdata.z[0] # Spacing in the z dimension + # Define new coordinates + new_x = np.concatenate(([mnhdata.x[0] - dx], mnhdata.x, [mnhdata.x[-1] + dx, mnhdata.x[-1] + 2 * dx])) + new_y = np.concatenate(([mnhdata.y[0] - dy], mnhdata.y, [mnhdata.y[-1] + dy, mnhdata.y[-1] + 2 * dy])) + new_z = np.concatenate(([mnhdata.z[0] - dz], mnhdata.z, [mnhdata.z[-1] + dz, mnhdata.z[-1] + 2 * dz])) + # Reindex the dataset, filling new values with -9999 + mnhdata = mnhdata.reindex( + x=new_x, + y=new_y, + z=new_z, + method=None, # Explicitly admnhdata NaNs to new positions + fill_value=-9999 + ) + + # Set nodata value explicitly for all variables + for var in mnhdata.data_vars: + mnhdata[var].attrs['_FillValue'] = -9999 + mnhdata.attrs['nodata'] = -9999 + + dirLCP = '/data/paugam/MNH/Cat_PdV/lcp/' + lcp = xr.open_dataset(dirLCP+'lcp_pgd80_pdV.tif') + + altitude = lcp.isel(band=0).drop_vars(['band','spatial_ref'])['band_data'] # dataarray + altitude.attrs["type"] = "data" + altitude = altitude + + #MERDE + altitude.values[:,:] = 100 + + fuelmap = lcp.isel(band=3).drop_vars(['band','spatial_ref'])['band_data'] # dataarray + fuelmap.attrs["type"] = "fuel" + fuelmap = fuelmap + + ''' + #MERDE + center_x = 408389.97 + center_y = 4617415.01 + + # define the square size (2km x 2km = 2000m x 2000m) + square_size = 2000 # in meters + half_square_size = square_size / 2 + + # define the range for the square + x_min = center_x - half_square_size + x_max = center_x + half_square_size + y_min = center_y - half_square_size + y_max = center_y + half_square_size + + # assuming fuelmap is already defined as an xarray.dataarray + # if not, this would be the time to load the data (replace with your actual fuelmap data) + # fuelmap = xr.open_dataarray('your_fuelmap_file.nc') # example to load a file + + # make a copy of the fuelmap + fuelmap_copy = fuelmap.copy() + + # define the x and y indices within the square bounds + x_indices = (fuelmap.coords['x'] >= x_min) & (fuelmap.coords['x'] <= x_max) + y_indices = (fuelmap.coords['y'] >= y_min) & (fuelmap.coords['y'] <= y_max) + + # modify the values in fuelmap_copy inside the defined square region + fuelmap_copy.loc[{'x': x_indices, 'y': y_indices}] = 99 # set value 99 inside the square area + + # visualize the modified fuelmap_copy to ensure the path is added correctly + #fuelmap_copy.plot() + fuelmap = fuelmap_copy + ''' + #MERDE + fuelmap.values[:,:] = 145 + + heatFlux = fuelmap.copy() + heatFlux[:,:]=0 + heatFlux.attrs["type"] = "flux" + heatFlux.attrs["indices"] = 0 + heatFlux.attrs["model1name"] = "heatFluxBasic" + + #moisture = fuelmap.copy() + #moisture[:,:]=0.1 + #moisture.attrs["type"] = "data" + + vaporFlux = fuelmap.copy() + vaporFlux[:,:]=1 + vaporFlux.attrs["type"] = "flux" + vaporFlux.attrs["indices"] = 1 + vaporFlux.attrs["model1name"] = "vaporFluxBasic" + + #windU = fuelmap.copy() + #windU[:,:]=0 + #windU.attrs["type"] = "data" + + #windV = fuelmap.copy() + #windV[:,:]=4 + #windV.attrs["type"] = "data" + + ds = xr.Dataset({"fuel": fuelmap, "heatFlux": heatFlux, "vaporFlux": vaporFlux, } ) + #ds = xr.Dataset({"altitude":altitude, "fuel": fuelmap, "heatFlux": heatFlux, "vaporFlux": vaporFlux, } ) + #ds = xr.Dataset({"altitude":altitude, "windU":windU, "windV":windV, "fuel": fuelmap, "heatFlux": heatFlux, "vaporFlux": vaporFlux, "moisture": moisture}) + #ds = xr.Dataset({"altitude":altitude, "fuel": fuelmap, "heatFlux": heatFlux, "vaporFlux": vaporFlux, "moisture": moisture}) + ds = ds.rio.write_crs(lcp.rio.crs) + ds = ds.rio.reproject(4326) + + latmean = MNH_domain['lat2d'].mean() + lonmean = MNH_domain['lon2d'].mean() + distance_m = 20 # 1000 meters + lat_offset, lon_offset =tools.offset_in_degrees(distance_m, latmean) + + nlon =int(np.ceil((MNH_domain['lon2d'][0,:].max()-MNH_domain['lon2d'][0,:].min())/lon_offset)) + nlat =int(np.ceil((MNH_domain['lat2d'][:,0].max()-MNH_domain['lat2d'][:,0].min())/lat_offset)) + + print('dimension mydata.nc: lon lat') + print(nlon, nlat) + lon1d = np.linspace(MNH_domain['lon2d'][0,:].min(), MNH_domain['lon2d'][0,:].max(),nlon) + lat1d = np.linspace(MNH_domain['lat2d'][:,0].min(), MNH_domain['lat2d'][:,0].max(),nlat) + + #to MNH domain + ds = ds.interp(x=lon1d, y=lat1d, method='nearest') + + + # then we do some cleaning + ds = ds.rename({'x':'DIMX'}) + ds = ds.rename({'y':'DIMY'}) + ds = ds.expand_dims(DIMZ=1) + ds = ds.expand_dims(DIMT=1) + ds = ds.transpose("DIMT", "DIMZ", "DIMY", "DIMX") + + ds['fuel'] = ds.fuel.astype('int32') + #ds['altitude'] = ds.altitude.astype('float64') + + date_ = datetime.strptime(MNH_domain['date'], "%Y-%m-%d_%H:%M:%S") + + + #add moisture + dirFWI = '/data/paugam/MNH/Cat_PdV/006_mnhSolo/Postproc/FWI/' + fwi = xr.open_dataset(dirFWI + 'FCAST_model3_fwi.nc') + crs = tools.load_crs_from_prj(dirFWI + 'FCAST_model3_fwi.prj') + fwi.rio.write_crs(crs,inplace=True) + + dx_fwi = fwi.x[1]-fwi.x[0] + factor = np.round(dx_fwi/distance_m,0) + lon1d_fwi = np.linspace(MNH_domain['lon2d'][0,:].min(), MNH_domain['lon2d'][0,:].max(),int(nlon/factor)) + lat1d_fwi = np.linspace(MNH_domain['lat2d'][:,0].min(), MNH_domain['lat2d'][:,0].max(),int(nlat/factor)) + + fwi_ll = fwi.rio.reproject(4326) + fwi_ = fwi_ll.interp(x=lon1d_fwi, y=lat1d_fwi, method='nearest') + fwi_ = fwi_.rename({'time':'mnh_t','y':'mnh_y','x':'mnh_x', }) + + #ds['mnh_x'] = (['mnh_x'], lon1d_fwi) + #ds['mnh_y'] = (['mnh_y'], lat1d_fwi) + #ds['mnh_z'] = (['mnh_z'], [0]) + #ds['mnh_t'] = (['mnh_t'], fwi_.mnh_t.values) + #ds['mnh_t'] = (['mnh_t'], [0]) + + #ds['moisture'] = 1.e-2* fwi_.FMC.isel(mnh_t = slice(0,1)).expand_dims(dim=["mnh_z",]) + moisture_ = 1.e-2* fwi_.FMC.expand_dims(dim=["mnh_z",]).transpose('mnh_t', 'mnh_z', 'mnh_y', 'mnh_x') + moisture_ = moisture_.astype(np.float32).fillna(-9999) + + #take the mean + start_time = pd.to_datetime("2022-07-17T11:00:00.000000000") + start_time_s = (start_time - date_).total_seconds() + + end___time = pd.to_datetime("2022-07-17T18:00:00.000000000") + end___time_s = (end___time - date_).total_seconds() + ds['moisture'] = moisture_.sel(mnh_t=slice(start_time, end___time)).mean(dim="mnh_t").expand_dims(dim=["mnh_t",]) + ds['moisture'].attrs["type"] = "data" + + #MERDE + ds['moisture'].values = (np.zeros(ds['moisture'].values.shape)+0.06).astype(np.float32) + + # add wind + dirMNH = '/data/paugam/MNH/Cat_PdV/006_mnhSolo/Postproc/outputFile/Netcdf/' + mnh = xr.open_dataset(dirMNH + 'FCAST_model3_utm.nc',decode_times=False) + crs = tools.load_crs_from_prj(dirMNH + 'FCAST_model3_utm.prj') + mnh.rio.write_crs(crs,inplace=True) + + dx_mnh = mnh.x[1]-mnh.x[0] + factor = np.round(dx_mnh/distance_m,0) + lon1d_mnh = np.linspace(MNH_domain['lon2d'][0,:].min(), MNH_domain['lon2d'][0,:].max(),int(nlon/factor)) + lat1d_mnh = np.linspace(MNH_domain['lat2d'][:,0].min(), MNH_domain['lat2d'][:,0].max(),int(nlat/factor)) + + #mnhu = mnh.u + #mnh_ll = mnhu.rio.reproject(4326) + #mnh_u = mnh_ll.u.interp(x=lon1d_mnh, y=lat1d_mnh, method='nearest') + + reprojected_u = [] + reprojected_v = [] + + # Iterate over each time step + for t in range(mnh.u.shape[0]): # mnhu.shape[0] is the number of time steps + # Reproject the u variable for the current time slice (t) + reprojectedu_slice = mnh.u.isel(time=t).rio.reproject(4326) + reprojectedv_slice = mnh.v.isel(time=t).rio.reproject(4326) + + # Append the reprojected slice to the list + reprojected_u.append(reprojectedu_slice) + reprojected_v.append(reprojectedv_slice) + + # After loop, stack the reprojected slices back together into a new xarray + mnhu_ll = xr.concat(reprojected_u, dim="time") + mnhv_ll = xr.concat(reprojected_v, dim="time") + + mnh_u = mnhu_ll.interp(x=lon1d_mnh, y=lat1d_mnh, method='nearest') + mnh_v = mnhv_ll.interp(x=lon1d_mnh, y=lat1d_mnh, method='nearest') + + mnh_u = mnh_u.rename({'time':'mnh_t','z':'mnh_z', 'y':'mnh_y', 'x':'mnh_x', }) + mnh_v = mnh_v.rename({'time':'mnh_t','z':'mnh_z', 'y':'mnh_y', 'x':'mnh_x', }) + + + + #UNCOMMENT BELOW MERDE + #update altitude with the one of MNH + #mnh_alt = mnh['altitude'].isel(z=0) + #mnh_alt_ll = mnh_alt.rio.reproject(4326) + #mnh_alt_reproj = mnh_alt_ll.interp(x=lon1d_mnh, y=lat1d_mnh, method='nearest') + #ds['altitude'] = mnh_alt_reproj.rename({'y':'mnh_y', 'x':'mnh_x', }).drop_vars(['mnh_x','mnh_y','z']).astype('float32').expand_dims(dim=["mnh_t","mnh_z"]).fillna(-9999) + #ds['altitude'].attrs["type"] = "data" + + + # add extra attr + ds = ds.assign_coords(domdim=("domdim", [1])) # Adding "domdim" as a new dimension + ds["parameters"] = xr.DataArray( + data=["A"], # You can replace "" with actual parameter strings + dims=["domdim"], + attrs={ + "type": "parameters", + "refYear": np.int32(date_.year), + "refDay": np.int32(date_.timetuple().tm_yday), + "duration": "None", + "projection": "None", + "projectionproperties": "None", + }, + ) + ds["domain"] = xr.DataArray( + data=["A"], # Replace with actual data if needed + dims=["domdim"], + attrs={ + "type": "domain", + "NEx": np.float32(MNH_domain['NEx']), + "NEy": np.float32(MNH_domain['NEy']), + "NEz": np.float32(0), + "SWx": np.float32(MNH_domain['SWx']), + "SWy": np.float32(MNH_domain['SWy']), + "SWz": np.float32(0), + "Lx": np.float32(MNH_domain['Lx']), + "Ly": np.float32(MNH_domain['Ly']), + "Lz": np.float32(0), + "t0": np.float32(0), #(ds.mnh_t[0].values -np.datetime64(date_)).item()*1.e-9 , + "Lt": np.float32(np.inf), #(ds.mnh_t[-1].values-np.datetime64(date_)).item()*1.e-9, # Use numpy's representation of infinity + #"t0": np.float32( (ds.mnh_t[0].values -np.datetime64(date_)).item()*1.e-9 ), + #"Lt": np.float32( (ds.mnh_t[-1].values - ds.mnh_t[0].values).item()*1.e-9 ), # Use numpy's representation of infinity + "epsg": lcp.rio.crs.to_epsg(), # Use numpy's representation of infinity + }, + ) + + + #if ignition shp file present in the ForeFire directory, we give the coordinate to input in the Init.ff file + try: + ingitionFile = glob.glob('./ForeFire/ignition/pointIgnition_t0/ignition-pdv.shp')[0] + print('ignition file found') + igni = gpd.read_file(ingitionFile) + igni_utm = igni.to_crs(lcp.rio.crs.to_epsg()) + print(igni) + ref_ll = shapely.Point((MNH_domain['lon2d'][0,0],MNH_domain['lat2d'][0,0]) ) + ref_ll = gpd.GeoDataFrame([{'geometry': ref_ll}], crs="EPSG:4326") + ref_utm = ref_ll.to_crs(lcp.rio.crs.to_epsg()) + + xi = igni_utm.geometry.x - ref_utm.geometry.x + MNH_domain['SWx'] + yj = igni_utm.geometry.y - ref_utm.geometry.y + MNH_domain['SWy'] + + print( "startFire[loc=({:.1f},{:.1f},0.);t=20]".format(xi[0], + yj[0] ) ) + + except: + print('no ignition file') + + try: + referenceFile = glob.glob('./ForeFire/ignition/referencePoints/*.shp') + if len(referenceFile)>0: + print('ref file found') + for referenceFile_ in referenceFile: + print(referenceFile__.split('/')[-1]) + igni = gpd.read_file(referenceFile_) + igni_utm = igni.to_crs(lcp.rio.crs.to_epsg()) + + ref_ll = shapely.Point((MNH_domain['lon2d'][0,0],MNH_domain['lat2d'][0,0]) ) + ref_ll = gpd.GeoDataFrame([{'geometry': ref_ll}], crs="EPSG:4326") + ref_utm = ref_ll.to_crs(lcp.rio.crs.to_epsg()) + + xi = igni_utm.geometry.x - ref_utm.geometry.x + MNH_domain['SWx'] + yj = igni_utm.geometry.y - ref_utm.geometry.y + MNH_domain['SWy'] + + print( "[loc=({:.1f},{:.1f},0.);t=20]".format(xi[0], + yj[0] ) ) + + except: + print('no reference file for point on map') + + + + + windU = mnh_u.isel(mnh_z=slice(0,1), ).fillna(-9999).astype(np.float32) + windU_m = windU.sel(mnh_t=slice(start_time_s, end___time_s)).mean(dim="mnh_t").expand_dims(dim=["mnh_t",]) + windU_m = windU_m.drop_vars('mnh_z') + windU_m = windU_m.drop_vars('mnh_y') + windU_m = windU_m.drop_vars('mnh_x') + ds['windU'] = windU_m + ds['windU'].attrs["type"] = "data" + + #MERDE + #ds['windU'].values = (np.zeros(ds['windU'].values.shape)-0.302).astype(np.float32) + ds['windU'].values = (np.zeros(ds['windU'].values.shape)).astype(np.float32) + + windV = mnh_v.isel(mnh_z=slice(0,1), ).fillna(-9999).astype(np.float32) + windV_m = windV.sel(mnh_t=slice(start_time_s, end___time_s)).mean(dim="mnh_t").expand_dims(dim=["mnh_t",]) + windV_m = windV_m.drop_vars('mnh_z') + windV_m = windV_m.drop_vars('mnh_y') + windV_m = windV_m.drop_vars('mnh_x') + ds['windV'] = windV_m + ds['windV'].attrs["type"] = "data" + + #MERDE + ds['windV'].values = (np.zeros(ds['windV'].values.shape)+1.9207).astype(np.float32) + + ds = ds.drop_vars('DIMX') + ds = ds.drop_vars('DIMY') + ds = ds.drop_vars('domdim') + ds = ds.drop_vars('spatial_ref') + ds = ds.drop_vars('mnh_x') + ds = ds.drop_vars('mnh_y') + #ds = ds.drop_vars('mnh_t') + + #and save + print('save ./ForeFire/mydata.nc') + if os.path.isfile('./ForeFire/mydata.nc'): + os.remove('./ForeFire/mydata.nc') + ds.to_netcdf('./ForeFire/mydata.nc') + diff --git a/tools/preprocessing/create_data_bmap_nc.py b/tools/preprocessing/create_data_bmap_nc.py new file mode 100755 index 00000000..b16f2060 --- /dev/null +++ b/tools/preprocessing/create_data_bmap_nc.py @@ -0,0 +1,1240 @@ +from __future__ import print_function +from __future__ import division +from builtins import zip +from builtins import range +from past.utils import old_div +import sys +import os +import numpy as np +import matplotlib as mpl +mpl.use('Qt5Agg') +import matplotlib.pyplot as plt +import matplotlib.path as mpath +import matplotlib.patches as mpatches +import pdb +#import imp +from scipy import io +import glob +import netCDF4 +import f90nml +import datetime +import socket +import warnings +warnings.filterwarnings("ignore",".*GUI is implemented.*") +from mpl_toolkits.axes_grid1 import make_axes_locatable +import pickle +import argparse +from osgeo import gdal,osr,ogr +from scipy import interpolate +import importlib +import pandas + + +path_src = os.environ['PATH_SRC_PYTHON_LOCAL'] + +try: + path_ForeFire = os.environ["PATH_FOREFIRE"] +except: + print('PATH_FOREFIRE is not defined. stop here') + sys.exit() + +sys.path.append(path_src+'/Regrid/') +#sys.path.append(path_ForeFire) +#sys.path.append(path_ForeFire+'/swig/') +#sys.path.append(path_ForeFire+'tools/') +sys.path.append(path_ForeFire+'tools/preprocessing/') +try: + import pyforefire as forefire + flag_run_ff = True +except ImportError: + print('could not find ForeFire API') + print('check variable path_ForeFire in path_ForeFire/tools/create_data_bmap_nc.py') + print('currently, path_ForeFire = ', path_ForeFire) + print('flag_run_ff set to False') + flag_run_ff = False +import regrid +importlib.reload(regrid) +import parametersProperties as attribute + + +################################################## +def ensure_dir(f): + d = os.path.dirname(f) + if not os.path.exists(d): + os.makedirs(d) + + +################################################## +def read_MNH_domain( expr_name, dir_mnh = '../02_mnh/', domainNumber=None, specificInputMNHFile=None): + + if specificInputMNHFile == None: + expr_numb = '.{:d}.'.format(domainNumber) if (domainNumber is not None) else '' + outfiles_mnh = sorted(list(set(glob.glob(dir_mnh+expr_name[:5]+expr_numb+'*.nc'))-\ + set(glob.glob(dir_mnh+expr_name[:5]+expr_numb+'*.OUT.*.nc')))+\ + glob.glob(dir_mnh+expr_name[:5]+expr_numb+'*.nc4')) + print(outfiles_mnh) + if len(outfiles_mnh) == 0: + print('missing mnh files') + pdb.set_trace() + outfiles_mnh_inputTime = outfiles_mnh[0] + outfiles_mnh_inputGrid = outfiles_mnh[1] + + else: + outfiles_mnh_inputGrid = specificInputMNHFile + outfiles_mnh_inputTime = specificInputMNHFile + + print('get time ref from :', outfiles_mnh_inputTime) + print('get grid from :', outfiles_mnh_inputGrid) + + + #read geometry of the MNH domain + nc = netCDF4.Dataset(outfiles_mnh_inputGrid,'r') # get the 001 file + + if domainNumber is not None: + expr_numb = '.{:d}.'.format(1) + outfiles_mnh_dad = sorted(list( set(glob.glob(dir_mnh+expr_name[:5]+expr_numb+'*.nc'))-\ + set(glob.glob(dir_mnh+expr_name[:5]+expr_numb+'*.OUT.*.nc')))+\ + glob.glob(dir_mnh+expr_name[:5]+expr_numb+'*.nc4')) + ncDad = netCDF4.Dataset(outfiles_mnh_dad[1],'r') + + MNH_properties= {} + + try: + if domainNumber is not None: + MNH_properties['lat00'] = np.array(ncDad.variables['LAT'])[0,0] + MNH_properties['lon00'] = np.array(ncDad.variables['LON'])[0,0] + ncDad.close() + else: + MNH_properties['lat00'] = np.array(nc.variables['LAT'])[0,0] + MNH_properties['lon00'] = np.array(nc.variables['LON'])[0,0] + except: + print('no ref point found, assume IDEAL RUN') + MNH_properties['lat00'] = 'ideal' + + DeltaY = nc.variables['YHAT'][1]-nc.variables['YHAT'][0] + DeltaX = nc.variables['XHAT'][1]-nc.variables['XHAT'][0] + + if nc.variables['MASDEV'][0].data==53: + MNH_properties['nx'] = len(nc.dimensions[u'X']) + MNH_properties['ny'] = len(nc.dimensions[u'Y']) + MNH_properties['nz'] = len(nc.dimensions[u'Z']) + else: + MNH_properties['nx'] = len(nc.dimensions[u'ni']) + MNH_properties['ny'] = len(nc.dimensions[u'nj']) + MNH_properties['nz'] = len(nc.dimensions[u'level_w']) + + MNH_properties['SWx'] = nc.variables['XHAT'][0]+0. + MNH_properties['SWy'] = nc.variables['YHAT'][0]+0. + MNH_properties['SWz'] = 0 + MNH_properties['NEx'] = nc.variables['XHAT'][-1]+DeltaX + MNH_properties['NEy'] = nc.variables['YHAT'][-1]+DeltaY + MNH_properties['NEz'] = 0 + + MNH_properties['Lx'] = nc.variables['XHAT'][-1]+DeltaX-nc.variables['XHAT'][0] + MNH_properties['Ly'] = nc.variables['YHAT'][-1]+DeltaY-nc.variables['YHAT'][0] + MNH_properties['Lz'] = 0 + nc.close() + + nc = netCDF4.Dataset(outfiles_mnh_inputTime,'r') + if nc.variables['MASDEV'][0].data==53: + MNH_properties['date'] = datetime.datetime(*nc.variables['DTCUR__TDATE'][:]).strftime("%Y-%m-%d_%H:%M:%S") + else: + MNH_properties['date'] = '_'.join(nc.variables['DTCUR'].units.split(' ')[2:4]) + + + if nc.variables['MASDEV'][0].data==53: + MNH_properties['t0'] = nc.variables['DTCUR__TIME'][0] + else: + MNH_properties['t0'] = nc.variables['DTCUR'][0] + MNH_properties['Lt'] = np.Inf + + nc.close() + + + return MNH_properties + + +################################################### +def get_dimension_name_and_value(key,nc_variable, nco): + + dim_var = np.array( list(nc_variable[key].keys()) ) + dim_var_ = np.array([dimension[:-1] for dimension in nc_variable[key] ]) + idx_dim_var = np.where(dim_var_ == 'dim') + + dimensions_dimX = tuple(sorted(dim_var[idx_dim_var])[::-1]) + dimensions_name = [] + for dim in dimensions_dimX: + dimensions_name.append(nc_variable[key][dim]) + + nco.createVariable(key,nc_variable[key]['type'], dimensions_name ) + + dimensions_value = [ nco.dimensions[dim] for dim in dimensions_name ] + + return dimensions_name, dimensions_value + + +################################################### +def merge_two_dicts(x, y): + '''Given two dicts, merge them into a new dict as a shallow copy.''' + z = x.copy() + z.update(y) + return z + + +################################# +def getLocationFromLine(line): + llv = line.split("loc=(") + if len(llv) < 2: + return None + llr = llv[1].split(","); + if len(llr) < 3: + return None + return (float(llr[0]),float(llr[1])) + + +################################# +def dist(a,b): + return np.sqrt(np.power(a[0]-b[0],2)+np.power(a[1]-b[1],2)) + + +################################# +def printToPathe(linePrinted): + fronts = linePrinted.split("FireFront") + pathes = [] + for front in fronts[1:]: + nodes =front.split("FireNode")[1:] + if len(nodes) > 0: + Path = mpath.Path + + codes = [] + verts = [] + lastNode = getLocationFromLine(nodes[0]) + firstNode = getLocationFromLine(nodes[0]) + + + codes.append(Path.MOVETO) + verts.append(firstNode) + + for node in nodes[:]: + newNode = getLocationFromLine(node) + codes.append(Path.LINETO) + verts.append(newNode) + lastNode = newNode + + + codes.append(Path.LINETO) + verts.append(firstNode) + + pathes.append(mpath.Path(verts, codes)) + + return pathes; + + +################################################### +def getDomainExtent(line): + print(line) + llv = line.split("sw=(") + llr = llv[1].split("ne=("); + return( float( llr[0].split(",")[0]), float(llr[1].split(",")[0]), float(llr[0].split(",")[1]), float(llr[1].split(",")[1]) ) + + +################################################### +if __name__ == '__main__': +################################################### + + print('check env variable:') + print('ntask =', os.environ['ntask']) + + + parser = argparse.ArgumentParser(description='create data for MNH et forefire') + parser.add_argument('-i','--input', help='Input MNH dir',required=False) + parser.add_argument('-d','--domainNumber', help='domain number',required=False) + parser.add_argument('-sf','--specificFile', help='to use in case of nested domain',required=False) + args = parser.parse_args() + + ################### ... start input param + + #ff param + input_ff_param_file = './Inputs/ff-param.nam' + + ################### ... end input param + + #create dir + ensure_dir('./Inputs/') + ensure_dir('./ForeFire/') + ensure_dir('./ForeFire/Outputs/') + ensure_dir('./MODEL1/') + + dir_input = args.input + if dir_input is None: dir_input = '../02_mnh/' + domainNumber = int(args.domainNumber) if (args.domainNumber is not None) else None + + #read MNH info + MNH_namelist = f90nml.read(dir_input+'/EXSEG1.nam') + expr_name = MNH_namelist['NAM_CONF']['CEXP'] + MNH_domain = read_MNH_domain(expr_name, dir_mnh = dir_input, domainNumber = domainNumber, specificInputMNHFile=args.specificFile) + ''' + try: + + except IOError: + print '************** set test run ****************' + expr_name = 'test' + MNH_domain= {} + MNH_domain['nx'] = 10 + MNH_domain['ny'] = 10 + MNH_domain['nz'] = 10 + + MNH_domain['SWx'] = 0 + MNH_domain['SWy'] = 0 + MNH_domain['SWz'] = 0 + MNH_domain['NEx'] = 100 + MNH_domain['NEy'] = 100 + MNH_domain['NEz'] = 100 + + MNH_domain['Lx'] = None + MNH_domain['Ly'] = None + MNH_domain['Lz'] = None + + MNH_domain['t0'] = 0. + MNH_domain['Lt'] = np.Inf + + MNH_domain['date'] = datetime.datetime(1970,1,1).strftime("%Y-%m-%d_%H:%M:%S") + ''' + + #read FF info + nmls = f90nml.read('./Inputs/ff-param.nam') + + #ronan ff info + ronan_param = nmls['FF_RONAN'] + print('############') + if not(ronan_param['flag_freshStart']): + print(' this is not a fresh start') + print(' init mod is :', ronan_param['init_mode']) + print('############') + + + if ronan_param['init_mode'] == 'ideal': + print('initialization in ideal mode') + + flag_rotation = False + shiftTheta = -999 + rotce = -999 + rotcn = -999 + + #ideal fix FireLine + x_location_fireLine = ronan_param['x_location_fireLine'] + length_fireLine = ronan_param['length_fireLine'] + depth_fireLine = ronan_param['depth_fireLine'] + time_start_fireLine = ronan_param['time_start_fireLine'] + residence_time = ronan_param['residence_time'] + burning_time = ronan_param['burning_time'] + + moisture = ronan_param['moisture'] + nominalHeatFlux_f = ronan_param['nominalHeatFlux_flaming'] + nominalHeatFlux_s = ronan_param['nominalHeatFlux_smoldering'] + + if type(time_start_fireLine) is float: + nbre_fire_line = 1 + x_location_fireLine = [x_location_fireLine] + length_fireLine = [length_fireLine] + depth_fireLine = [depth_fireLine] + time_start_fireLine = [time_start_fireLine] + residence_time = [residence_time] + burning_time = [burning_time] + moisture = [moisture] + nominalHeatFlux_f = [nominalHeatFlux_f] + nominalHeatFlux_s = [nominalHeatFlux_s] + else: + nbre_fire_line = len(time_start_fireLine) + flag_extra_contour = ronan_param['flag_extra_contour'] #extra parameter + shiftx = 0 + shifty = 0 + shiftz = 0 + + + if ronan_param['init_mode'] == 'fromOverHeadImg': + + firedata = np.load(ronan_param['init_file']) + firedata = firedata.view(np.recarray) + if MNH_domain['lat00'] == 'ideal': + shiftx = firedata.grid_e[0,0] - ronan_param['init_file_xshift'] + shifty = firedata.grid_n[0,0] - ronan_param['init_file_yshift'] + else: + prj_info = open(ronan_param['init_file'].replace('.npy','.prj')).readline() + wgs84 = osr.SpatialReference( ) # Define a SpatialReference object + wgs84.ImportFromEPSG( 4326 ) # And set it to WGS84 using the EPSG code + utm = osr.SpatialReference() + utm.ImportFromWkt(prj_info) + print('********* WARNING **********: need to change grid UTM zone in next obs data, hard reset to Zone 36N to correct error') + utm.SetUTM(36) + conv_ll2utm = osr.CoordinateTransformation(wgs84, utm) + conv_utm2ll = osr.CoordinateTransformation(utm,wgs84) + + e_00, n_00, _ = conv_ll2utm.TransformPoint(MNH_domain['lon00'], MNH_domain['lat00']) + + shiftx = e_00 + shifty = n_00 + #shiftx = firedata.grid_e[0,0] - ronan_param['init_file_xshift'] - e_00 + #shifty = firedata.grid_n[0,0] - ronan_param['init_file_yshift'] - n_00 + shiftz = 0 + + firedata.grid_e = firedata.grid_e - shiftx + firedata.grid_n = firedata.grid_n - shifty + + #apply rotation + if 'init_file_thetashift' in ronan_param.keys(): + flag_rotation = True + shiftTheta = ronan_param['init_file_thetashift'] + + idx = np.where(firedata.arrivalTime>0) + rotce, rotcn = firedata.grid_e[idx].mean(), firedata.grid_n[idx].mean() + + firedata.grid_e -= rotce + firedata.grid_n -= rotcn + + rot = np.zeros([2,2]) + theta = 3.14*shiftTheta/180. + rot[0,0] = np.cos(theta) + rot[0,1] = - np.sin(theta) + rot[1,0] = np.sin(theta) + rot[1,1] = np.cos(theta) + + nn = firedata.grid_e.flatten().shape[0] + ptsRot = np.zeros([nn,2]) + for ii_, [xx,yy] in enumerate( zip(firedata.grid_e.flatten(), firedata.grid_n.flatten())): + ptsRot[ii_] = rot.dot(np.array([xx,yy])) + firedata.grid_e = ptsRot[:,0].reshape(firedata.grid_e.shape) + firedata.grid_n = ptsRot[:,1].reshape(firedata.grid_e.shape) + + firedata.grid_e += rotce + firedata.grid_n += rotcn + + else: + flag_rotation = False + shiftTheta = -999 + rotce = -999 + rotcn = -999 + + + dx = np.sqrt( (firedata.grid_e[1,0]-firedata.grid_e[0,0])**2+ (firedata.grid_n[1,0]-firedata.grid_n[0,0])**2 ) + dy = np.sqrt( (firedata.grid_e[0,1]-firedata.grid_e[0,0])**2+ (firedata.grid_n[0,1]-firedata.grid_n[0,0])**2 ) + firedata_reso = np.ones_like(firedata.grid_e)*dx*dy + + #flag no fire pixel as -999 + idx = np.where((firedata.fre_f+firedata.fre_s)<=0) + firedata.fre_f[idx] = 0 + firedata.fre_s[idx] = 0 + firedata.residenceTime[idx] = -999 + firedata.burningTime[idx] = -999 + firedata.arrivalTime[idx] = -999 + + flag_extra_contour = False # not use here + + if flag_run_ff != None: + flag_run_ff = ronan_param['flag_run_ff'] + evaporation_time = ronan_param['evaporation_time'] + + #general ff info + ff_param = nmls['FF_PARAM'] + minimalPropagativeFrontDepth = ff_param['minimalPropagativeFrontDepth'] + spatialIncrement = ff_param['spatialIncrement'] + perimeterResolution = ff_param['perimeterResolution'] + outputsUpdate = ff_param['outputsUpdate'] + bmapOutputUpdate = ff_param['bmapOutputUpdate'] + # end time of the bmap. if only forcing from the bmap use large value > simulation time + InitTime = ff_param['InitTime'] + + # dimension + nx = MNH_domain['nx'] + ny = MNH_domain['ny'] + + #load attribute + for key in ['Lx','Ly','Lz','NEx','NEy','NEz','SWx','SWy','SWz','t0','Lt']: + attribute.domain[key] = MNH_domain[key] + #attribute.parameters['date'] = MNH_domain['date'] + date_mnh = datetime.datetime.strptime(MNH_domain['date'], '%Y-%m-%d_%H:%M:%S') + #time_mnh = attribute.domain['t0'] + attribute.parameters['refYear'] = date_mnh.year + #attribute.parameters['refMonth'] = date_mnh.month + attribute.parameters['refDay'] = date_mnh.timetuple().tm_yday #date_mnh.day + + #'projection','duration' are not defined here + + #add domain length to attribute + Lx = 1.*(attribute.domain['NEx']-attribute.domain['SWx']) + Ly = 1.*(attribute.domain['NEy']-attribute.domain['SWy']) + attribute.domain['Lx'] = Lx + attribute.domain['Ly'] = Ly + + # create atmo grid + ################## + atmo_dx = old_div(Lx,nx) + atmo_dy = old_div(Ly,ny) + y_center = .5 * (Ly - 2 * atmo_dy) #of the domain without the grid pt on the side + atmo_grid_x = np.arange(attribute.domain['SWx'],attribute.domain['NEx']+atmo_dx, atmo_dx) + atmo_grid_y = np.arange(attribute.domain['SWy'],attribute.domain['NEy']+atmo_dy, atmo_dy) + #2D grid + atmo_grid_y2d, atmo_grid_x2d = np.meshgrid(atmo_grid_y, atmo_grid_x) + + + # create bmap grid + ################## + #this need to be defined here for bmap matrix + BMapsResolution = max(old_div(spatialIncrement,np.sqrt(2)), minimalPropagativeFrontDepth) + print('BMapsResolution = ', BMapsResolution) + + nbre_firePt_per_atmoCell_x = (int(old_div(atmo_dx,BMapsResolution)) ) + bmap_dx = old_div(atmo_dx, nbre_firePt_per_atmoCell_x) + bmap_grid_x = np.arange(attribute.domain['SWx'],attribute.domain['NEx']+bmap_dx, bmap_dx) + + nbre_firePt_per_atmoCell_y = (int(old_div(atmo_dy,BMapsResolution)) ) + bmap_dy = old_div(atmo_dy, nbre_firePt_per_atmoCell_y) + bmap_grid_y = np.arange(attribute.domain['SWy'],attribute.domain['NEy']+bmap_dy, bmap_dy) + + #2D grid + bmap_grid_y2d, bmap_grid_x2d = np.meshgrid(bmap_grid_y, bmap_grid_x) + + if ronan_param['init_mode'] == 'ideal': + #create firedata array from input info + firedata_grid_x = np.arange(attribute.domain['SWx'],attribute.domain['NEx']+ronan_param['firedata_res'], ronan_param['firedata_res']) + firedata_grid_y = np.arange(attribute.domain['SWy'],attribute.domain['NEy']+ronan_param['firedata_res'], ronan_param['firedata_res']) + firedata_grid_y2d, firedata_grid_x2d = np.meshgrid(firedata_grid_y, firedata_grid_x) + firedata = np.zeros([firedata_grid_y2d.shape[0]-1,firedata_grid_y2d.shape[1]-1],dtype=np.dtype([('grid_e', '= xlocfire) &\ + (firedata.grid_e < xlocfire+depth_fireLine[i_fireLine]) &\ + (firedata.grid_n >= y_center-.5*length_fireLine[i_fireLine]) &\ + (firedata.grid_n < y_center+.5*length_fireLine[i_fireLine]) ) + + firedata.residenceTime[idx] = residence_time[i_fireLine] + firedata.burningTime[idx] = burning_time[i_fireLine] + firedata.moisture[idx] = moisture[i_fireLine] + firedata.arrivalTime[idx] = time_start_fireLine[i_fireLine] + firedata.fre_f[idx] = old_div(nominalHeatFlux_f[i_fireLine]*ff_param['fractionRadiation_f'],(1.-ff_param['fractionRadiation_f']))\ + *ronan_param['firedata_res']*ronan_param['firedata_res']*residence_time[i_fireLine] * 1.e-6 # MJ + firedata.fre_s[idx] = old_div(nominalHeatFlux_s[i_fireLine]*ff_param['fractionRadiation_s'],(1.-ff_param['fractionRadiation_s']))\ + *ronan_param['firedata_res']*ronan_param['firedata_res']*(burning_time[i_fireLine]-residence_time[i_fireLine]) * 1.e-6 # MJ + + + #match pixels from firedata grid and atmo grid + ####################### + regrid_name = 'fireData2atmo' + wkdir = './' + gridin_polygon, gridin_xx, gridin_yy = regrid.get_polygon_from_grid(firedata.grid_e, firedata.grid_n, flag_add_1extra_rowCol=True) + gridout_polygon, gridout_xx, gridout_yy = regrid.get_polygon_from_grid(atmo_grid_x2d, atmo_grid_y2d ) + + in_outA_grid_list, outA_in_grid_list = regrid.grid2grid_pixel_match(regrid_name, \ + gridin_xx, gridin_yy, gridin_polygon,None, \ + gridout_polygon, \ + wkdir, \ + flag_parallel=ronan_param['flag_parallel'], \ + flag_freshstart=ronan_param['flag_freshStart']) + + + #init param + windU_in = 10 # cte wind speed + windV_in = 0 + mydatanc = './Inputs/mydata_'+expr_name+'.nc' + mybmapnc = './Inputs/mybmap_'+expr_name+'.nc' + + #set model name + attribute.heatFlux['model0name'] = ronan_param['heatFluxModel'] + attribute.vaporFlux['model1name'] = ronan_param['vaporFluxModel'] + attribute.scalarFlux['model2name'] = ronan_param['scalarFluxModel'] + + + #--------------------------------------------- + #create data nc file for ForeFire Layer Models + #--------------------------------------------- + print('write data file: ', mydatanc.split('/')[-1]) + fnameout = mydatanc + nc_dimension = {'DIMX': -999, 'DIMY': -999, 'DIMZ': -999, 'DIMT': -999, 'domdim': -999} + + + var4dim = {'dim1': 'DIMX', 'dim2': 'DIMY', 'dim3': 'DIMZ', 'dim4': 'DIMT', 'type': 'float' } + nc_variable_data = {'vaporFlux' : var4dim, \ + 'heatFlux' : var4dim, \ + 'scalarFlux' : var4dim, \ + 'fuel' : var4dim, \ + 'windU' : var4dim, \ + 'windV' : var4dim, \ + 'altitude' : var4dim, \ + #fire regime control time + 'FromObs_residenceTime' : var4dim,\ + 'FromObs_burningTime' : var4dim,\ + 'FromObs_evaporationTime' : var4dim,\ + #sensible heat and Vapor flux + 'FromObs_NominalHeatFlux_flaming' : var4dim,\ + 'FromObs_NominalHeatFlux_smoldering' : var4dim,\ + 'FromObs_VaporFlux_evaporation' : var4dim,\ + #Emission Factor + 'FromObs_EFScalar_flaming': var4dim,\ + 'FromObs_EFScalar_smoldering': var4dim,\ + #radiation Fraction + 'FromObs_radiationFraction_flaming' : var4dim,\ + 'FromObs_radiationFraction_smoldering' : var4dim,\ + #convsersion factor + 'FromObs_conversionFactor' : var4dim,\ + #misc + 'parameters': {'dim1': 'domdim', 'type': 'c'}, + 'domain' : {'dim1': 'domdim', 'type': 'c'} } + + nc_attribute = {'vaporFlux': attribute.vaporFlux, + 'heatFlux' : attribute.heatFlux,\ + 'scalarFlux' : attribute.scalarFlux,\ + 'fuel' : attribute.fuel, \ + 'windU' : attribute.windU, \ + 'windV' : attribute.windV, \ + 'altitude' : attribute.altitude, \ + #fire regime control time + 'FromObs_residenceTime' : attribute.FromObs_residenceTime, \ + 'FromObs_burningTime' : attribute.FromObs_burningTime, \ + 'FromObs_evaporationTime' : attribute.FromObs_evaporationTime, \ + #sensible heat & vapor8flux + 'FromObs_NominalHeatFlux_flaming' : attribute.FromObs_NominalHeatFlux_flaming, \ + 'FromObs_NominalHeatFlux_smoldering' : attribute.FromObs_NominalHeatFlux_smoldering, \ + 'FromObs_VaporFlux_evaporation' : attribute.FromObs_VaporFlux_evaporation, \ + #Emission Factor + 'FromObs_EFScalar_flaming': attribute.FromObs_EFScalar_flamning, \ + 'FromObs_EFScalar_smoldering': attribute.FromObs_EFScalar_smoldering, \ + #radiation fraction + 'FromObs_radiationFraction_flaming': attribute.FromObs_radiationFraction_flaming, \ + 'FromObs_radiationFraction_smoldering': attribute.FromObs_radiationFraction_smoldering, \ + #conversion factor + 'FromObs_conversionFactor': attribute.FromObs_conversionFactor, \ + #misc + 'parameters': attribute.parameters, + 'domain' : attribute.domain} + + + #open nc file + nco = io.netcdf_file(fnameout, 'w') + + #set up dimension + for dim in list(nc_dimension.keys()): + if dim == 'DIMX': + nco.createDimension(dim, nx ) + elif dim == 'DIMY': + nco.createDimension(dim, ny ) + elif dim == 'DIMZ': + nco.createDimension(dim, 1 ) + elif dim == 'DIMT': + nco.createDimension(dim, 1 ) + elif dim == 'domdim': + nco.createDimension(dim, 1 ) + else: + print('should not be here') + pdb.set_trace() + + #set sensible heat flux at the resoluton of mnh + #------------- + _, dimensions_value = get_dimension_name_and_value('FromObs_residenceTime',nc_variable_data, nco) # get dimension + residenceTime_out, _ = regrid.map_data(outA_in_grid_list,dimensions_value[2:][::-1],firedata.residenceTime,flag='average') #[::-1] as in nc dimension are stored in inverse order + burningTime_out, _ = regrid.map_data(outA_in_grid_list,dimensions_value[2:][::-1],firedata.burningTime,flag='average') + evaporationTime_out = np.zeros_like(burningTime_out) + ronan_param['evaporation_time'] + + #afArea_in = np.where(firedata.fre>0,firedata_reso,np.zeros_like(firedata.fre)) + #afArea_out, _ = regrid.map_data(outA_in_grid_list,dimensions_value[2:][::-1],afArea_in,flag='sum',gridReso_in=firedata_reso) + + fre_f_out, _ = regrid.map_data(outA_in_grid_list,dimensions_value[2:][::-1],firedata.fre_f,flag='sum',gridReso_in=firedata_reso) + fre_f_out *= 1.e6 # J + + + fre_s_out, _ = regrid.map_data(outA_in_grid_list,dimensions_value[2:][::-1],firedata.fre_s,flag='sum',gridReso_in=firedata_reso) + fre_s_out *= 1.e6 # J + + ''' + #remove pixel with low burningTime and residenceTime + idx = np.where( (residenceTime_out<1) ) #& (fre_f_out>0)) + fre_f_out[idx] = 0 + residenceTime_out[idx] = 0 + + idx = np.where( ((burningTime_out-residenceTime_out)<1) )#& ((fre_f_out+fre_s_out)>0)) # smouldering time < 1 + fre_f_out[idx] = 0 + fre_s_out[idx] = 0 + burningTime_out[idx] = 0 + residenceTime_out[idx] = 0 + ''' + + #clip value of residence and burning time that do not fit fre + idx = np.where((residenceTime_out<0)|(burningTime_out<=residenceTime_out)) + for ii,jj in zip(idx[0],idx[1]): + if (fre_f_out[ii,jj] > 0) & (fre_s_out[ii,jj] > 0) : + fre_s_out[ii,jj] += fre_f_out[ii,jj] # pass flaming fre to smouldering + fre_f_out[ii,jj] = 0 + residenceTime_out[ii,jj] = 0 #. + elif (fre_f_out[ii,jj] > 0) & (fre_s_out[ii,jj] == 0) : + residenceTime_out[ii,jj] = residenceTime_out[np.where(residenceTime_out>0)].max() + burningTime_out[ii,jj] = residenceTime_out[ii,jj] + elif (fre_f_out[ii,jj] == 0) & (fre_s_out[ii,jj] > 0) : + burningTime_out[ii,jj] = burningTime_out[np.where(burningTime_out>0)].max() + residenceTime_out[ii,jj] = 0. + elif (fre_f_out[ii,jj] == 0) & (fre_s_out[ii,jj] == 0) : + burningTime_out[ii,jj] = 0. + residenceTime_out[ii,jj] = 0. + + + #mass comsumption + mass_comsumption_out = ff_param['conversionFactor'] * 1.e-6 * (fre_f_out + fre_s_out) # kg + + + + #flaming sensible heat flux + print('compute flaming senible heat flux') + idx = np.where(fre_f_out>0) + sensible_heat_flux_f = np.zeros_like(fre_f_out) + sensible_heat_flux_f[idx] = old_div((old_div((1.-ff_param['fractionRadiation_f']),ff_param['fractionRadiation_f']))\ + *fre_f_out[idx],(residenceTime_out[idx]* atmo_dx * atmo_dy)) # W/m2 + if len(np.where(sensible_heat_flux_f<0)[0])>0: + print('issue in flaming sensible heat flux. found <0 values') + sys.exit() + if sensible_heat_flux_f.max() > 1.e10: + pdb.set_trace() + + + #patch to remove high flux + #idx = np.where(sensible_heat_flux_f>7.e5) + #sensible_heat_flux_f[idx]=0 + #residenceTime_out[idx]=0 + + + #smoldering sensible heat flux + print('compute smoldering senible heat flux') + idx = np.where(fre_s_out>0) + sensible_heat_flux_s = np.zeros_like(fre_s_out) + sensible_heat_flux_s[idx] = old_div((old_div((1.-ff_param['fractionRadiation_s']),ff_param['fractionRadiation_s']))\ + *fre_s_out[idx],((burningTime_out[idx]-residenceTime_out[idx])* atmo_dx * atmo_dy)) # W/m2 + if len(np.where(sensible_heat_flux_s<0)[0])>0: + print('issue in smoldering sensible heat flux. found <0 values') + sys.exit() + + #patch to remove high flux + #idx = np.where(sensible_heat_flux_s>=4.e5) + #sensible_heat_flux_s[idx]=0 + #burningTime_out[idx]= residenceTime_out[idx] + + + #Vapor flux during evaporation + print('compute vapor flux') + moisture_out, _ = regrid.map_data(outA_in_grid_list,dimensions_value[2:][::-1],firedata.moisture,flag='average') + moisture_out = np.where(moisture_out<0,np.zeros_like(moisture_out),moisture_out) + + vapor_flux = np.zeros_like(moisture_out) + idx = np.where(moisture_out*mass_comsumption_out > 0) + vapor_flux[idx] = old_div((1.e-2*moisture_out[idx]*mass_comsumption_out[idx]),(evaporationTime_out[idx]* atmo_dx * atmo_dy)) # kg/m2 + + + #create variables + #-------------- + for key in nc_variable_data: + dimensions_name, dimensions_value = get_dimension_name_and_value(key,nc_variable_data, nco) + #nco.variables[key][:] = np.zeros(dimensions_value) + + if key == 'windU': + nco.variables[key][:] = np.zeros(dimensions_value) + windU_in + elif key == 'windV': + nco.variables[key][:] = np.zeros(dimensions_value) + windV_in + elif key == 'vaporFlux': + nco.variables[key][:] = np.zeros(dimensions_value) + 1 + elif key == 'heatFlux': + nco.variables[key][:] = np.zeros(dimensions_value) + 0 + elif key == 'scalarFlux': + nco.variables[key][:] = np.zeros(dimensions_value) + 2 + + elif key == 'FromObs_residenceTime': + print(' map residence time:') + nco.variables[key][:] = np.array(residenceTime_out.T,dtype=np.float32).reshape(dimensions_value) + + elif key == 'FromObs_evaporationTime': + print(' ********** WARNING ***********: map evaporation time is set constant') + nco.variables[key][:] = np.array(evaporationTime_out.T,dtype=np.float32).reshape(dimensions_value) + + elif key == 'FromObs_burningTime': + print(' map burning time:') + nco.variables[key][:] = np.array(burningTime_out.T,dtype=np.float32).reshape(dimensions_value) + + elif key == 'FromObs_NominalHeatFlux_flaming': + print(' map flaming nominal heat flux:') + nco.variables[key][:] = np.array(sensible_heat_flux_f.T,dtype=np.float32).reshape(dimensions_value) + + elif key == 'FromObs_NominalHeatFlux_smoldering': + print(' map smoldering nominal heat flux:') + nco.variables[key][:] = np.array(sensible_heat_flux_s.T,dtype=np.float32).reshape(dimensions_value) + + elif key == 'FromObs_VaporFlux_evaporation': + print(' map Vapor flux:') + nco.variables[key][:] = np.array(vapor_flux.T,dtype=np.float32).reshape(dimensions_value) + #nco.variables[key][:] = np.array(np.zeros_like(vapor_flux).T,dtype=np.float32).reshape(dimensions_value) + dimensions_value_EVAP = dimensions_value + + elif key == 'FromObs_EFScalar_flaming': + print(' ********** WARNING ***********: EF scalar flaming is set constant') + nco.variables[key][:] = np.zeros(dimensions_value) + 0.01 + + elif key == 'FromObs_EFScalar_smoldering': + print(' ********** WARNING ***********: EF scalar smoldering is set constant') + nco.variables[key][:] = np.zeros(dimensions_value) + 0.02 + + elif key == 'FromObs_radiationFraction_flaming': + print(' map fraction radiation for smoldering:') + nco.variables[key][:] = np.zeros(dimensions_value) + ff_param['fractionRadiation_f'] + + elif key == 'FromObs_radiationFraction_smoldering': + print(' map fraction radiation for smoldering:') + nco.variables[key][:] = np.zeros(dimensions_value) + ff_param['fractionRadiation_s'] + + elif key == 'FromObs_conversionFactor': + print(' map consersion factor:') + nco.variables[key][:] = np.zeros(dimensions_value) + ff_param['conversionFactor'] + + elif key == 'fuel': + print(' map fuel:') + nco.variables[key][:] = np.zeros(dimensions_value) + 1 + + else: + nco.variables[key][:] = np.zeros(dimensions_value) + + for attvar in nc_attribute[key]: + nco.variables[key]._attributes[attvar] = nc_attribute[key][attvar] + + + #close nc file + nco.close() + + + #--------------------------------------------- + # create bmap nc file to force fire propapation + #--------------------------------------------- + print('write bmap file: ', mybmapnc.split('/')[-1]) + fnameout = mybmapnc + nc_dimension = {'DIMX': -999, 'DIMY': -999, 'C_DIMX': -999, 'C_DIMY': -999, 'domdim': -999} + + var2dim_atm = {'dim1': 'C_DIMX', 'dim2': 'C_DIMY', 'type': 'int32' } + var2dim_fire = {'dim1': 'DIMX', 'dim2': 'DIMY', 'type': 'float' } + + nc_variable_bmap = {'arrival_time_of_front' : var2dim_fire, \ + 'cell_active' : var2dim_atm, + 'domain' : {'dim1': 'domdim', 'type': 'c'} } + + nc_attribute = {'arrival_time_of_front': None, + 'cell_active' : None, \ + 'domain' : merge_two_dicts(attribute.domain,attribute.parameters)} + + + + #sys.exit() + #open nc file + nco = io.netcdf_file(fnameout, 'w') + + #set up dimension + for dim in list(nc_dimension.keys()): + if dim == 'C_DIMX': + nco.createDimension(dim, nx ) + elif dim == 'C_DIMY': + nco.createDimension(dim, ny ) + elif dim == 'DIMX': + nco.createDimension(dim, nbre_firePt_per_atmoCell_x * (nx)) + elif dim == 'DIMY': + nco.createDimension(dim, nbre_firePt_per_atmoCell_y * (ny)) + elif dim == 'domdim': + nco.createDimension(dim, 1 ) + else: + print('should not be here') + pdb.set_trace() + + + #match pixels from firedata grid and atmo grid + ####################### + regrid_name = 'fireData2ff' + wkdir = './' + gridin_polygon, firein_polygon, gridin_xx, gridin_yy = regrid.get_polygon_from_grid(firedata.grid_e, firedata.grid_n, \ + flag_add_1extra_rowCol=True, firedata=firedata) + gridout_polygon, gridout_xx, gridout_yy = regrid.get_polygon_from_grid(bmap_grid_x2d, bmap_grid_y2d) + + in_outFF_grid_list, outFF_in_grid_list = regrid.grid2grid_pixel_match(regrid_name, \ + gridin_xx, gridin_yy, gridin_polygon, firein_polygon,\ + gridout_polygon, \ + wkdir, \ + flag_parallel=ronan_param['flag_parallel'], \ + flag_freshstart=ronan_param['flag_freshStart']) + + #regrid evaporation time if needed + if evaporationTime_out.size != dimensions_value_EVAP[0]*dimensions_value_EVAP[1]: + grid_x_ = .5*(gridout_xx[:-1,:-1]+gridout_xx[1:,1:]) + grid_y_ = .5*(gridout_yy[:-1,:-1]+gridout_yy[1:,1:]) + points = np.dstack( ( .5*(atmo_grid_x2d[:-1,:-1]+atmo_grid_x2d[1:,1:]).flatten(), .5*(atmo_grid_y2d[:-1,:-1]+atmo_grid_y2d[1:,1:]).flatten() ) )[0] + values = evaporationTime_out.flatten() + evaporationTime_out_regrid = interpolate.griddata(points, values, (grid_x_, grid_y_), method='nearest') + else: + evaporationTime_out_regrid = evaporationTime_out + + + #create variables + for key in nc_variable_bmap: + dimensions_name, dimensions_value = get_dimension_name_and_value(key,nc_variable_bmap, nco) + nco.variables[key][:] = np.zeros(dimensions_value) + + if key == 'arrival_time_of_front': + print('map arrival time to bmap:') + arrivalTime_out, _ = regrid.map_data(outFF_in_grid_list,dimensions_value[::-1],firedata.arrivalTime,flag='average') + idx = np.where(arrivalTime_out>=0) + arrivalTime_out[idx] -= evaporationTime_out_regrid[idx] # time when evaporation start + + #adjust arrival time to residence time map + #idx_ = np.where(residenceTime_out<=0) + #arrivalTime_out[idx_] = -999 + + if ronan_param['init_mode'] == 'fromOverHeadImg': + print('********** WARNING ***********: applying shift on arrivalTime') + idx = np.where(arrivalTime_out!=-999) + shifTime = 0 if 'ignitiontime' not in list(ronan_param.keys()) else \ + -1*(attribute.domain['t0']+arrivalTime_out[idx].min())+ronan_param['ignitiontime'] + arrivalTime_out[idx] = (arrivalTime_out[idx] + attribute.domain['t0'] + shifTime) + else: + shifTime = 0. + + nco.variables[key][:] = np.array(arrivalTime_out.T,dtype=np.float32).reshape(dimensions_value) + + elif key == 'cell_active': + nco.variables[key][:] = np.ones(dimensions_value,dtype=int) + + else: + nco.variables[key][:] = np.zeros(dimensions_value) + + if nc_attribute[key] is not None: + for attvar in nc_attribute[key]: + nco.variables[key]._attributes[attvar] = nc_attribute[key][attvar] + + + #close nc file + nco.close() + + print('') + print('Time Info:') + print('ref MNH time : ', attribute.domain['t0']) + print('shift time : ', shifTime) + print('min firedata arrival time: ', firedata.arrivalTime[np.where(firedata.arrivalTime>=0)].min()) + idx = np.where(arrivalTime_out!=-999) + print('nbre points = {:d}'.format(idx[0].shape[0])) + print('Final arrivalTime min={:.1f} max={:.1f}'.format(arrivalTime_out[idx].min(),arrivalTime_out[idx].max())) + print('') + #save shift in space and time to combine later the original 2D firescene and the atmospheric domain + atmosDomainReference = {'time': attribute.domain['t0'] + shifTime, + 'x': shiftx, + 'y': shifty, + 'z': shiftz, + 'flag_rotation': flag_rotation, + 'shiftTheta' : shiftTheta, + 'rotce' : rotce, + 'rotcn' : rotcn, + } + with open( 'Inputs/'+'atmosDomainReference.p', "wb" ) as f : + pickle.dump( atmosDomainReference,f ) + + #save "_out" array for use in runff + out = np.zeros_like(fre_f_out,dtype=np.dtype([('fre_f',float),('fre_s',float), + ('burningTime',float),('residenceTime',float)])) + out = out.view(np.recarray) + out.fre_f = fre_f_out + out.fre_s = fre_s_out + out.burningTime = burningTime_out + out.residenceTime = residenceTime_out + + np.save('Inputs/'+'out_{:s}.npy'.format(expr_name), out) + + #--------------------------------------------- + # set ForeFire simulation + #--------------------------------------------- + sys.exit() + if flag_run_ff: + ff = forefire.ForeFire() + + + #set param + ff.setString('ForeFireDataDirectory','Inputs') + ff.setString('fireOutputDirectory','ForeFire/Outputs') + ff.setInt('outputsUpdate',outputsUpdate) + + ff.setString('NetCDFfile', mydatanc.split('/')[-1]) + ff.setString('fluxNetCDFfile',mydatanc.split('/')[-1]) + ff.setString('fuelsTableFile','fuels.ff') + ff.setString('BMapFiles', mybmapnc.split('/')[-1]) + + ff.setDouble("spatialIncrement",spatialIncrement) + ff.setDouble("perimeterResolution",perimeterResolution) + ff.setDouble("minimalPropagativeFrontDepth",minimalPropagativeFrontDepth) + + if ronan_param['heatFluxModel'] == 'heatFluxNominal': + ff.setDouble("nominalHeatFlux", ff_param['nominalHeatFlux']) + if ronan_param['vaporFluxModel'] == 'vaporFluxNominal': + ff.setDouble("nominalVaporFlux", ff_param['nominalVaporFlux']) + if ronan_param['scalarFluxModel'] == 'scalarFluxNominal': + ff.setDouble("nominalScalarFlux", ff_param['nominalScalarFlux']) + + #ff.setDouble("nominalHeatFlux",nominalHeatFlux) + #ff.setDouble("nominalVaporFlux",nominalVaporFlux) + #ff.setDouble("burningDuration",burningDuration) + + ff.setDouble("bmapOutputUpdate",bmapOutputUpdate) + ff.setInt("defaultHeatType",0) + ff.setInt("defaultscalarType",2) + ff.setInt("defaultVaporType",1) + + #ff.setInt('bmapLayer',1) + ff.setDouble("InitTime",InitTime) + + #set domain + ff.setInt("atmoNX",nx) + ff.setInt("atmoNY",ny) + ff.execute("FireDomain[sw=(%f,%f,0.);ne=(%f,%f,0.);t=%f]"%(attribute.domain['SWx'],attribute.domain['SWy'],\ + attribute.domain['NEx'],attribute.domain['NEy'], \ + attribute.domain['t0'])) + + extentLocal= getDomainExtent(ff.execute("print[]").split("\n")[0]); + + + #set propagation model + ff.addLayer("propagation","TroisPourcent","propagationModel") + + #residenceTime = np.zeros((nx,ny,1)) + 60 + #ff.addScalarLayer("double","residenceTime",0 , 0, 0,extentLocal[1]-extentLocal[0], extentLocal[3]-extentLocal[2] , 0, residenceTime) + + fuelmap = ff.getDoubleArray("fuel").astype("int32") + ff.addIndexLayer("table","fuel", extentLocal[0], extentLocal[2],0, extentLocal[1]-extentLocal[0], extentLocal[3]-extentLocal[2], 0, fuelmap) + + print("resolution of bmap is ", ff.getString("bmapResolution")) + + #set fire line + if flag_extra_contour: + print('add an extra coutour') + ff.execute("\tFireFront[t={%f}]".format(attribute.domain['t0'])) + ff.execute("\t\tFireNode[loc=(090,100,0.);vel=(0.,0.,0.);t=%f]".format(attribute.domain['t0'])) + ff.execute("\t\tFireNode[loc=(090,300,0.);vel=(0.,0.,0.);t=%f]".format(attribute.domain['t0'])) + ff.execute("\t\tFireNode[loc=(100,300,0.);vel=(0.,0.,0.);t=%f]".format(attribute.domain['t0'])) + ff.execute("\t\tFireNode[loc=(100,100,0.);vel=(0.,0.,0.);t=%f]".format(attribute.domain['t0'])) + + sys.exit() + #--------------------------------------------- + # run ForeFire simulation + #--------------------------------------------- + plt.figure(1) + plt.ion() + pathes = [] + step = 50 + time_duration = 2000 + N_step = int(old_div(time_duration, step)) + 1 + #step = 10 + #N_step = 80/step + flux_out_ff_history = [] + fre_ff = np.zeros_like(fre_f_out) + h2o_ff = np.zeros_like(fre_f_out) + scalar_ff = np.zeros_like(fre_f_out) + chf_ff = np.zeros_like(fre_f_out) + for i in np.arange(1,N_step): + + ff_time = attribute.domain['t0'] + i*step + + print("goTo[t=%f]"%(ff_time), end=' ') + ff.execute("goTo[t=%f]"%(ff_time)) + + + #HeatFLux + Hflux2d = ff.getDoubleArray('heatFlux')[0,0,:,:] * 1.e-6 # MW/m2 (sensible heat flux) + flux_out_ff = (Hflux2d * atmo_dx * atmo_dy).sum() #MW + + frp_ff_add = old_div(ff_param['fractionRadiation_f'],(1.-ff_param['fractionRadiation_f'])) * Hflux2d * atmo_dx * atmo_dy # MW FRP + fre_ff_add = frp_ff_add * step # MJ FRE + + fre_ff += (fre_ff_add).T + chf_ff += (Hflux2d * step).T + + #Scalar Flux + Sflux2d = ff.getDoubleArray('scalarFlux')[0,0,:,:] + scalar_ff += (Sflux2d * atmo_dx * atmo_dy * step).T + + #Vapor Flux + Vflux2d = ff.getDoubleArray('vaporFlux')[0,0,:,:] + h2o_ff += (Vflux2d * atmo_dx * atmo_dy * step).T + + ''' + fig = plt.figure(figsize=(15,8)) + ax = plt.subplot(131) + ax.imshow(Sflux2d.T,origin='lower',interpolation='nearest') + ax.set_title('Scalar') + bx = plt.subplot(132) + bx.imshow(Hflux2d.T,origin='lower',interpolation='nearest') + bx.set_title('Heat') + cx = plt.subplot(133) + cx.imshow(Vflux2d.T,origin='lower',interpolation='nearest') + cx.set_title('Vapor') + plt.show() + ''' + + if ronan_param['init_mode'] == 'ideal': + idx_fire_line_f = np.where( (np.array(time_start_fireLine) <= ff_time) & + (np.array(time_start_fireLine)+np.array(residence_time) >= ff_time) ) + idx_fire_line_s = np.where( (np.array(time_start_fireLine)+np.array(residence_time) < ff_time) & + (np.array(time_start_fireLine)+np.array(burning_time) >= ff_time) ) + flux_expected = 0 + for idx_ in idx_fire_line_f[0]: + flux_expected += depth_fireLine[idx_]*length_fireLine[idx_]*nominalHeatFlux_f[idx_] * 1.e-6 + for idx_ in idx_fire_line_s[0]: + flux_expected += depth_fireLine[idx_]*length_fireLine[idx_]*nominalHeatFlux_s[idx_] * 1.e-6 + + pathes += printToPathe( ff.execute("print[]")) + if flux_expected != 0: + flux_out_ff_history.append(old_div(flux_out_ff,flux_expected)) + else: + if flux_out_ff < 1.e-6: + flux_out_ff_history.append(1) + else: + flux_out_ff_history.append(0) + + print('| flux out of ff = ', flux_out_ff, '| expected = ', flux_expected, '| ratio = ', flux_out_ff_history[-1]) + + if ronan_param['init_mode'] == 'fromOverHeadImg': + print('| flux out of ff = ', flux_out_ff , 1.e-6*fre_ff_add.max()) + + + plt.clf() + ax = plt.subplot(111) + ax.imshow(np.ma.masked_where(Hflux2d<=0,Hflux2d).T,origin='lower',interpolation='nearest',vmin=0.001,vmax=1.e-1) + ax.set_title('t={:.1f}'.format(ff_time)) + #plt.show() + plt.draw() + plt.pause(.001) + + if np.isnan(fre_ff.sum()): + pdb.set_trace() + + + plt.ioff() + plt.close() + + + print('conservation FRE') + print(fre_ff.sum()) + print((fre_f_out+fre_s_out)[np.where(fre_f_out+fre_s_out>0)].sum()*1.e-6) + print((firedata.fre_f+firedata.fre_s)[np.where((firedata.fre_f+firedata.fre_s>0))].sum()) + print((old_div(ff_param['fractionRadiation_f'],(1.-ff_param['fractionRadiation_f'])))*(sensible_heat_flux_f * residenceTime_out * atmo_dx * atmo_dy).sum() * 1.e-6 + \ + (old_div(ff_param['fractionRadiation_s'],(1.-ff_param['fractionRadiation_s'])))*(sensible_heat_flux_s * (burningTime_out-residenceTime_out) * atmo_dx * atmo_dy).sum() * 1.e-6) + + print('') + print('conservation Vapor') + print(h2o_ff.sum()) + print((mass_comsumption_out * 1.e-2 * moisture_out ).sum()) + + print('') + print('conservation Scalar') + print(scalar_ff.sum()) + print(ff_param['conversionFactor'] * 1.e-6 * ( 0.01*fre_f_out + .02*fre_s_out).sum()) + + print('') + print('mean nominal flux') + fuelparams = pandas.read_csv('./Inputs/fuels.ff', sep=';') + firedurationmap = np.zeros_like(fuelmap) + for fuelidx,tau0,sd in zip(fuelparams['Index'], fuelparams['Tau0'], fuelparams['sd']): + idx_ = np.where(fuelmap==fuelidx) + firedurationmap[idx_] = tau0/sd + fireduration_mean = firedurationmap[0,0,:,:][np.where(chf_ff>0)].mean() + + print('fireduration_mean', fireduration_mean) + + meanNominalHeatFlux = chf_ff.sum()/(np.where(chf_ff>0)[0].shape[0]) / fireduration_mean + meanNominalVaporFlux = h2o_ff.sum()/(np.where(chf_ff>0)[0].shape[0]) / (atmo_dx * atmo_dy) / fireduration_mean + meanNominalScalarFlux = scalar_ff.sum()/(np.where(chf_ff>0)[0].shape[0]) / (atmo_dx * atmo_dy) / fireduration_mean + + print('meanNominalHeatFlux (W/m2)', meanNominalHeatFlux*1.e6) + print('meanNominalVaporFlux', meanNominalVaporFlux) + print('meanNominalScalarFlux', meanNominalScalarFlux) + + lines=[] + lines.append( 'meanNominalHeatFlux (W/m2) = {:.8f}\n'.format( meanNominalHeatFlux*1.e6) ) + lines.append( 'meanNominalVaporFlux (kg/m2) = {:.8f}\n'.format( meanNominalVaporFlux) ) + lines.append( 'meanNominalScalarFlux (-/m2) = {:.8f}\n'.format( meanNominalScalarFlux) ) + + with open('equivalentNominalFlux.txt','w') as f: + f.writelines(lines) + + + if ronan_param['init_mode'] == 'ideal': + ax = plt.subplot(111) + ax.plot(step*np.arange(1,N_step),np.array(flux_out_ff_history)) + ax.set_ylabel('ratio: heat_ff / heat_expected') + ax.set_xlabel('time (s)') + plt.show() + + if ronan_param['init_mode'] == 'fromOverHeadImg': + fig = plt.figure(figsize=(12,6)) + ax = plt.subplot(131) + im = ax.imshow(np.ma.masked_where(fre_ff<=0,fre_ff).T,origin='lower',interpolation='nearest') + ax.set_title('fre_ff = \n rad_frac x flux_conv_FF') + divider = make_axes_locatable(ax) + cbaxes = divider.append_axes("bottom", size="5%", pad=0.05) + cbar = fig.colorbar(im ,cax = cbaxes,orientation='horizontal') + cbar.set_label('FRE (MW)',labelpad=10) + + ax = plt.subplot(132) + im = ax.imshow(np.ma.masked_where( (fre_f_out+fre_s_out)<=0,1.e-6*(fre_f_out+fre_s_out)).T,origin='lower',interpolation='nearest') + ax.set_title('fre_obs \n from Obs at resolution of MNH') + divider = make_axes_locatable(ax) + cbaxes = divider.append_axes("bottom", size="5%", pad=0.05) + cbar = fig.colorbar(im ,cax = cbaxes,orientation='horizontal') + cbar.set_label('FRE (MW)',labelpad=10) + + ax = plt.subplot(133) + ratio = np.zeros_like(fre_f_out) + idx = np.where(fre_f_out+fre_s_out>0) + ratio[idx] = old_div(fre_ff[idx],((fre_f_out+fre_s_out)[idx]*1.e-6)) + im = ax.imshow(np.ma.masked_where((fre_f_out+fre_s_out)<=0,ratio).T,origin='lower',interpolation='nearest') + ax.set_title('ratio fre_ff/fre_obs') + divider = make_axes_locatable(ax) + cbaxes = divider.append_axes("bottom", size="5%", pad=0.05) + cbar = fig.colorbar(im ,cax = cbaxes,orientation='horizontal') + cbar.set_label('ratio (-)',labelpad=10) + + fig.savefig(expr_name+'_fre_conservation.png') + plt.show() + + + sys.exit() + #--------------------------------------------- + # plot ForeFire perimeter + #--------------------------------------------- + fig, ax = plt.subplots() + + #tab = np.transpose(ff.getDoubleArray("BMap"))[0] + #CS = ax.imshow(tab, origin='lower', cmap=plt.cm.gray, interpolation='nearest',\ + # extent=(0,tab.shape[1]*np.float(ff.getString("bmapResolution")),0,tab.shape[0]*np.float(ff.getString("bmapResolution")))) + #cbar = plt.colorbar(CS) + #cbar.ax.set_ylabel('v') + cmap = mpl.cm.get_cmap('jet') + for i, path in enumerate(pathes): + rgba = cmap(old_div(1.*i,(N_step-1))) + patch = mpatches.PathPatch(path,edgecolor=rgba, facecolor='none', alpha=1) + ax.add_patch(patch) + + ax.grid() + ax.axis('equal') + ax.set_xlim(attribute.domain['SWx']+atmo_dx,attribute.domain['NEx']-atmo_dx) + ax.set_ylim(attribute.domain['SWy']+atmo_dy,attribute.domain['NEy']-atmo_dy) + plt.show() + + diff --git a/tools/preprocessing/create_data_tools.py b/tools/preprocessing/create_data_tools.py new file mode 100644 index 00000000..b6a5f66e --- /dev/null +++ b/tools/preprocessing/create_data_tools.py @@ -0,0 +1,146 @@ +from __future__ import print_function +from __future__ import division +from builtins import zip +from builtins import range +from past.utils import old_div +import sys +import os +import numpy as np +import matplotlib as mpl +mpl.use('Qt5Agg') +import matplotlib.pyplot as plt +import matplotlib.path as mpath +import matplotlib.patches as mpatches +import pdb +#import imp +from scipy import io +import glob +import netCDF4 +import f90nml +import datetime +import socket +import warnings +warnings.filterwarnings("ignore",".*GUI is implemented.*") +from mpl_toolkits.axes_grid1 import make_axes_locatable +import pickle +import argparse +from osgeo import gdal,osr,ogr +from scipy import interpolate +import importlib +import pandas +import math +import pyproj + +########################## +def load_crs_from_prj(prj_file_path): + # Check if the .prj file exists + if not os.path.exists(prj_file_path): + raise FileNotFoundError(f"The file {prj_file_path} does not exist.") + + # Read the WKT from the .prj file + with open(prj_file_path, 'r') as prj_file: + wkt = prj_file.read().strip() + + # Create a CRS object from the WKT + crs = pyproj.CRS(wkt) + + return crs + + +################################################## +def ensure_dir(f): + d = os.path.dirname(f) + if not os.path.exists(d): + os.makedirs(d) + + +################################################## +def read_MNH_domain( expr_name, dtoutput_mnh, dir_mnh = '../02_mnh/', domainNumber=None, specificInputMNHFile=None): + + if specificInputMNHFile == None: + expr_numb = '.{:d}.'.format(domainNumber) if (domainNumber is not None) else '' + outfiles_mnh = sorted(list(set(glob.glob(dir_mnh+expr_name[:5]+expr_numb+'*.nc'))-\ + set(glob.glob(dir_mnh+expr_name[:5]+expr_numb+'*.OUT.*.nc')))+\ + glob.glob(dir_mnh+expr_name[:5]+expr_numb+'*.nc4')) + #print(outfiles_mnh) + if len(outfiles_mnh) == 0: + print('missing mnh files') + pdb.set_trace() + outfiles_mnh_inputTime = outfiles_mnh[0] + outfiles_mnh_inputGrid = outfiles_mnh[1] + + else: + outfiles_mnh_inputGrid = specificInputMNHFile + outfiles_mnh_inputTime = specificInputMNHFile + + print('get time ref from :', outfiles_mnh_inputTime) + print('get grid from :', outfiles_mnh_inputGrid) + + + #read geometry of the MNH domain + nc = netCDF4.Dataset(outfiles_mnh_inputGrid,'r') # get the 001 file + + MNH_properties= {} + + MNH_properties['lat00'] = np.array(nc.variables['LAT'])[0,0] + MNH_properties['lon00'] = np.array(nc.variables['LON'])[0,0] + MNH_properties['lat2d'] = np.array(nc.variables['LAT']) + MNH_properties['lon2d'] = np.array(nc.variables['LON']) + + DeltaY = nc.variables['YHAT'][1]-nc.variables['YHAT'][0] + DeltaX = nc.variables['XHAT'][1]-nc.variables['XHAT'][0] + + MNH_properties['nx'] = len(nc.dimensions[u'ni']) + MNH_properties['ny'] = len(nc.dimensions[u'nj']) + MNH_properties['nz'] = len(nc.dimensions[u'level_w']) + + MNH_properties['SWx'] = nc.variables['XHAT'][0]+0. + MNH_properties['SWy'] = nc.variables['YHAT'][0]+0. + MNH_properties['SWz'] = 0 + MNH_properties['NEx'] = nc.variables['XHAT'][-1]+DeltaX + MNH_properties['NEy'] = nc.variables['YHAT'][-1]+DeltaY + MNH_properties['NEz'] = 0 + + MNH_properties['Lx'] = float(nc.variables['XHAT'][-1]+DeltaX-nc.variables['XHAT'][0]) + MNH_properties['Ly'] = float(nc.variables['YHAT'][-1]+DeltaY-nc.variables['YHAT'][0]) + MNH_properties['Lz'] = 0 + nc.close() + + try: + nc = netCDF4.Dataset(outfiles_mnh_inputTime,'r') + MNH_properties['date'] = '_'.join(nc.variables['DTCUR'].units.split(' ')[2:4]) + MNH_properties['t0'] = float(nc.variables['DTCUR'][0]) + except: + nc = netCDF4.Dataset(outfiles_mnh_inputGrid,'r') + MNH_properties['date'] = '_'.join(nc.variables['DTCUR'].units.split(' ')[2:4]) + MNH_properties['t0'] = float(nc.variables['DTCUR'][0]-dtoutput_mnh) + + MNH_properties['Lt'] = np.Inf + + + nc.close() + + + return MNH_properties + +####################################### +def offset_in_degrees(distance_m, lat): + """ + Convert a distance in meters to an offset in degrees for both latitude and longitude. + + Parameters: + distance_m (float): Distance in meters. + lat (float): Latitude in degrees. + + Returns: + lat_offset (float): Latitude offset in degrees. + lon_offset (float): Longitude offset in degrees. + """ + # Earth's radius approximation for latitude degrees + lat_offset = distance_m / 111320 + + # Longitude offset depends on latitude + lon_offset = distance_m / (111320 * math.cos(math.radians(lat))) + + return lat_offset, lon_offset + diff --git a/tools/preprocessing/ffToGeo.py b/tools/preprocessing/ffToGeo.py new file mode 100755 index 00000000..643da170 --- /dev/null +++ b/tools/preprocessing/ffToGeo.py @@ -0,0 +1,200 @@ +from __future__ import print_function +from __future__ import division +from builtins import range +from builtins import object +from past.utils import old_div +import os +import sys +from pyproj import Proj, transform +import json +import datetime +import pdb +import geopandas as gpd +import shapely +import pyproj +import geojson +from functools import partial +import copy + +#http://dev.firecaster.valabre.fr/forefireAPI.php?command=ignition&longitude=42.348570491488999&date_simulation=2020-01-04T09%3A37%3A09Z&date_eclosion=2020-02-04T09%3A37%3A09Z&latitude=9.017890207687699&pathdate=2020-01-07_T10_47_15Z + +#http://dev.firecaster.valabre.fr/forefireAPI.php?command=sequence&duree_h=4&vent_direction_degre=20&vent_vitesse_kmh=50&temperature=10&hygrometrie=10&pathdate=2020-01-07_T10_47_15Z + + +# $duree_h=$_REQUEST['duree_h']; +# 11 $vent_direction_degre=$_REQUEST['vent_direction_degre']; +# 12 $vent_vitesse_kmh=$_REQUEST['vent_vitesse_kmh']; +# 13 $hygrometrie=$_REQUEST['hygrometrie']; +# 14 $temperature=$_REQUEST['temperature']; +# 15 $pathdate=$_REQUEST['pathdate']; + +################################# +def projGeometry2WGS84(feature,projin): + + # the shapely projection function + datum_wgs84 = pyproj.Proj(init='EPSG:4326') + projection_wm_func = partial(pyproj.transform, projin, datum_wgs84) + + # project GPS lat/lon coordinates to web mercator for each polygon in the geojson + polygon = shapely.geometry.shape(feature['geometry']) + return [point for polygon in shapely.ops.transform(projection_wm_func, polygon) for point in polygon.exterior.coords[:-1]] + +################################# +class FFGEO(object): + + def __init__(self, name, proj): + self.json = {"type": "FeatureCollection", + "name": "Front name:%s"%name, + "features": [] } + self.proj = proj + self.fronts = [] + + def addFront(self,ff): + + def parse(dataString): + + def proj(point,inProj = None, outProj = None): + if inProj is None or outProj is None: + return point + return transform(inProj,outProj,point[0],point[1]) + + def isPoint(element): + if len(element) == 2: + if isinstance(element[0],float) and isinstance(element[1],float) : + return True + return False + + def getLocationFromLine(line,pattern="loc=("): + llv = line.split(pattern) + if len(llv) < 2: + return None + llr = llv[1].split(","); + if len(llr) < 3: + return None + return (float(llr[0]),float(llr[1])) + ''' + def printToPolygons(self, linePrinted, inProj = None, outProj = None, level=1): + if level > 8: + return + + #pointsMap + + fronts = linePrinted.split("\n%sFireFront"%('\t'*level)) + + #out = [] + + if len(fronts)>0: + nodes = fronts[0].split("FireNode") + if len(nodes) > 1: + pointsMap = [] + for node in nodes[1:]: + pointsMap.append(proj(getLocationFromLine(node),inProj, outProj)) + pointsMap.append(proj(getLocationFromLine(nodes[1]),inProj, outProj)) + self.fronts.append([level, pointsMap]) + print level + + for subline in fronts[1:]: + self.fronts.append(printToPolygons(self,subline,inProj, outProj, level+1)) + ''' + def printToPolygons(linePrinted, inProj = None, outProj = None, level=1): + if level > 8: + return + fronts = linePrinted.split("\n%sFireFront"%('\t'*level)) + + out = [] + if len(fronts)>0: + nodes = fronts[0].split("FireNode") + if len(nodes) > 1: + pointsMap = [] + for node in nodes[1:]: + pointsMap.append(proj(getLocationFromLine(node),inProj, outProj)) + pointsMap.append(proj(getLocationFromLine(nodes[1]),inProj, outProj)) + out.append([level,pointsMap]) + + for subline in fronts[1:]: + out.append(printToPolygons(subline,inProj, outProj, level+1)) + + return out + + def splitPoly(self, treePoly, dadLevel=999, featId_here=0): + + sousFront = [] + for poly in treePoly[0]: + + if (len(poly)==2) & (isinstance(poly[0], int) ): + levelHere = old_div(poly[0],2) + if levelHere > dadLevel: + self.out.append([]) + self.featId += 1 + featId_here = self.featId + #print levelHere, dadLevel, self.featId, featId_here + self.out[featId_here].append(poly[1]) + dadLevel=levelHere + + else: + splitPoly(self, [poly], dadLevel=dadLevel, featId_here=featId_here) + + + treePoly = printToPolygons(dataString) + self.out = [[]] + self.featId = 0 + splitPoly(self, treePoly) + + return self.out + + + fronts_polygons = parse(ff.execute("print[]")) + + feature_template = { "type": "Feature", + "properties": { "area": 999, + "date": 777, + "fill":"#ffffff", + "fill-opacity":0, + "stroke": "#555555", + "stroke-opacity": 1, + "stroke-width":"2", + }, + "geometry": { "type": "MultiPolygon", + }, + } + + for front_coords in fronts_polygons: + + feature = copy.deepcopy(feature_template) + try: + feature['properties']['area'] = shapely.geometry.Polygon(front_coords[0]).area*1.e-4 + feature['properties']['date'] =ff.getString("ISOdate") + feature['geometry']['coordinates'] = [front_coords] + except: + pdb.set_trace() + self.json['features'].append(feature) + + #re-compute geopanda dataframe from geojson + try: + self.pd = gpd.GeoDataFrame.from_features(self.json["features"]) + self.pd['date'] = [datetime.datetime.strptime(date, '%Y-%m-%dT%H:%M:%SZ') for date in self.pd['date'] ] + self.pd['timetofront'] = [ (self.pd['date'].iloc[-1]-date).total_seconds() for date in self.pd['date'] ] + self.pd = self.pd.sort_values(by='timetofront', ascending=True) + self.pd.loc[self.pd['timetofront']>4.*3600, 'timetofront'] = 4.*3600 + except: + pdb.set_trace() + + def dumgeojson(self,filename): + + # project geojson to 4326 + ffgeojson = dict(self.json) + for i in range(len(ffgeojson['features'])): + ffgeojson['features'][i]['geometry']['coordinates'] = [[projGeometry2WGS84(ffgeojson['features'][i],self.proj)]] + + #and save + with open(filename, 'w') as f: + geojson.dump(ffgeojson, f) + + +if __name__ == '__main__': + + ff = forefire.PLibForeFire() + ffFronts = ffgeojson('test') + ffFronts.addFront(ff) + + print(ffFronts.json) diff --git a/tools/preprocessing/forefire_helper.py b/tools/preprocessing/forefire_helper.py new file mode 100644 index 00000000..a4fd6580 --- /dev/null +++ b/tools/preprocessing/forefire_helper.py @@ -0,0 +1,512 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.path as mpath +import matplotlib.patches as mpatches +import matplotlib.colors as mcolors +import math +import xarray as xr +import pyforefire as forefire +from datetime import datetime +import time + +def standardRothermelFuelTable(): + return """Index;Rhod;Rhol;Md;Ml;sd;sl;e;Sigmad;Sigmal;stoch;RhoA;Ta;Tau0;Deltah;DeltaH;Cp;Cpa;Ti;X0;r00;Blai;me +111;500.0;500.0;0.15;0.5;2400.0;5700.0;0;0.0;0;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +112;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +121;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +122;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +123;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +124;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +131;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +132;720.0;720.0;0.08;1;5544.0;4766.0;0;0.89;1.79;8.3;1.0;300.0;75590.0;2300000.0;1.6E7;2000.0;1100.0;600.0;0.3;2.5E-5;4.0;0.3 +133;720.0;720.0;0.08;1;5544.0;4766.0;0;0.89;1.79;8.3;1.0;300.0;75590.0;2300000.0;1.6E7;2000.0;1100.0;600.0;0.3;2.5E-5;4.0;0.3 +141;512.0;512.0;0.08;1;6562;5905;0.46;0.22;0.25;8.3;1.0;300.0;75590.0;2300000.0;1.86E7;1800.0;1000.0;600.0;0.3;2.5E-05;4.0;0.3 +142;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +211;500.0;500.0;0.13;0.5;2400.0;5700.0;1;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +212;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +213;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +221;500.0;500.0;0.13;0.5;2400.0;5700.0;0.5;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +222;500.0;500.0;0.13;0.5;2400.0;5700.0;2;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +223;500.0;500.0;0.13;0.5;2400.0;5700.0;2;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +231;500.0;500.0;0.13;0.5;2400.0;5700.0;2;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +241;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +242;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +243;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +244;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +311;512.0;512.0;0.13;0.5;4922.0;2460.0;0.3;0.9;0.67;8.3;1.0;300.0;70000.0;2300000.0;1.86E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +312;512.0;512.0;0.13;0.5;4922.0;2460.0;0.3;0.9;0.67;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +313;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +321;512.0;512.0;0.13;0.5;5905.0;5250.0;0.46;0.54;0.11;8.3;1.0;300.0;70000.0;2300000.0;1.86E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +322;512.0;512.0;0.13;0.5;2460.0;5250.0;1.2;0.8;0.65;8.3;1.0;300.0;70000.0;2300000.0;1.86E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +323;512.0;512.0;0.13;0.5;2460.0;5250.0;1.8;0.8;0.65;8.3;1.0;300.0;70000.0;2300000.0;1.86E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +324;512.0;512.0;0.13;0.5;2460.0;5250.0;1.8;0.8;0.65;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +331;500.0;500.0;1.6;2;2400.0;5700.0;0;10;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +332;500.0;500.0;10;10;2400.0;5700.0;0;10;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +333;512.0;512.0;0.08;0.5;6560.0;5900.0;0.3;0.2;0.05;8.3;1.0;300.0;70000.0;2300000.0;1.86E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +334;512.0;512.0;0.08;0.5;6560.0;5900.0;0.3;0.2;0.05;8.3;1.0;300.0;70000.0;2300000.0;1.86E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +335;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +411;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +412;500.0;500.0;0.13;0.5;2400.0;5700.0;0.1;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +421;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +422;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +423;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +511;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +512;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +521;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +522;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +523;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3""" + + +def extendedRothermelFuelTable(): + return """Index;Rhod;Rhol;Md;Ml;sd;sl;e;Sigmad;Sigmal;stoch;RhoA;Ta;Tau0;Deltah;DeltaH;Cp;Cpa;Ti;X0;r00;Blai;me +0;563.0;522.0;0.1;1.0;6099.0;7273.0;0;0.764;0.352;8.3;1.0;300;70000;18169000.0;18167000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +1;563.0;522.0;0.1;1.0;6099.0;7273.0;0.24;0.764;0.352;8.3;1.0;300;70000;18169000.0;18167000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +2;614.0;613.0;0.1;1.0;4287.0;5738.0;0.4;1.378;0.174;8.3;1.0;300;70000;18727000.0;18727000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +3;614.0;613.0;0.1;1.0;4287.0;5738.0;0.4;1.378;0.174;8.3;1.0;300;70000;18727000.0;18727000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +4;613.0;538.0;0.1;1.0;4357.0;6524.0;0.19;1.286;0.085;8.3;1.0;300;70000;18677000.0;18677000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +5;626.0;600.0;0.1;1.0;4325.0;5844.0;0.6;1.393;0.201;8.3;1.0;300;70000;18802000.0;18802000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +6;562.0;474.0;0.1;1.0;6740.0;8195.0;0.57;1.326;0.166;8.3;1.0;300;70000;18941000.0;18941000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +7;658.0;651.0;0.1;1.0;4734.0;5733.0;0.15;1.415;0.541;8.3;1.0;300;70000;18472000.0;18466000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +8;446.0;513.0;0.1;1.0;7792.0;9072.0;0.78;0.492;0.023;8.3;1.0;300;70000;18587000.0;18587000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +9;467.0;543.0;0.1;1.0;6115.0;7224.0;0.285;0.855;0.174;8.3;1.0;300;70000;18474000.0;18474000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +10;674.0;612.0;0.1;1.0;4801.0;5928.0;0.45;1.525;0.709;8.3;1.0;300;70000;18280000.0;18277000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +11;653.0;582.0;0.1;1.0;4753.0;6569.0;0.75;1.096;1.105;8.3;1.0;300;70000;18226000.0;18221000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +12;596.0;586.0;0.1;1.0;3688.0;5551.0;0.475;1.346;0.077;8.3;1.0;300;70000;19050000.0;19050000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +13;438.0;488.0;0.1;1.0;7274.0;8453.0;0.38;1.053;0.321;8.3;1.0;300;70000;17842000.0;17842000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +14;634.0;570.0;0.1;1.0;5518.0;7704.0;0.2;1.293;0.402;8.3;1.0;300;70000;18507000.0;18507000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +15;667.0;636.0;0.1;1.0;5162.0;7145.0;0.3;1.457;0.519;8.3;1.0;300;70000;18851000.0;18851000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +16;563.0;531.0;0.1;1.0;4694.0;7370.0;0.285;1.18;0.143;8.3;1.0;300;70000;19235000.0;19235000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +17;599.0;531.0;0.1;1.0;5813.0;7370.0;0.285;1.18;0.143;8.3;1.0;300;70000;19041000.0;19041000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +18;565.0;531.0;0.1;1.0;4644.0;7370.0;0.285;1.18;0.143;8.3;1.0;300;70000;19169000.0;19169000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +19;565.0;531.0;0.1;1.0;4644.0;7370.0;0.285;1.18;0.143;8.3;1.0;300;70000;19169000.0;19169000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +20;678.0;573.0;0.1;1.0;6064.0;5456.0;0.475;1.604;0.241;8.3;1.0;300;70000;18992000.0;18992000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +21;597.0;442.0;0.1;1.0;7975.0;10000.0;0.09;0.566;0.008;8.3;1.0;300;70000;18887000.0;18887000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +22;551.0;545.0;0.1;1.0;5360.0;7065.0;1.0;1.283;0.261;8.3;1.0;300;70000;19021000.0;19021000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +23;565.0;531.0;0.1;1.0;4644.0;7370.0;0.285;1.18;0.143;8.3;1.0;300;70000;19169000.0;19169000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +24;659.0;609.0;0.1;1.0;6691.0;6974.0;0.2;0.923;0.469;8.3;1.0;300;70000;18920000.0;18920000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +25;491.0;482.0;0.1;1.0;5732.0;7599.0;0.77;1.287;0.247;8.3;1.0;300;70000;18802000.0;18802000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +26;597.0;442.0;0.1;1.0;7975.0;10000.0;0.09;0.566;0.008;8.3;1.0;300;70000;18887000.0;18887000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +27;590.0;471.0;0.1;1.0;7687.0;8345.0;0.475;0.259;0.139;8.3;1.0;300;70000;18715000.0;18715000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +28;565.0;491.0;0.1;1.0;4522.0;6577.0;0.8;1.343;0.125;8.3;1.0;300;70000;18866000.0;18864000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +29;559.0;481.0;0.1;1.0;4488.0;7793.0;0.38;1.139;0.106;8.3;1.0;300;70000;19207000.0;19207000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +30;533.0;527.0;0.1;1.0;3886.0;4358.0;0.42;1.0;0.251;8.3;1.0;300;70000;18508000.0;18499000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +31;625.0;601.0;0.1;1.0;5210.0;5747.0;0.475;1.188;1.006;8.3;1.0;300;70000;18922000.0;18919000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +32;625.0;601.0;0.1;1.0;5210.0;5747.0;0.475;1.188;1.006;8.3;1.0;300;70000;18922000.0;18919000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +33;500.0;442.0;0.1;1.0;8359.0;10000.0;0.27;0.469;0.0;8.3;1.0;300;70000;17129000.0;17129000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +34;621.0;600.0;0.1;1.0;5314.0;5947.0;0.665;1.103;0.811;8.3;1.0;300;70000;18853000.0;18847000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +35;502.0;442.0;0.1;1.0;8135.0;10000.0;0.3;0.482;0.0;8.3;1.0;300;70000;16931000.0;16931000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +36;592.0;547.0;0.1;1.0;5819.0;7229.0;0.57;1.209;0.606;8.3;1.0;300;70000;18426000.0;18422000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +37;502.0;442.0;0.1;1.0;8888.0;10000.0;0.3;0.482;0.0;8.3;1.0;300;70000;17103000.0;17103000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +38;624.0;603.0;0.1;1.0;5341.0;5776.0;1.2;1.533;1.387;8.3;1.0;300;70000;18912000.0;18907000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +39;618.0;584.0;0.1;1.0;6104.0;6996.0;0.76;0.549;0.499;8.3;1.0;300;70000;19075000.0;19029000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +40;618.0;584.0;0.1;1.0;6104.0;6996.0;0.76;0.549;0.499;8.3;1.0;300;70000;19075000.0;19029000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +41;618.0;584.0;0.1;1.0;6104.0;6996.0;0.76;0.549;0.499;8.3;1.0;300;70000;19075000.0;19029000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +42;762.0;761.0;0.1;1.0;6216.0;6985.0;0.728;0.809;1.069;8.3;1.0;300;70000;19130000.0;18581000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +43;644.0;595.0;0.1;1.0;6667.0;6781.0;0.95;0.426;0.624;8.3;1.0;300;70000;19036000.0;18986000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +44;577.0;550.0;0.1;1.0;5307.0;6592.0;0.665;0.836;0.344;8.3;1.0;300;70000;19046000.0;19021000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +45;577.0;550.0;0.1;1.0;5307.0;6592.0;0.665;0.217;0.344;8.3;1.0;300;70000;19046000.0;19021000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +46;544.0;486.0;0.1;1.0;5101.0;7518.0;0.57;0.109;0.159;8.3;1.0;300;70000;19295000.0;19295000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +47;623.0;651.0;0.1;1.0;4102.0;5421.0;0.2;1.349;0.12;8.3;1.0;300;70000;18632000.0;18632000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +48;615.0;606.0;0.1;1.0;4545.0;4814.0;0.6;1.009;1.126;8.3;1.0;300;70000;18771000.0;18760000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +49;591.0;584.0;0.1;1.0;4441.0;4857.0;1.05;1.576;0.818;8.3;1.0;300;70000;19322000.0;19322000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +50;661.0;617.0;0.1;1.0;6252.0;7700.0;0.18;0.321;0.647;8.3;1.0;300;70000;17630000.0;17529000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +51;831.0;792.0;0.1;1.0;6758.0;7320.0;0.2;0.407;1.134;8.3;1.0;300;70000;19144000.0;18486000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +52;442.0;442.0;0.1;1.0;10000.0;10000.0;0.3;0.04;0.092;8.3;1.0;300;70000;16298000.0;16298000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +53;554.0;442.0;0.1;1.0;7236.0;10000.0;0.21;0.128;0.065;8.3;1.0;300;70000;17428000.0;17428000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +54;436.0;464.0;0.1;1.0;6435.0;8759.0;0.285;0.78;0.102;8.3;1.0;300;70000;18396000.0;18396000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +55;554.0;442.0;0.1;1.0;7236.0;10000.0;0.14;0.118;0.043;8.3;1.0;300;70000;17428000.0;17428000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +62;563.0;522.0;0.1;1.0;6099.0;7273.0;0;0.764;0.352;8.3;1.0;300;70000;18169000.0;18167000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +75;614.0;613.0;0.25;1.0;4287.0;5738.0;0.1;1.378;0.174;8.3;1.0;300;70000;18727000.0;18727000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +73;613.0;538.0;0.1;1.0;4357.0;6524.0;0.19;1.286;0.085;8.3;1.0;300;70000;18677000.0;18677000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +82;614.0;613.0;0.1;1.0;4287.0;5738.0;0.4;1.378;0.174;8.3;1.0;300;70000;18727000.0;18727000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +83;563.0;522.0;0.05;1.0;6099.0;7273.0;0.4;0.764;0.352;8.3;1.0;300;70000;18169000.0;18167000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +102;626.0;600.0;0.1;1.0;4325.0;5844.0;0.6;1.393;0.201;8.3;1.0;300;70000;18802000.0;18802000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +103;626.0;600.0;0.1;1.0;4325.0;5844.0;0.6;1.393;0.201;8.3;1.0;300;70000;18802000.0;18802000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +104;626.0;600.0;0.1;1.0;4325.0;5844.0;0.6;1.393;0.201;8.3;1.0;300;70000;18802000.0;18802000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +105;563.0;522.0;0.1;1.0;6099.0;7273.0;0;0.764;0.352;8.3;1.0;300;70000;18169000.0;18167000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +106;563.0;522.0;0.1;1.0;6099.0;7273.0;0;0.764;0.352;8.3;1.0;300;70000;18169000.0;18167000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +111;500.0;500.0;0.15;0.5;2400.0;5700.0;0;0.0;0;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +112;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +121;563.0;522.0;0.1;1.0;6099.0;7273.0;0;0.764;0.352;8.3;1.0;300;70000;18169000.0;18167000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +123;563.0;522.0;0.1;1.0;6099.0;7273.0;0;0.764;0.352;8.3;1.0;300;70000;18169000.0;18167000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +124;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +131;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +132;720.0;720.0;0.08;1;5544.0;4766.0;0;0.89;1.79;8.3;1.0;300.0;75590.0;2300000.0;1.6E7;2000.0;1100.0;600.0;0.3;2.5E-5;4.0;0.3 +133;720.0;720.0;0.08;1;5544.0;4766.0;0;0.89;1.79;8.3;1.0;300.0;75590.0;2300000.0;1.6E7;2000.0;1100.0;600.0;0.3;2.5E-5;4.0;0.3 +141;720.0;720.0;0.08;1;5544.0;4766.0;0.5;0.89;1.79;8.3;1.0;300.0;75590.0;2300000.0;1.6E7;2000.0;1100.0;600.0;0.3;2.5E-5;4.0;0.3 +142;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +162;563.0;522.0;0.1;1.0;6099.0;7273.0;0;0.764;0.352;8.3;1.0;300;70000;18169000.0;18167000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +211;500.0;500.0;0.13;0.5;2400.0;5700.0;1;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +212;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +213;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +221;500.0;500.0;0.13;0.5;2400.0;5700.0;0.5;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +222;500.0;500.0;0.13;0.5;2400.0;5700.0;2;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +223;500.0;500.0;0.13;0.5;2400.0;5700.0;2;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +231;500.0;500.0;0.13;0.5;2400.0;5700.0;2;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +241;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +242;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +243;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +244;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +255;563.0;522.0;0.1;1.0;6099.0;7273.0;0;0.764;0.352;8.3;1.0;300;70000;18169000.0;18167000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3 +311;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +312;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +313;500.0;500.0;0.13;0.5;2400.0;5700.0;1.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +321;500.0;500.0;0.13;0.5;2400.0;5700.0;0.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +322;500.0;500.0;0.13;0.5;2400.0;5700.0;2.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +323;500.0;500.0;0.08;0.4;6000.0;5000.0;2;0.4;1.8;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +324;500.0;500.0;0.14;0.5;2400.0;5700.0;1.6;10;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +331;500.0;500.0;1.6;2;2400.0;5700.0;0;10;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +332;500.0;500.0;10;10;2400.0;5700.0;0;10;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +333;500.0;500.0;2.13;2.5;2400.0;5700.0;0.6;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +334;500.0;500.0;2.13;2.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +335;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +411;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +412;500.0;500.0;0.13;0.5;2400.0;5700.0;0.1;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +421;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +422;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +423;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +511;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +512;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +521;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +522;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3 +523;500.0;500.0;0.13;0.5;2400.0;5700.0;0;0.6;1.28;8.3;1.0;300.0;70000.0;2300000.0;1.5E7;1800.0;1000.0;600.0;0.3;2.5E-5;4.0;0.3""" + + +def testAnnFuelTable(): + return """Index;Rhod;sd;Ta +0;0.24;0.68;0.097 +1;0.95;0.73;0.44 +2;0.82;0.36;0.050 +3;0.77;0.11;0.05 +4;0.50;0.92;0.22""" + +def genAltitudeMap(slope_coef, sub_sim_shape, data_resolution): + """ + Generate a matrix of altitudes given a slope coefficient + """ + slope = np.linspace(0, 1, sub_sim_shape[1]) + slope = np.repeat(slope, sub_sim_shape[0]) + slope = slope.reshape(sub_sim_shape[1], sub_sim_shape[0]).T + return slope * slope_coef * (data_resolution / 5) + +# Functions definitions + +def get_multi_sub_domain_indices_from_location(x, y, originX, originY, domain_width, domain_height, shape_multisim): + """ + Used for retrieve indices of coordinates inside simulation matrix + """ + i = np.floor(((x - originX) / domain_width) * shape_multisim[0]) + j = np.floor(((y - originY) / domain_height) * shape_multisim[1]) + return int(i), int(j) + +def get_sub_domain_indices_from_location(x, y, originX, originY, domain_width, domain_height): + """ + Used for retrieve indices of coordinates inside simulation matrix + """ + i = np.floor(((x - originX) / domain_width)) + j = np.floor(((y - originY) / domain_height)) + return int(i), int(j) + +def maxDiff(a): + """ + Used for get the maximum difference along first axis of an array + """ + vmin = a[0] + dmax = 0 + for i in range(len(a)): + if (a[i] < vmin): + vmin = a[i] + elif (a[i] - vmin > dmax): + dmax = a[i] - vmin + return dmax + +def getLocationFromLine(line): + """ + Return the location of current node (line). + """ + llv = line.split("loc=(") + if len(llv) < 2: + return None + llr = llv[1].split(",") + if len(llr) < 3: + return None + return (float(llr[0]),float(llr[1])) + +def printToPathe(linePrinted): + """ + Compute the current results of simulation to pathes. + """ + fronts = linePrinted.split("FireFront") + pathes = [] + for front in fronts[1:]: + nodes = front.split("FireNode")[1:] + if len(nodes) > 0: + Path = mpath.Path + codes = [] + verts = [] + firstNode = getLocationFromLine(nodes[0]) + codes.append(Path.MOVETO) + verts.append(firstNode) + + for node in nodes[:]: + newNode = getLocationFromLine(node) + codes.append(Path.LINETO) + verts.append(newNode) + codes.append(Path.LINETO) + verts.append(firstNode) + pathes.append(mpath.Path(verts, codes)) + + return pathes + + + +import struct +import zlib + +def write_png_header(width, height): + # PNG file signature + png_signature = b'\x89PNG\r\n\x1a\n' + + # IHDR chunk: width, height, bit depth, color type, compression, filter, interlace + # Color type 4: grayscale with alpha, Bit depth 8 + ihdr_data = struct.pack(">IIBBBBB", width, height, 8, 4, 0, 0, 0) + return png_signature + create_chunk(b'IHDR', ihdr_data) + +def create_chunk(chunk_type, data): + # Chunk structure: length, type, data, CRC + chunk_length = struct.pack(">I", len(data)) + chunk_type = chunk_type + chunk_crc = struct.pack(">I", zlib.crc32(chunk_type + data) & 0xffffffff) + return chunk_length + chunk_type + data + chunk_crc + +def write_png_data(data, width, height): + # Modify to handle transparency for grayscale value 255 + raw_data = b'' + for y in range(height): + row_data = b'\x00' # Filter type 0 (None) + for x in range(width): + gray = data[y * width + x] + alpha = 0 if gray == 255 else 255 # Transparent if gray is 255 + row_data += struct.pack("BB", gray, alpha) # Grayscale value followed by alpha value + raw_data += row_data + compressed_data = zlib.compress(raw_data) # Compress the data as required by the PNG specification + return create_chunk(b'IDAT', compressed_data) + +def write_png_file(filename, data, width, height): + with open(filename, 'wb') as f: + # Write header + f.write(write_png_header(width, height)) + + # Write image data + f.write(write_png_data(data, width, height)) + + # Write IEND chunk + f.write(create_chunk(b'IEND', b'')) + + +from matplotlib.colors import ListedColormap + +def map_fuel_to_colors(fuelmap, fuel_list): + """ + Convert a fuel_map for use with colors. + The returned fuel_map values are replaced by indices of fuels from fuel_list. + """ + for i in range(len(fuelmap)): + for j in range(len(fuelmap[0])): + try: + fuelmap[i][j] = fuel_list.index(fuelmap[i][j]) + 1 + except ValueError: + fuelmap[i][j] = 0 + return fuelmap + +def fill_random(s, k, value_yes, value_no=0): + """Generate a randomly filled array.""" + a = np.random.random(size=s) + return np.where(a > k, value_yes, value_no) + + +def computeSpeed(atime): + """ + Computes the speed as the inverse of the gradient of arrival times. + 'inf' in the arrival times indicates that the point was never reached. + + Parameters: + atime (np.array): 2D array of arrival times. + + Returns: + np.array: 2D array of speeds, with the same shape as atime. + Returns 'inf' where the arrival time is 'inf', indicating no arrival. + """ + # Check where the times are infinite + inf_mask = np.isinf(atime) + + # Replace 'inf' with the maximum finite value in atime + max_finite = np.nanmax(np.where(inf_mask, np.nan, atime)) + atime_temp = np.where(inf_mask, max_finite, atime) + + # Compute the gradient in the x and y directions on the modified array + grad_y, grad_x = np.gradient(atime_temp) + + # Compute the magnitude of the gradient + grad_mag = np.sqrt(grad_x**2 + grad_y**2) + + # Speed is the inverse of the gradient magnitude + # Avoid division by zero by adding a small number in the denominator + speed = np.where(grad_mag == 0, np.inf, 1 / (grad_mag + 1e-10)) + + # Re-assign 'inf' to the speed array where the original arrival time was 'inf' + speed[inf_mask] = np.inf + + return speed + + +def plot_simulation(pathes, fuel_map, elevation_map, myExtents, scalMap = None): + """ + Used for plot 4 axis graph, with Heatflux, Fuels, Altitude plotted under simulation, + and Statistics for the last axis. + """ + # Create a figure with 2 axis (2 subplots) + fig, ax = plt.subplots(figsize=(10,7), dpi=120) + + # Get fuel_map matrix + if fuel_map is not None: + fuels = fuel_map + + fuel_types = { # ESA EU and COrine + 10: {'color': (0, 100, 0), 'description': 'Tree cover'}, + 20: {'color': (255, 187, 34), 'description': 'Shrubland'}, + 30: {'color': (255, 255, 76), 'description': 'Grassland'}, + 40: {'color': (240, 150, 255), 'description': 'Cropland'}, + 50: {'color': (250, 0, 0), 'description': 'Built-up'}, + 60: {'color': (180, 180, 180), 'description': 'Bare / sparse vegetation'}, + 70: {'color': (240, 240, 240), 'description': 'Snow and Ice'}, + 80: {'color': (0, 100, 200), 'description': 'Permanent water bodies'}, + 90: {'color': (0, 150, 160), 'description': 'Herbaceous wetland'}, + 95: {'color': (0, 207, 117), 'description': 'Mangroves'}, + 100: {'color': (250, 230, 160), 'description': 'Moss and lichen'}, + 0: {'color': (255, 255, 255), 'description': 'Clouds'}, + 62: {'color': (210, 0, 0),'description': 'Artificial surfaces and constructions'}, + 73: {'color': (253, 211, 39), 'description': 'Cultivated areas'}, + 75: {'color': (176, 91, 16), 'description': 'Vineyards'}, + 82: {'color': (35, 152, 0), 'description': 'Broadleaf tree cover'}, + 83: {'color': (8, 98, 0), 'description': 'Coniferous tree cover'}, + 102: {'color': (249, 150, 39), 'description': 'Herbaceous vegetation'}, + 103: {'color': (141, 139, 0), 'description': 'Moors and Heathland'}, + 104: {'color': (95, 53, 6), 'description': 'Sclerophyllous vegetation'}, + 105: {'color': (149, 107, 196), 'description': 'Marshes'}, + 106: {'color': (77, 37, 106), 'description': 'Peatbogs'}, + 121: {'color': (154, 154, 154), 'description': 'Natural material surfaces'}, + 123: {'color': (106, 255, 255),'description': 'Permanent snow covered surfaces'}, + 162: {'color': (20, 69, 249), 'description': 'Water bodies'}, + 255: {'color': (255, 255, 255), 'description': 'No data'}, + 111: {'color': (230, 0, 77), 'description': 'Urban fabric, Continuous urban fabric'}, + 112: {'color': (255, 0, 0), 'description': 'Urban fabric, Discontinuous urban fabric'}, + 121: {'color': (204, 77, 242), 'description': 'Industrial, commercial and transport units, Industrial or commercial units'}, + 122: {'color': (204, 0, 0), 'description': 'Industrial, commercial and transport units, Road and rail networks and associated land'}, + 123: {'color': (230, 204, 204), 'description': 'Industrial, commercial and transport units, Port areas'}, + 124: {'color': (230, 204, 230), 'description': 'Industrial, commercial and transport units, Airports'}, + 131: {'color': (166, 0, 204), 'description': 'Mine, dump and construction sites, Mineral extraction sites'}, + 132: {'color': (166, 77, 0), 'description': 'Mine, dump and construction sites, Dump sites'}, + 133: {'color': (255, 77, 255), 'description': 'Mine, dump and construction sites, Construction sites'}, + 141: {'color': (255, 166, 255), 'description': 'Artificial, non-agricultural vegetated areas, Green urban areas'}, + 142: {'color': (255, 230, 255), 'description': 'Artificial, non-agricultural vegetated areas, Sport and leisure facilities'}, + 211: {'color': (255, 255, 168), 'description': 'Agricultural areas, Arable land, Non-irrigated arable land'}, + 212: {'color': (255, 255, 0), 'description': 'Agricultural areas, Arable land, Permanently irrigated land'}, + 213: {'color': (230, 230, 0), 'description': 'Agricultural areas, Arable land, Rice fields'}, + 221: {'color': (230, 128, 0), 'description': 'Agricultural areas, Permanent crops, Vineyards'}, + 222: {'color': (242, 166, 77), 'description': 'Agricultural areas, Permanent crops, Fruit trees and berry plantations'}, + 223: {'color': (230, 166, 0), 'description': 'Agricultural areas, Permanent crops, Olive groves'}, + 231: {'color': (230, 230, 77), 'description': 'Agricultural areas, Pastures'}, + 241: {'color': (255, 230, 166), 'description': 'Agricultural areas, Heterogeneous agricultural areas, Annual crops associated with permanent crops'}, + 242: {'color': (255, 230, 77), 'description': 'Agricultural areas, Heterogeneous agricultural areas, Complex cultivation patterns'}, + 243: {'color': (230, 204, 77), 'description': 'Agricultural areas, Heterogeneous agricultural areas, Land principally occupied by agriculture, with significant areas of natural vegetation'}, + 244: {'color': (242, 204, 166), 'description': 'Agricultural areas, Heterogeneous agricultural areas, Agro-forestry areas'}, + 311: {'color': (128, 255, 0), 'description': 'Forest and semi natural areas, Forests, Broad-leaved forest'}, + 312: {'color': (0, 166, 0), 'description': 'Forest and semi natural areas, Forests, Coniferous forest'}, + 313: {'color': (77, 255, 0), 'description': 'Forest and semi natural areas, Forests, Mixed forest'}, + 321: {'color': (204, 242, 77), 'description': 'Forest and semi natural areas, Scrub and/or herbaceous vegetation associations, Natural grasslands'}, + 322: {'color': (166, 255, 128), 'description': 'Forest and semi natural areas, Scrub and/or herbaceous vegetation associations, Moors and heathland'}, + 323: {'color': (166, 230, 77), 'description': 'Forest and semi natural areas, Scrub and/or herbaceous vegetation associations, Sclerophyllous vegetation'}, + 324: {'color': (166, 242, 0), 'description': 'Forest and semi natural areas, Scrub and/or herbaceous vegetation associations, Transitional woodland-shrub'}, + 331: {'color': (230, 230, 230), 'description': 'Forest and semi natural areas, Open spaces with little or no vegetation, Beaches, dunes, sands'}, + 332: {'color': (204, 204, 204), 'description': 'Forest and semi natural areas, Open spaces with little or no vegetation, Bare rocks'}, + 333: {'color': (204, 255, 204), 'description': 'Forest and semi natural areas, Open spaces with little or no vegetation, Sparsely vegetated areas'}, + 334: {'color': (0, 0, 0), 'description': 'Forest and semi natural areas, Open spaces with little or no vegetation, Burnt areas'}, + 335: {'color': (166, 230, 204), 'description': 'Forest and semi natural areas, Open spaces with little or no vegetation, Glaciers and perpetual snow'}, + 411: {'color': (166, 166, 255), 'description': 'Wetlands, Inland wetlands, Inland marshes'}, + 412: {'color': (77, 77, 255), 'description': 'Wetlands, Inland wetlands, Peat bogs'}, + 421: {'color': (204, 204, 255), 'description': 'Wetlands, Maritime wetlands, Salt marshes'}, + 422: {'color': (230, 230, 255), 'description': 'Wetlands, Maritime wetlands, Salines'}, + 423: {'color': (166, 166, 230), 'description': 'Wetlands, Maritime wetlands, Intertidal flats'}, + 511: {'color': (0, 204, 242), 'description': 'Water bodies, Inland waters, Water courses'}, + 512: {'color': (128, 242, 230), 'description': 'Water bodies, Inland waters, Water bodies'}, + 521: {'color': (0, 255, 166), 'description': 'Water bodies, Marine waters, Coastal lagoons'}, + 522: {'color': (166, 255, 230), 'description': 'Water bodies, Marine waters, Estuaries'}, + 523: {'color': (230, 242, 255), 'description': 'Water bodies, Marine waters, Sea and ocean'} + } + # Normalize the colors to the [0, 1] range expected by matplotlib + for key in fuel_types: + rgb = fuel_types[key]['color'] + fuel_types[key]['color'] = tuple([x/255.0 for x in rgb]) + + # Generate the list of colors + colors = [fuel_types[key]['color'] for key in sorted(fuel_types.keys())] + + # Create the colormap + lcmap = ListedColormap(colors) + + CS = ax.imshow(fuels, cmap=lcmap, interpolation='nearest', origin='lower', extent=myExtents) + # norm = mcolors.BoundaryNorm(bounds, cmap.N) + plt.colorbar(CS) + + + if elevation_map is not None: + elevation = elevation_map#np.transpose(elevation_map.reshape(elevation_map.shape[0], elevation_map.shape[1], 1))[0] + + contour_levels = np.arange(np.min(elevation), np.max(elevation), 200) # Contours every 200m + ax.contour(elevation, levels=contour_levels, colors='white', origin='lower', extent=myExtents, linewidths=0.5, linestyles='solid') + + if scalMap is not None: + CS = ax.imshow(scalMap, origin='lower', extent=myExtents) + plt.colorbar(CS) + + # bounds = [0, 1] + # cmap = mpl.colors.ListedColormap(['brown', 'green']) + # norm = mpl.colors.BoundaryNorm(bounds, cmap.N) + # CS = ax.imshow(fuels, cmap=cmap, interpolation='nearest', origin='lower', extent=(0,domain_width,0,domain_height)) + # # plt.colorbar(plt.cm.ScalarMappable(cmap=cmap, norm=norm), cax=ax) + # plt.colorbar(CS, norm=norm) + + path_colors = ['red', 'orange', 'yellow', 'white'] + + # Plot current firefronts to the first 3 subplots + for p, path in enumerate(pathes): + patch = mpatches.PathPatch(path, edgecolor=path_colors[p % len(path_colors)], facecolor='none', alpha=1, lw=2) + ax.add_patch(patch) + ax.grid() + ax.axis('equal') + + ax.grid() + ax.axis('equal') + plt.show() + + + +def print_ff_script(file_path): + """ + Executes each line of a .ff script file using the ff.execute() method, + after stripping leading and trailing whitespace, including tabs. + """ + # Open the file with the script + with open(file_path, 'r') as file: + # Iterate over each line in the file + for line in file: + # Strip leading and trailing whitespace, including tabs, from the line + clean_line = line.strip() + # Check if the cleaned line is not empty + if clean_line: + # Execute the cleaned line with ff.execute, ensuring no leading/trailing whitespace + command = f'ff.execute("{clean_line}")' + print(command) # This is for demonstration; replace with actual ff.execute(clean_line) in use + # If you have the ff object with an execute method available, you would call it here: + # ff.execute(clean_line) diff --git a/tools/preprocessing/kmlDomain.py b/tools/preprocessing/kmlDomain.py index 26be3f1b..174aa140 100644 --- a/tools/preprocessing/kmlDomain.py +++ b/tools/preprocessing/kmlDomain.py @@ -1,4 +1,4 @@ - #!python +#!python import sys import netCDF4 import numpy as np @@ -7,10 +7,10 @@ #fnames = ["/Users/filippi_j/Volumes/orsu/firecaster/2023/nest150Ref/001_pgd/PGD_D2000mA.nested.nc","/Users/filippi_j/Volumes/orsu/firecaster/2023/nest150Ref/001_pgd/PGD_D400mA.nested.nc","/Users/filippi_j/Volumes/orsu/firecaster/2023/nest150Ref/001_pgd/PGD_D80mA.nested.nc"] #foutName = "/Users/filippi_j/data/2023/corbara20230727/domains.kml" -#if len(sys.argv) != 2: -# print("usage kmlDomain PGDFile1 2....") -#else : -# fnames = sys.argv[1:] +if len(sys.argv) != 2: + print("usage kmlDomain PGDFile1 2....") +else : + fnames = sys.argv[1:] def pgds_to_KML_nc4(fnames, fout): dbo={} diff --git a/tools/preprocessing/parametersProperties.py b/tools/preprocessing/parametersProperties.py new file mode 100755 index 00000000..ef96c246 --- /dev/null +++ b/tools/preprocessing/parametersProperties.py @@ -0,0 +1,57 @@ +heatFlux = {'type':'flux', \ + 'indices': [0], \ + 'model0name': 'na', \ + } +vaporFlux= {'type':'flux', \ + 'indices': [1], \ + 'model1name': 'na',\ + } +scalarFlux = {'type':'flux', \ + 'indices': [2], \ + 'model2name': 'na', \ + } + +fuel = {'type' : "fuel"} +windU = {'type' : "data"} +windV = {'type' : "data"} +altitude = {'type' : "data"} +#fire regime control time +FromObs_evaporationTime = {'type' : "data"} +FromObs_residenceTime = {'type' : "data"} +FromObs_burningTime = {'type' : "data"} +#sensible heat & vapor flux +FromObs_VaporFlux_evaporation = {'type' : "data"} +FromObs_NominalHeatFlux_flaming = {'type' : "data"} +FromObs_NominalHeatFlux_smoldering = {'type' : "data"} +#Emission Factor +FromObs_EFScalar_flamning = {'type' : "data"} +FromObs_EFScalar_smoldering = {'type' : "data"} +#radiation fraction +FromObs_radiationFraction_flaming = {'type' : "data"} +FromObs_radiationFraction_smoldering = {'type' : "data"} +#conversion factor +FromObs_conversionFactor = {'type' : "data"} + +domain = {'type': 'domain', + 'NEx' : None , + 'NEy' : None, + 'NEz' : None, + 'SWx' : None, + 'SWy' : None, + 'SWz' : None, + 'Lx' : None, + 'Ly' : None, + 'Lz' : None, + 't0' : None, + 'Lt' : None + } + +parameters= {'type' : "parameters" , + #'date' : None , #"2000-01-01_00:00:00", + 'refYear' : None, #2000 , + #'refMonth' : None, #1, + 'refDay' : None, #1, + 'duration' : None, #360000 , + 'projection' : None, #"OPENMAP" , + 'projectionproperties' : None, #"41.551998,8.828396,41.551998,8.828396" , + } diff --git a/tools/preprocessing/run_data.py b/tools/preprocessing/run_data.py new file mode 100644 index 00000000..a8197db6 --- /dev/null +++ b/tools/preprocessing/run_data.py @@ -0,0 +1,218 @@ +import numpy as np +import matplotlib as mpl +import matplotlib.pyplot as plt +import xarray as xr +import rioxarray +import pyforefire as forefire +import os +import sys +import argparse +import importlib +#import ffToGeoJson +import pyproj +import glob +import rasterio +from skimage import exposure +import cv2 +import pdb + +#homebrewed +import create_data_tools as tools +from ffToGeo import FFGEO +from satpy import Scene +from satpy.writers import get_enhanced_image + +from forefire_helper import * + +################################################### +def getDomainExtent(line): + print(line) + llv = line.split("sw=(") + llr = llv[1].split("ne=("); + return( float( llr[0].split(",")[0]), float(llr[1].split(",")[0]), float(llr[0].split(",")[1]), float(llr[1].split(",")[1]) + ) +################################################### +if __name__ == '__main__': +################################################### + + importlib.reload(tools) + print('check env variable:') + print('ntask =', os.environ['ntask']) + + + #create dir + tools.ensure_dir('./ForeFire/') + tools.ensure_dir('./ForeFire/Outputs/') + + mydatanc = './ForeFire/mydata_ll.nc' + s2filedir = '/data/paugam/Data/ElPontDeVilamora/Sentinel/S2/S2A_MSIL1C_20220716T103641_N0400_R008_T31TDG_20220716T155924.SAFE/' + + #load fuelmap + #----------- + mydata = xr.open_dataset(mydatanc) + xxmnh = np.linspace(mydata.domain.attrs['SWx'], mydata.domain.attrs['NEx'], mydata.sizes['DIMX']) + yymnh = np.linspace(mydata.domain.attrs['SWy'], mydata.domain.attrs['NEy'], mydata.sizes['DIMY']) + + #fuelmap_f = xr.DataArray( + # mydata.fuel[0,0], + # coords={"x": xx, "y": yy}, + # dims=["y", "x"], # Order of dimensions matches the data array shape + # ) + fuelmap_ll = mydata.fuel[0,0,:,:] + #fuelmap_f.name = 'fuel map' + #fuelmap_f = fuelmap_f.rio.write_crs('epsg:{:d}'.format(mydata.domain.attrs['epsg'])) + #fuelmap_ll = fuelmap_f.rio.reproject(4326) + + ''' + N = len(np.unique(fuelmap_ll.values)) -1 # Number of colors + colors = np.random.rand(N, 3) #plt.cm.hsv(np.linspace(0, 1, N)) + cmap_fm = mcolors.ListedColormap(colors) + colors[0,:] = [1,0,0] + + # Define the boundaries of each color segment + bounds = np.unique(fuelmap_ll.values)[1:] # from 1 to 10 + bounds = np.append(bounds, 9999) + norm_fm = mcolors.BoundaryNorm(bounds, cmap_fm.N) + + #set no data to not be ploted + nodata_value = -9999 + fuelmap_ll = fuelmap_ll.where(fuelmap_ll != nodata_value) + nx,ny = fuelmap_ll.shape + #plot fuelmap + mpl.rcParams['figure.subplot.left'] = .0 + mpl.rcParams['figure.subplot.right'] = 1. + mpl.rcParams['figure.subplot.top'] = 1. + mpl.rcParams['figure.subplot.bottom'] = .0 + fig = plt.figure(figsize=(10,10*ny/nx)) + fuelmap_ll.plot(cmap=cmap_fm,norm=norm_fm,add_colorbar=False) + + fig.savefig('fuelmap.png',dpi=200) + plt.close(fig) + + #load RGB sentinel2 + ################### + postfire_rgb_path = s2filedir+'preFireRGB-s2.tif' + if ~os.path.isfile(postfire_rgb_path): + # Open b2, b3 and b4 + band2=rasterio.open(glob.glob(s2filedir+'GRANULE/L1C_T31TDG_A036899_20220716T104613/IMG_DATA/'+'*B02.jp2')[0]) + band3=rasterio.open(glob.glob(s2filedir+'GRANULE/L1C_T31TDG_A036899_20220716T104613/IMG_DATA/'+'*B03.jp2')[0]) + band4=rasterio.open(glob.glob(s2filedir+'GRANULE/L1C_T31TDG_A036899_20220716T104613/IMG_DATA/'+'*B04.jp2')[0]) + + # Extract and update the metadata + meta = band2.meta + meta.update({"count": 3}) + + # Write the natural colour composite image with metadata + postfire_rgb_path = s2filedir+'preFireRGB-s2.tif' + + with rasterio.open(postfire_rgb_path, 'w', **meta) as dest: + dest.write(band4.read(1),1) + dest.write(band3.read(1),2) + dest.write(band2.read(1),3) + + # Transpose and rescale the image + s2 = rioxarray.open_rasterio(postfire_rgb_path) + s2 = s2.rio.reproject(fuelmap_ll.rio.crs) + s2 =s2.rename({'x':'DIMX'}) + s2 =s2.rename({'y':'DIMY'}) + s2_pgd3 = s2.interp(DIMX=fuelmap_ll.DIMX, DIMY=fuelmap_ll.DIMY) + + # Assuming 'image' is your 3-band RGB image + def equalize_histogram(image): + # Split the image into its color channels + b, g, r = cv2.split(image) + + # Apply histogram equalization on each channel + b_eq = cv2.equalizeHist((255*b).astype(np.uint8)) + g_eq = cv2.equalizeHist((255*g).astype(np.uint8)) + r_eq = cv2.equalizeHist((255*r).astype(np.uint8)) + + # Merge the equalized channels back + image_eq = cv2.merge((b_eq, g_eq, r_eq)) + + return image_eq + + ir = 0 + ig = 1 + ib = 2 + #image = np.array([s2_pgd3[0],s2_pgd3[1],s2_pgd3[2]]).transpose(1,2,0) + image = np.array([s2_pgd3[ir][::-1,:]/s2_pgd3[ir].max(), + s2_pgd3[ig][::-1,:]/s2_pgd3[ig].max(), + s2_pgd3[ib][::-1,:]/s2_pgd3[ib].max()] + ).transpose(1,2,0) + + #p2, p98 = np.percentile(image, (2, 98)) + #image_contrast_stretched = np.clip((image - p2) / (p98 - p2), 0, 1) + + # Equalize the histogram of the image + image_equalized = equalize_histogram(image) + + # Plot the resulting image + ny,nx,_ = image_equalized.shape + fig = plt.figure(figsize=(10,10*ny/nx)) + plt.imshow(np.transpose(image_equalized,[1,0,2])) + fig.savefig('s2RGB.png',dpi=200) + plt.close(fig) + ''' + ################### ... run forefire + + ff = forefire.ForeFire() + outputsUpdate = 10 + spatialIncrement = 2 + perimeterResolution = 10 + minimalPropagativeFrontDepth = 10 + bmapOutputUpdate = 10 + InitTime = 0 + + #set param + ff.setString('ForeFireDataDirectory','ForeFire') + ff.setString('fireOutputDirectory','ForeFire/Outputs') + ff.setInt('outputsUpdate',outputsUpdate) + + ff.setString('fuelsTableFile','fuels.ff') # end dur + ff["propagationModel"]="Andrews" + + ff.setDouble("spatialIncrement",spatialIncrement) + ff.setDouble("perimeterResolution",perimeterResolution) + ff.setDouble("minimalPropagativeFrontDepth",minimalPropagativeFrontDepth) + + ff.setDouble("bmapOutputUpdate",bmapOutputUpdate) + ff.setInt("defaultHeatType",0) + ff.setInt("defaultVaporType",1) + + ff.setDouble("InitTime",InitTime) + ff.execute("loadData[mydata.nc;2022-07-16T00:00:00Z]") + #ff.execute("startFire[loc=(60390,76090,0.);t=0]") + ff.execute("startFire[loc=(59601.3,77044.1,0.);t=39600]") + print(ff.execute("print[]")) + + ff.setString('dumpMode','FF') + + N_step = 10 + step = 20 + t0 = 39000 + pathes = [] + ffgeo = FFGEO('pdv',4326) + ffgeo.addFront(ff) + + for i in np.arange(1,N_step): + + ff_time = t0 + i*step + + print("goTo[t=%f]"%(ff_time), end='\n') + ff.execute("goTo[t=%f]"%(ff_time)) + + pathes += printToPathe(ff.execute("print[]")) + #print(pathes) + ffgeo.addFront(ff) + + + #fig = plt.figure() + #ax = plt.subplot(111) + #ax.pcolormesh(xx,yy,mydata.fuel[0,0]) + + plt.figure() + ax = plt.subplot() + ax.pcolormesh(xxmnh,yymnh,mydata.fuel[0,0]) + ffgeo.pd.plot(ax=ax,facecolor='None') + plt.show() diff --git a/tools/preprocessing/run_data_simple.py b/tools/preprocessing/run_data_simple.py new file mode 100644 index 00000000..d639f4dc --- /dev/null +++ b/tools/preprocessing/run_data_simple.py @@ -0,0 +1,64 @@ +import numpy as np +import pyforefire as forefire +#from ffToGeo import FFGEO +#from forefire_helper import * + +################################################### +if __name__ == '__main__': +################################################### + + ff = forefire.ForeFire() + outputsUpdate = 10 + spatialIncrement = 2 + perimeterResolution = 10 + minimalPropagativeFrontDepth = 10 + bmapOutputUpdate = 10 + InitTime = 0 + + #set param + ff.setString('ForeFireDataDirectory','ForeFire') + ff.setString('fireOutputDirectory','ForeFire/Outputs') + ff.setInt('outputsUpdate',outputsUpdate) + + ff.setString('fuelsTableFile','fuels.ff') # end dur + ff["propagationModel"]="Andrews" + + ff.setDouble("spatialIncrement",spatialIncrement) + ff.setDouble("perimeterResolution",perimeterResolution) + ff.setDouble("minimalPropagativeFrontDepth",minimalPropagativeFrontDepth) + + ff.setDouble("bmapOutputUpdate",bmapOutputUpdate) + ff.setInt("defaultHeatType",0) + ff.setInt("defaultVaporType",1) + + ff.setDouble("InitTime",InitTime) + ff.execute("loadData[mydata.nc;2022-07-16T01:30:00Z]") + ff.execute("startFire[loc=(59601.3,77044.1,0.);t=39600]") + print(ff.execute("print[]")) + + ff.setString('dumpMode','FF') + + N_step = 10 + step = 20 + t0 = 39000 + pathes = [] + #ffgeo = FFGEO('pdv',4326) + #ffgeo.addFront(ff) + + for i in np.arange(1,N_step): + + ff_time = t0 + i*step + + print("goTo[t=%f]"%(ff_time), end='\n') + ff.execute("goTo[t=%f]"%(ff_time)) + + #pathes += printToPathe(ff.execute("print[]")) + #print(pathes) + #ffgeo.addFront(ff) + + + #fig = plt.figure() + #ax = plt.subplot(111) + #ax.pcolormesh(xx,yy,mydata.fuel[0,0]) + +