diff --git a/src/HeatFluxBasicModel.cpp b/src/HeatFluxBasicModel.cpp index 94428cb..1c1ab9a 100644 --- a/src/HeatFluxBasicModel.cpp +++ b/src/HeatFluxBasicModel.cpp @@ -101,7 +101,7 @@ string HeatFluxBasicModel::getName(){ /* ****************** */ /* Model for the flux */ -/* ****************** */ +/* *** *************** */ double HeatFluxBasicModel::getValue(double* valueOf , const double& bt, const double& et, const double& at){ diff --git a/tools/preprocessing/extract_from_eu_land_cover.py b/tools/preprocessing/extract_from_eu_land_cover.py index d51eb2a..2c4164b 100644 --- a/tools/preprocessing/extract_from_eu_land_cover.py +++ b/tools/preprocessing/extract_from_eu_land_cover.py @@ -9,25 +9,26 @@ import rasterio from rasterio.windows import from_bounds from pyproj import Transformer -import matplotlib.pyplot as plt +import fiona import numpy as np -from pyproj import CRS +import os import osmnx as ox from rasterio.features import geometry_mask -from rasterio.warp import reproject, Resampling +from rasterio.warp import calculate_default_transform, reproject, Resampling + from rasterio.transform import from_origin import geopandas as gpd from rasterio.mask import mask from shapely.geometry import box import json -from PIL import Image -import rasterio -from rasterio.warp import calculate_default_transform, reproject, transform_bounds, Resampling -import numpy as np + from .ffToGeoJson import create_kml,generate_indexed_png_and_legend - + +from rasterio.features import rasterize +from shapely.geometry import mapping + attribute_widths_road_edge_full = { 'secondary': 0.8, 'track': 0.8, @@ -100,9 +101,16 @@ 'primary': 1.5, } +# Example usage +attribute_widths_waterway = { + 'river': 1.5, + 'stream': 0.5, + 'canal': 1.0, + 'drain': 0.7, + 'generic': 0.5 # Set a generic width for all water features +} + def extract_subregion(input_tif, output_tif, westI, southI, eastI, northI): - print(f"clipping {input_tif} into {output_tif} bounds W{westI}, S{southI}, E{eastI}, N{northI}") - print(f"NW({northI},{westI}), SW({southI},{westI}), SE({southI},{eastI}), NE({northI},{eastI})") with rasterio.open(input_tif) as src: # Transformer pour convertir les coordonnées lat/lon (WGS84) en coordonnées du système du raster @@ -214,15 +222,234 @@ def warp_and_clip_raster(input_file_path: str, warped_file_path: str, output_fil -def extract_roads_from_geotiff(west,south,east,north, road_shape_file): - print(f"Getting roads with bounds from {west,south,east,north} saving {road_shape_file}") +def extract_roads_from_geotiff_original(west,south,east,north, road_shape_file): try: G = ox.graph_from_bbox(north, south, east, west, network_type='all') except Exception as e: return str(e) ox.save_graph_shapefile(G, filepath=road_shape_file ) -def rasterize_shapefile(shapefile_path, ref_tif, output_path, attribute_widths=attribute_widths_road_Edge, default_width = 0.3, imbounds =None,code=62): + +def extract_roads_from_geotiff(west, south, east, north, road_shape_file, attribute_widths=attribute_widths_road_Edge): + + + # Define the road types to extract based on the attribute_widths_road_Edge dictionary keys + road_types = list(attribute_widths.keys()) + + # Create a custom filter to fetch only the roads of interest + custom_filter = f'["highway"~"{ "|".join(road_types) }"]' + + try: + # Define the bounding box + bbox = (north, south, east, west) + + # Retrieve the graph with the custom filter using the new bbox parameter + G = ox.graph_from_bbox(bbox=bbox, network_type='drive', custom_filter=custom_filter) + except Exception as e: + return str(e) + + if len(G.edges) <= 0: + return None + # You might want to further process the graph here to assign or use your custom attributes + # For example, adjusting edge widths as per your dictionary before saving + for u, v, key, data in G.edges(data=True, keys=True): + road_type = data.get('highway', None) + if isinstance(road_type, list): # Sometimes 'highway' can be a list of types + road_type = road_type[0] # Select the first type as the representative + width = attribute_widths.get(road_type, None) + if width is not None: + data['width'] = width + + # Save the graph to a shapefile + ox.save_graph_geopackage(G, filepath=road_shape_file) + + return road_shape_file + + + +def extract_waterways_from_geotiff(west, south, east, north, waterway_shape_file, attribute_widths): + + # Define the waterway types to extract based on the attribute_widths dictionary keys + waterway_types = list(attribute_widths.keys()) + + # Create a custom filter to fetch only the waterways of interest + custom_filter = f'["waterway"~"{ "|".join(waterway_types) }"]' + + + try: + bbox = (north, south, east, west) + # Retrieve the graph with the custom filter using the new 'bbox' parameter + G = ox.graph_from_bbox(bbox=bbox, network_type='none', custom_filter=custom_filter, retain_all=True) + except Exception as e: + return str(e) + + if len(G.edges) <= 0: + return None + + # You might want to further process the graph here to assign or use your custom attributes + # For example, adjusting widths as per your dictionary before saving + for u, v, key, data in G.edges(data=True, keys=True): + waterway_type = data.get('waterway', None) + if isinstance(waterway_type, list): # Sometimes 'waterway' can be a list of types + waterway_type = waterway_type[0] # Select the first type as the representative + width = attribute_widths.get(waterway_type, None) + if width is not None: + data['width'] = width + + # Save the graph to a shapefile + ox.save_graph_geopackage(G, filepath=waterway_shape_file) + return waterway_shape_file + +def rasterize_shapefiles_check(road_shapefile_path, water_shapefile_path, ref_tif, output_path, + roads_attribute_widths, water_attribute_widths, + default_width=0.3, imbounds=None, road_code=62, water_code=162): + + shapes = [] + + if os.path.exists(road_shapefile_path): + roads_gdf = gpd.read_file(road_shapefile_path) + for _, feature in roads_gdf.iterrows(): + if('highway' in feature.keys()): + highway_type = feature['highway'] + line_width = roads_attribute_widths.get(highway_type, default_width) + if line_width > 0: + buffered_geom = feature['geometry'].buffer(line_width / 2) # Adjust this buffer as needed + shapes.append((mapping(buffered_geom), road_code)) + else: + print(f"Warning: Road shapefile {road_shapefile_path} not found.") + + if os.path.exists(water_shapefile_path): + water_gdf = gpd.read_file(water_shapefile_path) + generic_water_width = water_attribute_widths.get('generic', default_width) # Assuming a generic type + for _, feature in water_gdf.iterrows(): + buffered_geom = feature['geometry'].buffer(generic_water_width / 2) # Adjust this buffer as needed + shapes.append((mapping(buffered_geom), water_code)) + else: + print(f"Warning: Water shapefile {water_shapefile_path} not found.") + + with rasterio.open(ref_tif) as refT: + ref_crs = refT.crs + width = refT.width + height = refT.height + transform = refT.transform # Directly use the transform from the reference raster + + # Create a base raster array + raster = refT.read(1) + refT.close() + + if shapes: + # Rasterize the shapes onto the raster array only if there are shapes to rasterize + rasterized = rasterize(shapes, out_shape=(height, width), fill=0, transform=transform, dtype=np.uint8) + raster[rasterized == road_code] = road_code + raster[rasterized == water_code] = water_code + + # Save the raster data to a new GeoTIFF file + with rasterio.open(output_path, 'w', driver='GTiff', width=width, height=height, count=1, + dtype=raster.dtype, crs=ref_crs, transform=transform) as dst: + dst.write(raster, 1) + + refT.close() + +def rasterize_shapefiles(road_shapefile_path, water_shapefile_path, ref_tif, output_path, + roads_attribute_widths, water_attribute_widths, + default_width=0.3, imbounds=None, road_code=62, water_code=162): + + roads_gdf = None + water_gdf = None + + if road_shapefile_path is not None: + if os.path.exists(road_shapefile_path): + layers = fiona.listlayers(road_shapefile_path) + if "edges" in layers: + roads_gdf = gpd.read_file(road_shapefile_path,layer="edges") + print(f"Loaded road 'edges' layer with {len(roads_gdf)} shapes.") + + + if water_shapefile_path is not None: + if os.path.exists(water_shapefile_path): + layers = fiona.listlayers(water_shapefile_path) + if "edges" in layers: + water_gdf = gpd.read_file(water_shapefile_path,layer="edges") + print(f"Loaded water stream 'edges' layer with {len(water_gdf)} shapes.") + + # Use specified bounds if provided, otherwise derive bounds from both GeoDataFrames + if imbounds is not None: + bounds = imbounds + else: + print("WARNING getting bounds from roads") + + bounds = [ + roads_gdf.total_bounds[0], # minx + roads_gdf.total_bounds[1], # miny + roads_gdf.total_bounds[2], # maxx + roads_gdf.total_bounds[3] # maxy + ]# + + x_min, y_min, x_max, y_max = bounds + + # Load the reference raster + with rasterio.open(ref_tif) as refT: + ref_crs = refT.crs + width = refT.width + height = refT.height + + # Calculate the resolution based on reference raster + resolutionx = (x_max - x_min) / width + resolutiony = (y_max - y_min) / height + + # Create a transformation from the bounds + transform = rasterio.transform.from_origin(x_min, y_max, resolutionx, resolutiony) + + # Read the existing raster data + raster = refT.read(1) + #raster *= 0 + refT.close() + + # Prepare shapes and associated codes for roads and waterways + shapes = [] + + if roads_gdf is not None: + print("using road shapefile :", road_shapefile_path) + for _, feature in roads_gdf.iterrows(): + if('highway' in feature.keys()): + highway_type = feature['highway'] + # print("YESSSSSSSSS HIGHWAY ") + line_width = roads_attribute_widths.get(highway_type, default_width) + if line_width > 0: + buffered_geom = feature['geometry'].buffer(line_width * resolutionx / 2) + shapes.append((mapping(buffered_geom), road_code)) + else: + print("NO HIGHWAY ", feature.keys()) + + if water_gdf is not None: + print("using water shapefile :", water_shapefile_path) + generic_water_width = water_attribute_widths.get('generic', default_width) # Assuming a generic type if no specific 'waterway' key + for _, feature in water_gdf.iterrows(): + buffered_geom = feature['geometry'].buffer(generic_water_width * resolutionx / 2) + shapes.append((mapping(buffered_geom), water_code)) + + + # Rasterize the shapes onto the raster array + if len(shapes) > 0: + print("Number of shapes to resterize :", len(shapes)) + rasterized = rasterize(shapes, out_shape=(height, width), fill=0, transform=transform, dtype=np.uint8) + raster[rasterized == road_code] = road_code + raster[rasterized == water_code] = water_code + + + # for _, feature in roads_gdf.iterrows(): + # highway_type = feature['highway'] + # line_width = 10#attribute_widths.get(highway_type, default_width) # Default line width is 1 pixel + # mask = geometry_mask([feature['geometry'].buffer(line_width * resolutionx / 2)], + # transform=transform, invert=True, out_shape=(height, width)) + # raster[mask] = road_code # Set pixel value to 255 (white) where there is a feature + # Save the raster data to a new GeoTIFF file + with rasterio.open(output_path, 'w', driver='GTiff', width=width, height=height, count=1, + dtype=raster.dtype, crs=ref_crs, transform=transform) as dst: + dst.write(raster, 1) + + +def rasterize_shapefile_original(shapefile_path, ref_tif, output_path, attribute_widths=attribute_widths_road_Edge, default_width = 0.3, imbounds =None,code=62): print(f"Rasterizing roads {shapefile_path} to {output_path} with {attribute_widths}") reference_properties = {} @@ -363,14 +590,14 @@ def rasterize_kml(kml_path, ref_tif, output_path, default_value=1, imbounds=None crs=ref_crs, transform=transform) as dst: dst.write(raster, 1) -def landcover_roads_to_fuel(S2GLC_tif,legend_file_path, WSEN, LBRT,output_dir,fuel_modifier=None, fuel_resolution = 10): +def landcover_roads_to_fuel(S2GLC_tif,legend_file_path, WSEN, LBRT,output_dir,fuel_modifier=None, fuel_resolution = 10, no_fuel_code = 62): fuel_indices_origin = f"{output_dir}/fuel_indices_S2GLC.tif" fuel_indices_warped = f"{output_dir}/fuel_indices_warped.tif" fuel_indices_warped_clipped = f"{output_dir}/fuel_indices_warped_clipped.tif" roads_shape = f"{output_dir}/roads_shape.shp" - + water_shape = f"{output_dir}/water_shape.shp" fuel_road_indices = f'{output_dir}/fuel.tif' fuel_mod_road_indices = f'{output_dir}/fuelMod.tif' @@ -385,15 +612,20 @@ def landcover_roads_to_fuel(S2GLC_tif,legend_file_path, WSEN, LBRT,output_dir,fu west, south, east, north =WSEN l,b,r,t = LBRT - - - - extract_subregion(S2GLC_tif, fuel_indices_origin, west, south, east, north) + warp_and_clip_raster(fuel_indices_origin,fuel_indices_warped,fuel_indices_warped_clipped ,west, south, east, north,target_width = int((r-l)/fuel_resolution),target_height = int((t-b)/fuel_resolution)) + + roads_shape = extract_roads_from_geotiff(west, south, east, north, roads_shape,attribute_widths=attribute_widths_road_Edge) - extract_roads_from_geotiff(west, south, east, north, roads_shape) - rasterize_shapefile(roads_shape,fuel_indices_warped_clipped,fuel_road_indices,imbounds=(west, south,east ,north),code=62) + + water_shape = extract_waterways_from_geotiff(west, south, east, north, water_shape, attribute_widths=attribute_widths_waterway) + + + rasterize_shapefiles(roads_shape,water_shape,fuel_indices_warped_clipped,fuel_road_indices,roads_attribute_widths=attribute_widths_road_Edge,water_attribute_widths=attribute_widths_waterway,imbounds=(west, south,east ,north), road_code=no_fuel_code, water_code=162) + + + print("Rasterized roads and water") reftif = fuel_road_indices @@ -402,6 +634,8 @@ def landcover_roads_to_fuel(S2GLC_tif,legend_file_path, WSEN, LBRT,output_dir,fu fiona.drvsupport.supported_drivers['KML'] = 'rw' rasterize_kml(fuel_modifier,fuel_road_indices,fuel_mod_road_indices) reftif = fuel_mod_road_indices + print("applied ",fuel_modifier) generate_indexed_png_and_legend(legend_file_path,reftif, fuel_png, fuel_cbar_png) create_kml(west, south, east, north, "FUEL", fuel_png,fuel_kml , pngcbarfile=fuel_cbar_png) + return fuel_png, fuel_kml, fuel_cbar_png diff --git a/tools/preprocessing/ffToGeoJson.py b/tools/preprocessing/ffToGeoJson.py index 1a39559..7fc1aa6 100644 --- a/tools/preprocessing/ffToGeoJson.py +++ b/tools/preprocessing/ffToGeoJson.py @@ -773,6 +773,8 @@ def normalize_rgb(color): # Find unique values in the TIFF data unique_values = np.unique(tif_data) + if len(unique_values) < 2: + print("WARNING ::: Found ",len(unique_values), " in fuel data") # Filter color_palette and class_names to include only unique values filtered_color_palette = {key: color_palette[key] for key in unique_values if key in color_palette} @@ -799,7 +801,8 @@ def normalize_rgb(color): ax.set_xlim(0, 4) # Adjust limits to prevent clipping ax.set_ylim(0, len(all_labels)) ax.set_aspect('auto') - plt.savefig(output_legend_png_path, transparent=True) + plt.tight_layout() + plt.savefig(output_legend_png_path, transparent=True,bbox_inches='tight') diff --git a/tools/preprocessing/kmlDomain.py b/tools/preprocessing/kmlDomain.py index 7d02d0d..26be3f1 100644 --- a/tools/preprocessing/kmlDomain.py +++ b/tools/preprocessing/kmlDomain.py @@ -12,7 +12,7 @@ #else : # fnames = sys.argv[1:] -def pgds_to_KML(fnames, fout): +def pgds_to_KML_nc4(fnames, fout): dbo={} kmlString = str(""" @@ -66,3 +66,37 @@ def pgds_to_KML(fnames, fout): with open(fout, "w") as text_file: text_file.write(kmlString) + + +def bounds_to_KML(boundaries, fout): + kmlString = """ + + + kmldomain.kml + """ + + w, s, e, n = boundaries + + kmlString += f""" + + MNH Bound in PGD + #doms + + 0 + + + + {w},{s},0 {e},{s},0 {e},{n},0 {w},{n},0 {w},{s},0 + + + + + + """ + + kmlString += """ + + """ + + with open(fout, "w") as text_file: + text_file.write(kmlString) diff --git a/tools/preprocessing/prealCF2Case.py b/tools/preprocessing/prealCF2Case.py index c95dab3..3f882f8 100644 --- a/tools/preprocessing/prealCF2Case.py +++ b/tools/preprocessing/prealCF2Case.py @@ -18,6 +18,8 @@ import netCDF4 as nc4 from PIL import Image from .ffToGeoJson import get_WSEN_LBRT_ZS_From_Pgd + + def addFieldToNcFile(ncfile, field, fieldname, typeName, dvartype): print("adding ", fieldname) @@ -63,71 +65,7 @@ def addFieldToNcFile(ncfile, field, fieldname, typeName, dvartype): variable.type = typeName; return variable - - -def FiretoNC(filename, domainProperties, parametersProperties, fuelModelMap, elevation=None, wind=None, fluxModelMap =None, bmap=None, cellMap=None): - - ncfile = nc4.Dataset(filename, 'w', format='NETCDF3_CLASSIC') - ncfile.version = "FF.1.0" - domain = ncfile.createVariable('domain', 'S1', ()) - domain.type = "domain" - for key, value in domainProperties.items(): - setattr(domain, key, value) - - parameters = ncfile.createVariable('parameters', 'S1', ()) - parameters.type = "parameters" - - if (parametersProperties is not None): - for key, value in parametersProperties.items(): - setattr(parameters, key, value) - - - if (fuelModelMap is not None): - addFieldToNcFile(ncfile, fuelModelMap, 'fuel', 'fuel', 'i4') - - if (elevation is not None): - addFieldToNcFile(ncfile, elevation, 'altitude', 'data', 'f4') - - if (wind is not None): - addFieldToNcFile(ncfile, wind, 'wind', 'data', 'f4') - # 2 types de vents, - # soit windU/windV statique (2 arrays, windU et WindV de forme: 1,1,NI,NJ) - # soit variable dans le temps, dynamique 2 arrays, windU et WindV de forme NT, 1, NI,NJ) - # soit champ potentiel nom wind de forme NAXES,NDIRCOEFFS,NI,NJ -# addFieldToNcFile(ncfile, wind["zonal"], 'windU', 'data', 'f4') - # addFieldToNcFile(ncfile, wind["meridian"], 'windV', 'data', 'f4') - - # windU = ncfile.createVariable('windU', 'f8', ('DIMT', 'DIMZ', 'DIMY', 'DIMX')) - # windU.type = "data" - # windU[0,0,:,:] = wind["zonal"] - # windV = ncfile.createVariable('windV', 'f8', ('DIMT', 'DIMZ', 'DIMY', 'DIMX')) - # windV.type = "data" - # windV[0,0,:,:] = wind["meridian"] - - if (fluxModelMap is not None): - numOfModels = 0 - for fMap in fluxModelMap: - numOfModels += len(fMap["table"]) - - for fMap in fluxModelMap: - fVar = addFieldToNcFile(ncfile, fMap["data"], fMap["name"], 'flux', 'i4') - #ncfile.createVariable(fMap["name"], 'i4', ('DIMT', 'DIMZ', 'DIMY', 'DIMX')) - #fVar.type = "flux" ; - for entry in fMap["table"].keys(): - setattr(fVar, "model%dname"%fMap["table"][entry], entry) - fVar.indices = np.array(list(fMap["table"].values()),dtype=('i4')) - #fVar[0,0,:,:] = fMap["data"] - - print("writing ", filename) - ncfile.sync() - ncfile.close() - -def PGD2Case(pgd_path, png_path, out_path, dateStartDom, fuel_test=None, gen_wind=None): - WSEN, LBRT, ZS = get_WSEN_LBRT_ZS_From_Pgd(pgd_path) - image2Case(WSEN, LBRT, ZS, png_path, out_path, dateStartDom, fuel_test=fuel_test, gen_wind=None) - - def image2Case(WSEN, LBRT, ZS, png_path, out_path, dateStartDom, fuel_test=None, gen_wind=None): @@ -198,12 +136,7 @@ def image2Case(WSEN, LBRT, ZS, png_path, out_path, dateStartDom, fuel_test=None, ni=3 nj=3 wind=np.full((2, 4, ni, nj), 0) - #wind params examples : - #Southerly 10m.s-1 trigger[wind;loc=(0.,0.,0.);vel=(0.0,10.0,0.);t=0.] - #Westerly 10m.s-1 trigger[wind;loc=(0.,0.,0.);vel=(10.0,0.0,0.);t=0.] - #Northerly 10m.s-1 trigger[wind;loc=(0.,0.,0.);vel=(0.0,-10.0,0.);t=0.] - #Easterly 10m.s-1 trigger[wind;loc=(0.,0.,0.);vel=(-10.0,0.0,0.);t=0.] - N = 8 # Nombre de directions + N = é # Nombre de directions angle_step = 2 * np.pi / N wind = np.zeros((2, N, ni, nj)) @@ -214,24 +147,7 @@ def image2Case(WSEN, LBRT, ZS, png_path, out_path, dateStartDom, fuel_test=None, v = 10 * np.sin(angle) # Composante V wind[0][i] = np.full((ni, nj), u) wind[1][i] = np.full((ni, nj), v) - # #for U - # #stherly - # wind[0][0] = np.full((ni, nj), 0) - # #westerly - # wind[0][1] = np.full((ni, nj), 10) - # #northrly - # wind[0][2] = np.full((ni, nj), 0) - # #eastrly - # wind[0][3] = np.full((ni, nj), -10) - - # #for V - # wind[1][0] = np.full((ni, nj), 10) - # #westerly - # wind[1][1] = np.full((ni, nj), 0) - # #northrly - # wind[1][2] = np.full((ni, nj), -10) - # #eastrly - # wind[1][3] = np.full((ni, nj), 0) + heatFluxModelMap = {} @@ -247,3 +163,67 @@ def image2Case(WSEN, LBRT, ZS, png_path, out_path, dateStartDom, fuel_test=None, FiretoNC(fout, domainProperties,parametersProperties,fuelMap,elevation=elevation, wind=wind, fluxModelMap = (heatFluxModelMap,vaporFluxModelMap)) + +def FiretoNC(filename, domainProperties, parametersProperties, fuelModelMap, elevation=None, wind=None, fluxModelMap =None, bmap=None, cellMap=None): + + ncfile = nc4.Dataset(filename, 'w', format='NETCDF3_CLASSIC') + ncfile.version = "FF.1.0" + domain = ncfile.createVariable('domain', 'S1', ()) + domain.type = "domain" + for key, value in domainProperties.items(): + setattr(domain, key, value) + + parameters = ncfile.createVariable('parameters', 'S1', ()) + parameters.type = "parameters" + + if (parametersProperties is not None): + for key, value in parametersProperties.items(): + setattr(parameters, key, value) + + + if (fuelModelMap is not None): + addFieldToNcFile(ncfile, fuelModelMap, 'fuel', 'fuel', 'i4') + + if (elevation is not None): + addFieldToNcFile(ncfile, elevation, 'altitude', 'data', 'f4') + + if (wind is not None): + addFieldToNcFile(ncfile, wind, 'wind', 'data', 'f4') + # 2 types de vents, + # soit windU/windV statique (2 arrays, windU et WindV de forme: 1,1,NI,NJ) + # soit variable dans le temps, dynamique 2 arrays, windU et WindV de forme NT, 1, NI,NJ) + # soit champ potentiel nom wind de forme NAXES,NDIRCOEFFS,NI,NJ +# addFieldToNcFile(ncfile, wind["zonal"], 'windU', 'data', 'f4') + # addFieldToNcFile(ncfile, wind["meridian"], 'windV', 'data', 'f4') + + # windU = ncfile.createVariable('windU', 'f8', ('DIMT', 'DIMZ', 'DIMY', 'DIMX')) + # windU.type = "data" + # windU[0,0,:,:] = wind["zonal"] + # windV = ncfile.createVariable('windV', 'f8', ('DIMT', 'DIMZ', 'DIMY', 'DIMX')) + # windV.type = "data" + # windV[0,0,:,:] = wind["meridian"] + + if (fluxModelMap is not None): + numOfModels = 0 + for fMap in fluxModelMap: + numOfModels += len(fMap["table"]) + + for fMap in fluxModelMap: + fVar = addFieldToNcFile(ncfile, fMap["data"], fMap["name"], 'flux', 'i4') + #ncfile.createVariable(fMap["name"], 'i4', ('DIMT', 'DIMZ', 'DIMY', 'DIMX')) + #fVar.type = "flux" ; + for entry in fMap["table"].keys(): + setattr(fVar, "model%dname"%fMap["table"][entry], entry) + fVar.indices = np.array(list(fMap["table"].values()),dtype=('i4')) + #fVar[0,0,:,:] = fMap["data"] + + + print("writing ", filename) + ncfile.sync() + ncfile.close() + +def PGD2Case(pgd_path, png_path, out_path, dateStartDom, fuel_test=None, gen_wind=None): + WSEN, LBRT, ZS = get_WSEN_LBRT_ZS_From_Pgd(pgd_path) + image2Case(WSEN, LBRT, ZS, png_path, out_path, dateStartDom, fuel_test=fuel_test, gen_wind=None) + + diff --git a/tools/preprocessing/upDirHDR.py b/tools/preprocessing/upDirHDR.py index 2cf5e1d..e8d16c7 100644 --- a/tools/preprocessing/upDirHDR.py +++ b/tools/preprocessing/upDirHDR.py @@ -126,7 +126,7 @@ def makeSubset(): data_xarray = dirHDRtoXarray(dirin, hdrin) print("data loaded") # Sélectionner le sous-ensemble - subset = data_xarray.sel(latitude=slice(30, 60.02), longitude=slice(-15, 35)) + subset = data_xarray.sel(latitude=slice(28, 60.02), longitude=slice(-15, 40)) #subset.plot() # Augmenter la résolution # Le facteur 4 signifie que chaque dimension sera 4 fois plus grande @@ -159,3 +159,4 @@ def makeSubset(): # Sauvegarder les données xarrayToDirHdr(resampled_xarray, dirout, hdrout) +makeSubset() \ No newline at end of file