From 04d219eef46f52d961f7211c139056f4b1b08ae7 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Tue, 23 Jul 2019 14:52:36 +0200 Subject: [PATCH 1/3] - Agriculture: use global tiff to make the land use auxiliary file --- hermesv3_bu/io_server/io_raster.py | 19 +++++++++------ hermesv3_bu/sectors/agricultural_sector.py | 27 +++++++++++++++++++++- 2 files changed, 38 insertions(+), 8 deletions(-) diff --git a/hermesv3_bu/io_server/io_raster.py b/hermesv3_bu/io_server/io_raster.py index 562921f..d004468 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 b257406..cb8fbac 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) -- GitLab From ac3541981bca1fdf53668c5cdb384c49fff9a22f Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Tue, 23 Jul 2019 15:00:21 +0200 Subject: [PATCH 2/3] PEP8 corrections --- hermesv3_bu/sectors/aviation_sector.py | 10 ++++++---- hermesv3_bu/sectors/point_source_sector.py | 1 - hermesv3_bu/sectors/recreational_boats_sector.py | 3 ++- hermesv3_bu/writer/default_writer.py | 1 + 4 files changed, 9 insertions(+), 6 deletions(-) diff --git a/hermesv3_bu/sectors/aviation_sector.py b/hermesv3_bu/sectors/aviation_sector.py index 4987cd7..003c3a3 100755 --- a/hermesv3_bu/sectors/aviation_sector.py +++ b/hermesv3_bu/sectors/aviation_sector.py @@ -411,7 +411,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 +472,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) diff --git a/hermesv3_bu/sectors/point_source_sector.py b/hermesv3_bu/sectors/point_source_sector.py index cf4fd45..d313d06 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 de75925..59177fc 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 91773f5..8076de5 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): -- GitLab From 75ddbd38473ce538e2d678738f36d5fd3ab7d5c4 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Tue, 23 Jul 2019 15:00:57 +0200 Subject: [PATCH 3/3] PEP8 corrections --- conf/hermes.conf | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/conf/hermes.conf b/conf/hermes.conf index 05ee139..a2bd0a3 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 -- GitLab