From fa0b502b72fc5f7ef28e0c9bc4bd00fafddfd931 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Fri, 17 Jan 2020 10:28:09 +0100 Subject: [PATCH 1/4] Solved bug on Agricultural Machinery sector. Solved bug on MONARCH writer --- hermesv3_bu/grids/grid.py | 2 +- .../sectors/agricultural_machinery_sector.py | 20 +++++++++---------- hermesv3_bu/writer/monarch_writer.py | 16 ++++++++++----- 3 files changed, 22 insertions(+), 16 deletions(-) diff --git a/hermesv3_bu/grids/grid.py b/hermesv3_bu/grids/grid.py index 6dea082..1f1e49b 100755 --- a/hermesv3_bu/grids/grid.py +++ b/hermesv3_bu/grids/grid.py @@ -243,5 +243,5 @@ class Grid(object): from cdo import Cdo cdo = Cdo() - cell_area = cdo.gridarea(input=input_file, returnArray='cell_area') + cell_area = cdo.gridarea(input=self.netcdf_path, returnArray='cell_area') self.shapefile['cell_area'] = cell_area.flatten() diff --git a/hermesv3_bu/sectors/agricultural_machinery_sector.py b/hermesv3_bu/sectors/agricultural_machinery_sector.py index 6d62458..6021ecd 100755 --- a/hermesv3_bu/sectors/agricultural_machinery_sector.py +++ b/hermesv3_bu/sectors/agricultural_machinery_sector.py @@ -11,7 +11,7 @@ import numpy as np from hermesv3_bu.sectors.agricultural_sector import AgriculturalSector from hermesv3_bu.io_server.io_shapefile import IoShapefile -from hermesv3_bu.tools.checker import check_files +from hermesv3_bu.tools.checker import check_files, error_exit class AgriculturalMachinerySector(AgriculturalSector): @@ -233,6 +233,7 @@ class AgriculturalMachinerySector(AgriculturalSector): return df.loc[:, ['MF']] # month_distribution = self.crop_distribution.loc[:, ['FID', 'timezone', 'geometry']].copy() dataframe = self.calcualte_yearly_emissions_by_nut_vehicle().reset_index() + dataframe['MF'] = dataframe.groupby('vehicle').apply( lambda x: get_mf(x, month) ) @@ -253,19 +254,19 @@ class AgriculturalMachinerySector(AgriculturalSector): aux = df.apply(lambda row: row * nut_emissions) return aux.loc[:, self.source_pollutants] - self.crop_distribution.reset_index(inplace=True) - self.crop_distribution[self.source_pollutants] = self.crop_distribution.groupby('NUT_code')['fraction'].apply( + crop_distribution = self.crop_distribution.reset_index().copy() + crop_distribution[self.source_pollutants] = crop_distribution.groupby('NUT_code')['fraction'].apply( lambda x: distribute_by_nut(x, dataframe.loc[int(x.name), self.source_pollutants]) ) - self.crop_distribution.drop(columns=['fraction', 'NUT_code'], inplace=True) - timezones = self.crop_distribution.groupby('FID')[['timezone']].first() - self.crop_distribution = self.crop_distribution.reset_index().groupby('FID').sum() + crop_distribution.drop(columns=['fraction', 'NUT_code'], inplace=True) + timezones = crop_distribution.groupby('FID')[['timezone']].first() + crop_distribution = crop_distribution.reset_index().groupby('FID').sum() - self.crop_distribution['timezone'] = timezones - self.crop_distribution.reset_index(inplace=True) + crop_distribution['timezone'] = timezones + crop_distribution.reset_index(inplace=True) self.logger.write_time_log('AgriculturalMachinerySector', 'distribute', timeit.default_timer() - spent_time) - return self.crop_distribution + return crop_distribution def add_dates(self, df_by_month): spent_time = timeit.default_timer() @@ -343,7 +344,6 @@ class AgriculturalMachinerySector(AgriculturalSector): for month in self.months.keys(): distribution_by_month[month] = self.calculate_monthly_emissions_by_nut(month) distribution_by_month[month] = self.distribute(distribution_by_month[month]) - self.crop_distribution = self.add_dates(distribution_by_month) self.crop_distribution.drop('date_utc', axis=1, inplace=True) self.crop_distribution = self.calculate_hourly_emissions() diff --git a/hermesv3_bu/writer/monarch_writer.py b/hermesv3_bu/writer/monarch_writer.py index 3e43f0d..270ddaf 100755 --- a/hermesv3_bu/writer/monarch_writer.py +++ b/hermesv3_bu/writer/monarch_writer.py @@ -95,11 +95,14 @@ class MonarchWriter(Writer): # From mol/h g/h to mol/m2.s g/m2.s emissions = emissions.divide(cell_area['cell_area'].mul(3600), axis=0, level='FID') - + print(self.pollutant_info) for pollutant, info in self.pollutant_info.iterrows(): if info.get('units') == "kg.s-1.m-2": - # From g.s-1.m-2 to kg.s-1.m-2 - emissions[[pollutant]] = emissions[[pollutant]].div(10**3) + try: + # From g.s-1.m-2 to kg.s-1.m-2 + emissions[pollutant] = emissions[pollutant].div(10**3) + except KeyError: + pass self.logger.write_time_log('MonarchWriter', '__init__', timeit.default_timer() - spent_time) return emissions @@ -188,7 +191,7 @@ class MonarchWriter(Writer): rlon[:] = self.grid.rlon # ========== POLLUTANTS ========== - for var_name in emissions.columns.values: + for var_name in self.pollutant_info.index: self.logger.write_log('\t\tCreating {0} variable'.format(var_name), message_level=3) # var = netcdf.createVariable(var_name, np.float64, ('time', 'lev',) + var_dim, @@ -200,7 +203,10 @@ class MonarchWriter(Writer): else: var = netcdf.createVariable(var_name, np.float64, ('time', 'lev',) + var_dim, zlib=True) - var_data = self.dataframe_to_array(emissions.loc[:, [var_name]]) + try: + var_data = self.dataframe_to_array(emissions.loc[:, [var_name]]) + except KeyError: + var_data = 0 var[:, :, self.rank_distribution[self.comm_write.Get_rank()]['y_min']: -- GitLab From e73f0b9d4d9a7852c3ea351635911cb53b64c2c4 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Mon, 20 Jan 2020 17:45:53 +0100 Subject: [PATCH 2/4] Solved bugs on rotated domains for MONARCH --- hermesv3_bu/io_server/io_netcdf.py | 4 ++-- hermesv3_bu/sectors/traffic_area_sector.py | 8 ++++---- hermesv3_bu/writer/monarch_writer.py | 2 +- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/hermesv3_bu/io_server/io_netcdf.py b/hermesv3_bu/io_server/io_netcdf.py index ed44f4a..56f57e8 100755 --- a/hermesv3_bu/io_server/io_netcdf.py +++ b/hermesv3_bu/io_server/io_netcdf.py @@ -150,12 +150,12 @@ class IoNetcdf(IoServer): try: var = nc.variables[var_name][i_time:i_time + (len(date_array)), j_min:j_max, i_min:i_max] except KeyError as e: - error_exit("{0} variable not found in {1} file.".format(str(e), path)) + error_exit("{0} variable not founddate_array in {1} file.".format(str(e), path)) nc.close() # That condition is fot the cases that the needed temperature is in a different NetCDF. while len(var) < len(date_array): - aux_date = date_array[len(var) + 1] + aux_date = date_array[len(var)] path = os.path.join(netcdf_dir, '{0}_{1}{2}.nc'.format(var_name, aux_date.year, str(aux_date.month).zfill(2))) # self.logger.write_log('Getting {0} from {1}'.format(var_name, path), message_level=2) diff --git a/hermesv3_bu/sectors/traffic_area_sector.py b/hermesv3_bu/sectors/traffic_area_sector.py index 44e9ae0..db74dc9 100755 --- a/hermesv3_bu/sectors/traffic_area_sector.py +++ b/hermesv3_bu/sectors/traffic_area_sector.py @@ -132,7 +132,7 @@ class TrafficAreaSector(Sector): self.logger.write_log("\t\tCalculating population percentage on source resolution", message_level=3) pop_by_nut2 = self.get_population_by_nut2(pop_by_nut_path) pop_shp['tot_pop'] = pop_shp['nut_code'].map(pop_by_nut2) - pop_shp['pop_percent'] = pop_shp['population'] / pop_shp['tot_pop'] + pop_shp['pop_per'] = pop_shp['population'] / pop_shp['tot_pop'] pop_shp.drop(columns=['tot_pop', 'population'], inplace=True) # 5th Calculate percent by destiny cell @@ -142,7 +142,7 @@ class TrafficAreaSector(Sector): pop_shp = self.spatial_overlays(pop_shp.reset_index(), self.grid.shapefile.reset_index()) pop_shp.drop(columns=['idx1', 'idx2', 'index'], inplace=True) pop_shp['src_inter_fraction'] = pop_shp.geometry.area / pop_shp['src_inter_fraction'] - pop_shp['pop_percent'] = pop_shp['pop_percent'] * pop_shp['src_inter_fraction'] + pop_shp['pop_per'] = pop_shp['pop_per'] * pop_shp['src_inter_fraction'] pop_shp.drop(columns=['src_inter_fraction'], inplace=True) pop_shp = IoShapefile(self.comm).gather_shapefile(pop_shp) @@ -280,8 +280,8 @@ class TrafficAreaSector(Sector): left_on='nut_code', right_on='nuts3_id', how='left') vehicle_by_cell.drop(columns=['nut_code', 'nuts3_id'], inplace=True) vehicle_by_cell[vehicle_list] = vehicle_by_cell[vehicle_list].multiply( - vehicle_by_cell['pop_percent'], axis='index') - vehicle_by_cell.drop(columns=['pop_percent'], inplace=True) + vehicle_by_cell['pop_per'], axis='index') + vehicle_by_cell.drop(columns=['pop_per'], inplace=True) vehicle_by_cell = IoShapefile(self.comm).gather_shapefile(vehicle_by_cell, rank=0) if self.comm.Get_rank() == 0: vehicle_by_cell = vehicle_by_cell.groupby('FID').sum() diff --git a/hermesv3_bu/writer/monarch_writer.py b/hermesv3_bu/writer/monarch_writer.py index 270ddaf..deb5d4c 100755 --- a/hermesv3_bu/writer/monarch_writer.py +++ b/hermesv3_bu/writer/monarch_writer.py @@ -93,9 +93,9 @@ class MonarchWriter(Writer): cell_area = None cell_area = self.comm_write.bcast(cell_area, root=0) + emissions = emissions.reset_index().groupby(['FID', 'layer', 'tstep']).sum() # From mol/h g/h to mol/m2.s g/m2.s emissions = emissions.divide(cell_area['cell_area'].mul(3600), axis=0, level='FID') - print(self.pollutant_info) for pollutant, info in self.pollutant_info.iterrows(): if info.get('units') == "kg.s-1.m-2": try: -- GitLab From 82c68b9cc309e514fcb69bddbd3ed65c8fbcd0be Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Wed, 12 Feb 2020 13:18:49 +0100 Subject: [PATCH 3/4] Intersected clip with the unary union of the domain --- hermesv3_bu/clipping/clip.py | 7 ++++--- hermesv3_bu/clipping/custom_clip.py | 14 ++++++++++---- hermesv3_bu/clipping/default_clip.py | 2 +- hermesv3_bu/clipping/shapefile_clip.py | 8 +++++--- 4 files changed, 20 insertions(+), 11 deletions(-) diff --git a/hermesv3_bu/clipping/clip.py b/hermesv3_bu/clipping/clip.py index 5b2ebb6..9addce1 100755 --- a/hermesv3_bu/clipping/clip.py +++ b/hermesv3_bu/clipping/clip.py @@ -34,10 +34,10 @@ def select_clip(comm, logger, auxiliary_path, clipping, grid): clip = DefaultClip(logger, auxiliary_path, grid) elif clipping[0] == os.path.sep: from hermesv3_bu.clipping.shapefile_clip import ShapefileClip - clip = ShapefileClip(logger, auxiliary_path, clipping) + clip = ShapefileClip(logger, auxiliary_path, clipping, grid) else: from hermesv3_bu.clipping.custom_clip import CustomClip - clip = CustomClip(logger, auxiliary_path, clipping) + clip = CustomClip(logger, auxiliary_path, clipping, grid) else: clip = None @@ -49,9 +49,10 @@ def select_clip(comm, logger, auxiliary_path, clipping, grid): class Clip(object): - def __init__(self, logger, auxiliary_path): + def __init__(self, logger, auxiliary_path, grid): spent_time = timeit.default_timer() self.logger = logger + self.grid = grid self.shapefile = None self.shapefile_path = os.path.join(auxiliary_path, 'clip', 'clip.shp') diff --git a/hermesv3_bu/clipping/custom_clip.py b/hermesv3_bu/clipping/custom_clip.py index f6f79b2..2026cfb 100755 --- a/hermesv3_bu/clipping/custom_clip.py +++ b/hermesv3_bu/clipping/custom_clip.py @@ -9,7 +9,7 @@ from hermesv3_bu.logger.log import Log class CustomClip(Clip): - def __init__(self, logger, auxiliary_path, points_str): + def __init__(self, logger, auxiliary_path, points_str, grid): """ Initialise the Custom Clip class @@ -24,7 +24,7 @@ class CustomClip(Clip): """ spent_time = timeit.default_timer() logger.write_log('Custom clip selected') - super(CustomClip, self).__init__(logger, auxiliary_path) + super(CustomClip, self).__init__(logger, auxiliary_path, grid) self.clip_type = 'Custom clip' self.shapefile = self.create_clip(points_str) self.logger.write_time_log('CustomClip', '__init__', timeit.default_timer() - spent_time) @@ -56,11 +56,17 @@ class CustomClip(Clip): if not ((lon_list[0] == lon_list[-1]) and (lat_list[0] == lat_list[-1])): lon_list.append(lon_list[0]) lat_list.append(lat_list[0]) - + geom = Polygon([[p.x, p.y] for p in [Point(xy) for xy in zip(lon_list, lat_list)]]) clip = gpd.GeoDataFrame( - geometry=[Polygon([[p.x, p.y] for p in [Point(xy) for xy in zip(lon_list, lat_list)]])], + geometry=[geom], crs={'init': 'epsg:4326'}) + border = gpd.GeoDataFrame(geometry=[self.grid.shapefile.unary_union], crs=self.grid.shapefile.crs) + geom = gpd.overlay(clip, border.to_crs(clip.crs), how='intersection').unary_union + clip = gpd.GeoDataFrame( + geometry=[geom], + crs={'init': 'epsg:4326'}) + clip.to_file(self.shapefile_path) else: clip = gpd.read_file(self.shapefile_path) diff --git a/hermesv3_bu/clipping/default_clip.py b/hermesv3_bu/clipping/default_clip.py index f05bda8..5cc1f99 100755 --- a/hermesv3_bu/clipping/default_clip.py +++ b/hermesv3_bu/clipping/default_clip.py @@ -24,7 +24,7 @@ class DefaultClip(Clip): """ spent_time = timeit.default_timer() logger.write_log('Default clip selected') - super(DefaultClip, self).__init__(logger, auxiliary_path) + super(DefaultClip, self).__init__(logger, auxiliary_path, grid) self.clip_type = 'Default clip' self.shapefile = self.create_clip(grid) self.logger.write_time_log('DefaultClip', '__init__', timeit.default_timer() - spent_time) diff --git a/hermesv3_bu/clipping/shapefile_clip.py b/hermesv3_bu/clipping/shapefile_clip.py index 88792ef..7d4eb46 100755 --- a/hermesv3_bu/clipping/shapefile_clip.py +++ b/hermesv3_bu/clipping/shapefile_clip.py @@ -10,7 +10,7 @@ from hermesv3_bu.tools.checker import error_exit class ShapefileClip(Clip): - def __init__(self, logger, auxiliary_path, clip_input_path): + def __init__(self, logger, auxiliary_path, clip_input_path, grid): """ Initialise the Shapefile Clip class @@ -25,7 +25,7 @@ class ShapefileClip(Clip): """ spent_time = timeit.default_timer() logger.write_log('Shapefile clip selected') - super(ShapefileClip, self).__init__(logger, auxiliary_path) + super(ShapefileClip, self).__init__(logger, auxiliary_path, grid) self.clip_type = 'Shapefile clip' self.shapefile = self.create_clip(clip_input_path) self.logger.write_time_log('ShapefileClip', '__init__', timeit.default_timer() - spent_time) @@ -46,7 +46,9 @@ class ShapefileClip(Clip): if not os.path.exists(os.path.dirname(self.shapefile_path)): os.makedirs(os.path.dirname(self.shapefile_path)) clip = gpd.read_file(clip_path) - clip = gpd.GeoDataFrame(geometry=[clip.unary_union], crs=clip.crs) + border = gpd.GeoDataFrame(geometry=[self.grid.shapefile.unary_union], crs=self.grid.shapefile.crs) + geom = gpd.overlay(clip, border.to_crs(clip.crs), how='intersection').unary_union + clip = gpd.GeoDataFrame(geometry=[geom], crs=clip.crs) clip.to_file(self.shapefile_path) else: error_exit(" Clip shapefile {0} not found.") -- GitLab From 11ddd1805ab50454dbc1e91670ee42272870df75 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Wed, 12 Feb 2020 13:33:34 +0100 Subject: [PATCH 4/4] PEP8 corrections --- hermesv3_bu/clipping/custom_clip.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hermesv3_bu/clipping/custom_clip.py b/hermesv3_bu/clipping/custom_clip.py index 2026cfb..e780dd2 100755 --- a/hermesv3_bu/clipping/custom_clip.py +++ b/hermesv3_bu/clipping/custom_clip.py @@ -66,7 +66,7 @@ class CustomClip(Clip): clip = gpd.GeoDataFrame( geometry=[geom], crs={'init': 'epsg:4326'}) - + clip.to_file(self.shapefile_path) else: clip = gpd.read_file(self.shapefile_path) -- GitLab