diff --git a/.gitignore b/.gitignore index 723ef36f4e4f32c4560383aa5987c575a30c6535..2849b446a4e038c0ba9ff3ff742a6c84c1aaffa4 100755 --- a/.gitignore +++ b/.gitignore @@ -1 +1,3 @@ -.idea \ No newline at end of file +.idea +.directory +*.pyc \ No newline at end of file diff --git a/conf/hermes.conf b/conf/hermes.conf index 05ee139264a51505d7a1a690a4436d2a0012a4b7..a2bd0a3acdeb2084405b8743f04a5f50f5c76c5e 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -3,13 +3,13 @@ log_level = 3 input_dir = /home/Earth/ctena/Models/hermesv3_bu_data data_path = /esarchive/recon output_dir = /scratch/Earth/HERMESv3_BU_OUT -output_name = HERMES__ships_new_2.nc +output_name = HERMES_.nc emission_summary = 0 start_date = 2016/11/29 00:00:00 # ----- end_date = start_date [DEFAULT] ----- # end_date = 2010/01/01 00:00:00 output_timestep_num = 25 -auxiliary_files_path = /scratch/Earth/HERMESv3_BU_aux/__test +auxiliary_files_path = /scratch/Earth/HERMESv3_BU_aux/_ erase_auxiliary_files = 0 @@ -115,10 +115,10 @@ vertical_description = /profiles/vertical/MONARCH_Global_48layers_ver # SECTORS #################################################################### [SECTOR MANAGEMENT] -writing_processors = 4 +writing_processors = 1 aviation_processors = 0 -shipping_port_processors = 2 +shipping_port_processors = 0 livestock_processors = 0 crop_operations_processors = 0 crop_fertilizers_processors = 0 @@ -204,7 +204,7 @@ livestock_hourly_profiles = /profiles/temporal/livestock/hourly_profi livestock_speciation_profiles = /profiles/speciation/livestock/speciation_profiles_base.csv [AGRICULTURAL] -land_uses_path = /ecmwf/clc/original_files/g250_clc12_v18_5a/g250_clc12_v18_5a_ES_CCAA.shp +land_uses_path = /ecmwf/clc/original_files/g250_clc12_v18_5a/g250_clc12_V18_5.tif land_use_by_nut_path = /agriculture/land_use_ccaa.csv crop_by_nut_path = /agriculture/crops_ha_2017.csv crop_from_landuse_path = /agriculture/map_crops_landuse.csv diff --git a/hermesv3_bu/io_server/io_raster.py b/hermesv3_bu/io_server/io_raster.py index 562921f62588b6cba3e891ff999af1b76a053e08..d00446821df997bfdd9c6452f2926f7723cac0b7 100755 --- a/hermesv3_bu/io_server/io_raster.py +++ b/hermesv3_bu/io_server/io_raster.py @@ -22,7 +22,7 @@ class IoRaster(IoServer): comm = MPI.COMM_WORLD super(IoRaster, self).__init__(comm) - def clip_raster_with_shapefile(self, raster_path, shape_path, clipped_raster_path): + def clip_raster_with_shapefile(self, raster_path, shape_path, clipped_raster_path, values=None, nodata=0): """ Clip a raster using given shapefile path. @@ -37,8 +37,8 @@ class IoRaster(IoServer): :param clipped_raster_path: Place to store the clipped raster. :type clipped_raster_path: str - :param rank: Rank who have to do the work. Default 0 - :type rank: int + :param values: List of data values to clip. + :type values: list :return: Path where is stored the clipped raster. :rtype: str @@ -57,7 +57,10 @@ class IoRaster(IoServer): geo = geo.to_crs(crs=data.crs.data) coords = getFeatures(geo) - out_img, out_transform = mask(data, shapes=coords, crop=True) + out_img, out_transform = mask(data, shapes=coords, crop=True, all_touched=True, nodata=nodata) + if values is not None: + out_img[~np.isin(out_img, values)] = nodata + out_meta = data.meta.copy() out_meta.update({ @@ -73,7 +76,7 @@ class IoRaster(IoServer): return clipped_raster_path - def clip_raster_with_shapefile_poly(self, raster_path, geo, clipped_raster_path, nodata=0): + def clip_raster_with_shapefile_poly(self, raster_path, geo, clipped_raster_path, values=None, nodata=0): """ Clip a raster using given shapefile. @@ -88,8 +91,8 @@ class IoRaster(IoServer): :param clipped_raster_path: Place to store the clipped raster. :type clipped_raster_path: str - :param rank: Rank who have to do the work. Default 0 - :type rank: int + :param values: List of data values to clip. + :type values: list :param nodata: Value for the no data elements. Default 0 :type nodata: float @@ -112,6 +115,8 @@ class IoRaster(IoServer): coords = get_features(geo) out_img, out_transform = mask(data, shapes=coords, crop=True, all_touched=True, nodata=nodata) + if values is not None: + out_img[~np.isin(out_img, values)] = nodata out_meta = data.meta.copy() out_meta.update( diff --git a/hermesv3_bu/sectors/agricultural_sector.py b/hermesv3_bu/sectors/agricultural_sector.py index b257406d61853668d6a28d5010ea7a08ac9995a6..cb8fbac85ff4f5274a39d343e0dd0b5227a7f5cf 100755 --- a/hermesv3_bu/sectors/agricultural_sector.py +++ b/hermesv3_bu/sectors/agricultural_sector.py @@ -11,6 +11,7 @@ from mpi4py import MPI from hermesv3_bu.sectors.sector import Sector from hermesv3_bu.io_server.io_shapefile import IoShapefile +from hermesv3_bu.io_server.io_raster import IoRaster from hermesv3_bu.logger.log import Log @@ -189,7 +190,7 @@ class AgriculturalSector(Sector): return land_uses_list - def get_land_use_src_by_nut(self, land_uses): + def get_land_use_src_by_nut_old(self, land_uses): spent_time = timeit.default_timer() df_land_use_with_nut = gpd.read_file(self.land_uses_path) @@ -204,6 +205,30 @@ class AgriculturalSector(Sector): self.logger.write_time_log('AgriculturalSector', 'get_land_use_src_by_nut', timeit.default_timer() - spent_time) return df_land_use_with_nut + def get_land_use_src_by_nut(self, land_uses): + spent_time = timeit.default_timer() + land_use_src_by_nut_path = os.path.join(self.auxiliary_dir, 'agriculture', 'land_uses', 'land_uses_src.shp') + if not os.path.exists(land_use_src_by_nut_path): + land_uses_clipped = IoRaster(self.comm_agr).clip_raster_with_shapefile_poly( + self.land_uses_path, self.clip.shapefile, + os.path.join(self.auxiliary_dir, 'agriculture', 'land_uses', 'land_uses_clip.tif'), values=land_uses) + + land_uses_shp = IoRaster(self.comm_agr).to_shapefile_serie(land_uses_clipped) + ccaa_shp = IoShapefile(self.comm_agr).read_shapefile_serial(self.nut_shapefile).to_crs(land_uses_shp.crs) + ccaa_shp.drop(columns=['NAME', 'ORDER06'], inplace=True) + ccaa_shp.rename(columns={'CODE': 'NUT'}, inplace=True) + land_use_src_by_nut = self.spatial_overlays(land_uses_shp, ccaa_shp, how='intersection') + land_use_src_by_nut.drop(columns=['idx1', 'idx2', 'CELL_ID'], inplace=True) + land_use_src_by_nut.rename(columns={'data': 'land_use'}, inplace=True) + land_use_src_by_nut['land_use'] = land_use_src_by_nut['land_use'].astype(np.int16) + land_use_src_by_nut.reset_index(inplace=True, drop=True) + IoShapefile(self.comm_agr).write_shapefile_serial(land_use_src_by_nut, land_use_src_by_nut_path) + else: + land_use_src_by_nut = IoShapefile(self.comm_agr).read_shapefile_serial(land_use_src_by_nut_path) + + self.logger.write_time_log('AgriculturalSector', 'get_land_use_src_by_nut', timeit.default_timer() - spent_time) + return land_use_src_by_nut + def get_tot_land_use_by_nut(self, land_uses): spent_time = timeit.default_timer() df = pd.read_csv(self.land_use_by_nut) diff --git a/hermesv3_bu/sectors/aviation_sector.py b/hermesv3_bu/sectors/aviation_sector.py index 4987cd7591078f798509b482709d556bd44874d2..a4e5f96b61907f8f0fdec6d48d0ceb2d0e52bd71 100755 --- a/hermesv3_bu/sectors/aviation_sector.py +++ b/hermesv3_bu/sectors/aviation_sector.py @@ -7,6 +7,7 @@ from hermesv3_bu.logger.log import Log import numpy as np import pandas as pd import geopandas as gpd +from warnings import warn from hermesv3_bu.sectors.sector import Sector @@ -135,12 +136,13 @@ class AviationSector(Sector): self.source_pollutants.append(poll) self.source_pollutants.remove('hc') - self.ef_dir = ef_dir + # self.ef_dir = ef_dir + self.ef_files = self.read_ef_files(ef_dir) airport_trajectories_shapefile = self.read_trajectories_shapefile( airport_trajectories_shapefile_path, airport_runways_corners_shapefile_path, airport_runways_shapefile_path) self.airport_list_full = None # only needed for master process - self.airport_list = self.get_airport_list(airport_list, airport_trajectories_shapefile) + self.airport_list = self.get_airport_list(airport_list, airport_trajectories_shapefile, operations_path) self.plane_list = plane_list full_airport_shapefile = gpd.read_file(airport_shapefile_path) @@ -164,6 +166,18 @@ class AviationSector(Sector): comm.Barrier() self.logger.write_time_log('AviationSector', '__init__', timeit.default_timer() - spent_time) + def read_ef_files(self, ef_path): + if self.comm.Get_rank() == 0: + ef_files = {} + for phase in PHASE_TYPE.keys(): + ef_files[phase] = pd.read_csv(os.path.join(ef_path, PHASE_EF_FILE[phase])) + else: + ef_files = None + + ef_files = self.comm.bcast(ef_files, root=0) + + return ef_files + def read_trajectories_shapefile(self, trajectories_path, runways_corners_path, runways_path): """ Create a shapefile with 2 geometries: trajectories & staring point @@ -279,11 +293,11 @@ class AviationSector(Sector): for index, aux_operations in operations.groupby(['airport_id', 'plane_id', 'operation']): if len(aux_operations) > 1: print index, len(aux_operations) - if self.plane_list is None: self.plane_list = list(np.unique(operations['plane_id'].values)) else: operations = operations.loc[operations['plane_id'].isin(self.plane_list), :] + if len(operations) == 0: raise NameError("The plane/s defined in the plane_list do not exist.") operations = operations.loc[operations['airport_id'].isin(self.airport_list), :] @@ -336,7 +350,7 @@ class AviationSector(Sector): return dataframe - def get_airport_list(self, conf_airport_list, airport_shapefile): + def get_airport_list(self, conf_airport_list, airport_shapefile, operations_file): """ Get the airport list from the involved airports on the domain. @@ -367,25 +381,32 @@ class AviationSector(Sector): if len(shp_airport_list) == 0: raise NameError("No airports intersect with the defined domain or the defined aiport/s in the " + "airport_list do no exist ") - max_len = len(shp_airport_list) + + airports_with_operations = np.unique(pd.read_csv(operations_file, usecols=['airport_id']).values) + + new_list = list(set(shp_airport_list) & set(airports_with_operations)) + if len(new_list) != len(shp_airport_list): + warn('{0} airports have no operations. Ignoring them.'.format( + list(set(new_list) - set(shp_airport_list)))) + + max_len = len(new_list) # Only for master (rank == 0) - self.airport_list_full = shp_airport_list + self.airport_list_full = new_list - shp_airport_list = [shp_airport_list[i * len(shp_airport_list) // self.comm.size: - (i + 1) * len(shp_airport_list) // self.comm.size] + new_list = [new_list[i * len(new_list) // self.comm.size: (i + 1) * len(new_list) // self.comm.size] for i in range(self.comm.size)] - for sublist in shp_airport_list: + for sublist in new_list: if len(sublist) == 0: raise ValueError("ERROR: The selected number of processors is to high. " + "The maximum number of processors accepted are {0}".format(max_len) + "(Maximum number of airports included in the working domain") else: - shp_airport_list = None + new_list = None - shp_airport_list = self.comm.scatter(shp_airport_list, root=0) + new_list = self.comm.scatter(new_list, root=0) self.logger.write_time_log('AviationSector', 'get_airport_list', timeit.default_timer() - spent_time) - return shp_airport_list + return new_list def calculate_airport_distribution(self, airport_shapefile): """ @@ -411,7 +432,8 @@ class AviationSector(Sector): os.makedirs(os.path.dirname(airport_distribution_path)) airport_shapefile.to_crs(self.grid_shp.crs, inplace=True) airport_shapefile['area'] = airport_shapefile.area - airport_distribution = self.spatial_overlays(airport_shapefile, self.grid_shp.reset_index(), how='intersection') + airport_distribution = self.spatial_overlays(airport_shapefile, self.grid_shp.reset_index(), + how='intersection') airport_distribution['fraction'] = airport_distribution.area / airport_distribution['area'] airport_distribution.drop(columns=['idx2', 'area', 'geometry', 'cons'], inplace=True) airport_distribution.rename(columns={'idx1': 'airport_id'}, inplace=True) @@ -471,10 +493,11 @@ class AviationSector(Sector): runway_shapefile.to_crs(self.grid_shp.crs, inplace=True) runway_shapefile['length'] = runway_shapefile.length # duplicating each runway by involved cell - runway_shapefile = gpd.sjoin(runway_shapefile, self.grid_shp.reset_index(), how="inner", op='intersects') + runway_shapefile = gpd.sjoin(runway_shapefile, self.grid_shp.reset_index(), how="inner", + op='intersects') # Adding cell geometry - runway_shapefile = runway_shapefile.merge(self.grid_shp.reset_index().loc[:, ['FID', 'geometry']], on='FID', - how='left') + runway_shapefile = runway_shapefile.merge(self.grid_shp.reset_index().loc[:, ['FID', 'geometry']], + on='FID', how='left') # Intersection between line (roadway) and polygon (cell) # runway_shapefile['geometry'] = runway_shapefile.apply(do_intersection, axis=1) runway_shapefile['mini_length'] = runway_shapefile.apply(get_intersection_length, axis=1) @@ -641,7 +664,8 @@ class AviationSector(Sector): Emission factor associated to phase and pollutant """ engine = self.planes_info.loc[df.name, 'engine_id'] - ef_dataframe = pd.read_csv(os.path.join(self.ef_dir, PHASE_EF_FILE[phase])) + # ef_dataframe = pd.read_csv(os.path.join(self.ef_dir, PHASE_EF_FILE[phase])) + ef_dataframe = self.ef_files[phase].reset_index() ef_dataframe.set_index('engine_id', inplace=True) df['EF'] = ef_dataframe.loc[engine, poll] return df.loc[:, ['EF']] @@ -680,6 +704,7 @@ class AviationSector(Sector): # Merging operations with airport geometry dataframe = pd.DataFrame(index=self.operations.xs(PHASE_TYPE[phase], level='operation').index) + dataframe = dataframe.reset_index().set_index('airport_id') dataframe = self.airport_shapefile.join(dataframe, how='inner') dataframe.index.name = 'airport_id' @@ -737,7 +762,8 @@ class AviationSector(Sector): """ Emission factor associated to phase and pollutant """ - ef_dataframe = pd.read_csv(os.path.join(self.ef_dir, PHASE_EF_FILE[phase])) + # ef_dataframe = pd.read_csv(os.path.join(self.ef_dir, PHASE_EF_FILE[phase])) + ef_dataframe = self.ef_files[phase].reset_index() ef_dataframe.set_index('plane_id', inplace=True) ef = ef_dataframe.loc['default', poll] return ef @@ -837,7 +863,8 @@ class AviationSector(Sector): Emission factor associated to phase and pollutant """ engine = self.planes_info.loc[df.name, 'apu_id'] - ef_dataframe = pd.read_csv(os.path.join(self.ef_dir, PHASE_EF_FILE[phase])) + # ef_dataframe = pd.read_csv(os.path.join(self.ef_dir, PHASE_EF_FILE[phase])) + ef_dataframe = self.ef_files[phase].reset_index() ef_dataframe.set_index('apu_id', inplace=True) try: df['EF'] = ef_dataframe.loc[engine, poll] diff --git a/hermesv3_bu/sectors/point_source_sector.py b/hermesv3_bu/sectors/point_source_sector.py index cf4fd45c6e12b07a53807f2cf91f657c9b671ef7..d313d068f532ec6b7113f1864c3e3f1b9a40144b 100755 --- a/hermesv3_bu/sectors/point_source_sector.py +++ b/hermesv3_bu/sectors/point_source_sector.py @@ -480,7 +480,6 @@ class PointSourceSector(Sector): self.plume_rise_pahts['temperature_sfc_dir'], 't2_{0}.nc'.format(self.date_array[0].replace(hour=0).strftime("%Y%m%d%H")))) - catalog = catalog.merge(meteo_xy, left_index=True, right_index=True) # ===== 3D Meteo variables ===== diff --git a/hermesv3_bu/sectors/recreational_boats_sector.py b/hermesv3_bu/sectors/recreational_boats_sector.py index de7592598999b194495384debfea30b217cd2dc0..59177fcfa3cf83a97ce0f5c2d911fd443383eb79 100755 --- a/hermesv3_bu/sectors/recreational_boats_sector.py +++ b/hermesv3_bu/sectors/recreational_boats_sector.py @@ -44,7 +44,8 @@ class RecreationalBoatsSector(Sector): src_density_map['data'] = src_density_map['data'] / src_density_map['data'].sum() src_density_map.to_crs(self.grid_shp.crs, inplace=True) src_density_map['src_inter_fraction'] = src_density_map.area - src_density_map = self.spatial_overlays(src_density_map, self.grid_shp.reset_index(), how='intersection') + src_density_map = self.spatial_overlays(src_density_map, self.grid_shp.reset_index(), + how='intersection') src_density_map['src_inter_fraction'] = src_density_map.area / src_density_map['src_inter_fraction'] src_density_map['data'] = src_density_map.loc[:, 'data'].multiply(src_density_map["src_inter_fraction"], diff --git a/hermesv3_bu/writer/default_writer.py b/hermesv3_bu/writer/default_writer.py index 91773f5851c042bdb47944511300aa2e703e9c6d..8076de535e6c966a1b14fd0a2a9cf2513611cb6b 100755 --- a/hermesv3_bu/writer/default_writer.py +++ b/hermesv3_bu/writer/default_writer.py @@ -8,6 +8,7 @@ import timeit from hermesv3_bu.logger.log import Log import time + class DefaultWriter(Writer): def __init__(self, comm_world, comm_write, logger, netcdf_path, grid, date_array, pollutant_info, rank_distribution, emission_summary=False):