From 038a9e5f06467fea2dc4e00e505f0fe4d7e1dc4f Mon Sep 17 00:00:00 2001 From: Michael Steiner Date: Wed, 12 May 2021 16:04:49 +0200 Subject: [PATCH] Add mapping for ICON model * Added class for ICON-grid in grids.py * Added class ICONGrid in grids.py * Extended the code for mapping to ICON output grid: - globally changed "cosmo_grid" to "output_grid" (except for the existing config-files) - created a new config file for the tno/ICON-case - added 2 new functions in utilities: "prepare_ICON_output_file" and "write_variable_ICON" - built in switches in __init__ and utilities where necessary to distinguish between COSMO and ICON grid - changed further variables where reasonable, e.g. cosmo_cell_x --> output_cell_x * Changed all output_grid-variables in the config-files * Gridding of the pointsources is much more efficient now looping only over cells around the closest vertex * central_longitude in Mellweide-projection is now pollon-dependent and Readme is updated with a link to the Readme in cases/ * cases/README.md is now a link --- README.md | 2 + cases/config_append.py | 4 +- cases/config_carbocount.py | 7 +- cases/config_d1_append.py | 4 +- cases/config_d1_swiss_art.py | 7 +- cases/config_d1_tno_art.py | 7 +- cases/config_edgar.py | 7 +- cases/config_icon_tno.py | 79 +++++++++++++ cases/config_tno.py | 15 +-- cases/test_swiss_cc.py | 7 +- cases/test_tno.py | 7 +- emiproc/__init__.py | 91 ++++++++------- emiproc/grids.py | 207 ++++++++++++++++++++++++++++++++++- emiproc/utilities.py | 181 ++++++++++++++++++++---------- 14 files changed, 498 insertions(+), 127 deletions(-) create mode 100644 cases/config_icon_tno.py diff --git a/README.md b/README.md index 86742a2..e9496b7 100644 --- a/README.md +++ b/README.md @@ -36,6 +36,8 @@ be obtained separately. $ python -m emiproc vp # for vertical profiles ``` +Further examples, including the generation of offline files and inventory merging, can be found in [cases/README.md](cases/README.md) + ## Gridded annual emissions Emissions are read from the inventory and projected onto the COSMO grid. diff --git a/cases/config_append.py b/cases/config_append.py index ea21c42..68a951e 100644 --- a/cases/config_append.py +++ b/cases/config_append.py @@ -1,7 +1,7 @@ import time import os -from emiproc.grids import COSMOGrid, TNOGrid +from emiproc.grids import COSMOGrid, TNOGrid, ICONGrid inv_1 = os.path.join('outputs', '{online}', 'tno.nc') inv_name_1 = '_TNO' @@ -12,7 +12,7 @@ # COSMO domain -cosmo_grid = COSMOGrid( +output_grid = COSMOGrid( nx=900, ny=600, dx=0.01, diff --git a/cases/config_carbocount.py b/cases/config_carbocount.py index 042ef54..ff6d386 100644 --- a/cases/config_carbocount.py +++ b/cases/config_carbocount.py @@ -1,12 +1,13 @@ import os import time -from emiproc.grids import COSMOGrid, SwissGrid +from emiproc.grids import COSMOGrid, SwissGrid, ICONGrid # inventory inventory = 'swiss-cc' -# model either "cosmo-art" or "cosmo-ghg" (affects the output units) +# model either "cosmo-art", "cosmo-ghg" or "icon" (affects the +# output units and handling of the output grid) model = 'cosmo-ghg' # path to input inventory @@ -85,7 +86,7 @@ varname_format = '{species}_{category}' # COSMO domain -cosmo_grid = COSMOGrid( +output_grid = COSMOGrid( nx=900, ny=600, dx=0.01, diff --git a/cases/config_d1_append.py b/cases/config_d1_append.py index de16e7d..1743b54 100644 --- a/cases/config_d1_append.py +++ b/cases/config_d1_append.py @@ -1,6 +1,6 @@ import time -from emiproc.grids import COSMOGrid, TNOGrid +from emiproc.grids import COSMOGrid, TNOGrid, ICONGrid @@ -14,7 +14,7 @@ # Output grid is European domain (rotated pole coordinates) -cosmo_grid = COSMOGrid( +output_grid = COSMOGrid( nx=192, ny=164, dx=0.12, diff --git a/cases/config_d1_swiss_art.py b/cases/config_d1_swiss_art.py index 60512fb..9c407b6 100755 --- a/cases/config_d1_swiss_art.py +++ b/cases/config_d1_swiss_art.py @@ -1,12 +1,13 @@ import os import time -from emiproc.grids import COSMOGrid, SwissGrid +from emiproc.grids import COSMOGrid, SwissGrid, ICONGrid # inventory inventory = 'swiss-art' -# model is "cosmo-art" or "cosmo-ghg" (affects the output units) +# model either "cosmo-art", "cosmo-ghg" or "icon" (affects the +# output units and handling of the output grid) model = 'cosmo-art' # add total emissions (only for swiss inventory) @@ -58,7 +59,7 @@ nx = 192 ny = 164 -cosmo_grid = COSMOGrid( +output_grid = COSMOGrid( nx=nx, ny=ny, dx=0.12, diff --git a/cases/config_d1_tno_art.py b/cases/config_d1_tno_art.py index 33c5b13..4ef7850 100644 --- a/cases/config_d1_tno_art.py +++ b/cases/config_d1_tno_art.py @@ -4,12 +4,13 @@ import os import time -from emiproc.grids import COSMOGrid, TNOGrid +from emiproc.grids import COSMOGrid, TNOGrid, ICONGrid # inventory inventory = 'TNO' -# model either "cosmo-art" or "cosmo-ghg" (affects the output units) +# model either "cosmo-art", "cosmo-ghg" or "icon" (affects the +# output units and handling of the output grid) model = 'cosmo-art' # input filename @@ -53,7 +54,7 @@ output_name = 'tno-art.nc' # Output grid is European domain (rotated pole coordinates) -cosmo_grid = COSMOGrid( +output_grid = COSMOGrid( nx=192, ny=164, dx=0.12, diff --git a/cases/config_edgar.py b/cases/config_edgar.py index 99b67eb..f6dd48f 100644 --- a/cases/config_edgar.py +++ b/cases/config_edgar.py @@ -1,11 +1,12 @@ -from emiproc.grids import COSMOGrid, EDGARGrid +from emiproc.grids import COSMOGrid, EDGARGrid, ICONGrid import os import time # inventory inventory = 'EDGAR' -# model either "cosmo-art" or "cosmo-ghg" (affect the output units) +# model either "cosmo-art", "cosmo-ghg" or "icon" (affects the +# output units and handling of the output grid) model = 'cosmo-ghg' # path to input inventory @@ -64,7 +65,7 @@ # Domain # CHE_Europe domain -cosmo_grid = COSMOGrid( +output_grid = COSMOGrid( nx=760, ny=610, dx=0.05, diff --git a/cases/config_icon_tno.py b/cases/config_icon_tno.py new file mode 100644 index 0000000..6b933c4 --- /dev/null +++ b/cases/config_icon_tno.py @@ -0,0 +1,79 @@ +import os +import time + +from emiproc.grids import COSMOGrid, TNOGrid, ICONGrid + +# inventory +inventory = 'TNO' + +# model either "cosmo-art", "cosmo-ghg" or "icon" (affects the +# output units and handling of the output grid) +model = 'icon' + +# path to input inventory +input_path = "/input/TNOMACC/TNO_GHGco/TNO_6x6_GHGco_v1_1/TNO_GHGco_v1_1_year2015.nc" + +# input grid +input_grid = TNOGrid(input_path) + +# input species +species = ['co2_ff', 'co2_bf'] + +# input categories +categories = [ + "A", + "B", + "C", + "D", + "E", + "F1", + "F2", + "F3", + "G", + "H", + "I", + "J", + "K", + "L", +] + +# mapping from input to output species (input is used for missing keys) +in2out_species = { + 'co2_ff': 'CO2', + 'co2_bf': 'CO2' +} + +# mapping from input to output species (input is used for missing keys) +in2out_category = {} + +# output variables are written in the following format using species and +# category after applying mapping as well as source_type (AREA or POINT) for +# TNO inventories +varname_format = '{species}_{category}' # not providing source_type will add up + # point and area sources + +# path to ICON output grid +icon_path = "/newhome/stem/git/C2SM-RCM/ICON_grids/MCH2km/domain1_DOM01.nc" + +# output ICON grid +output_grid = ICONGrid(icon_path) + +# output path and filename +output_path = os.path.join('outputs', '{online}') +output_name = "tno.nc" + +# resolution of shape file used for country mask +shpfile_resolution = "10m" + +# number of processes computing the mapping inventory->COSMO-grid +nprocs = 18 + +# metadata added as global attributes to netCDF output file +nc_metadata = { + "DESCRIPTION": "Gridded annual emissions", + "DATAORIGIN": "TNO", + "CREATOR": "Jean-Matthieu Haussaire", + "EMAIL": "jean-matthieu.haussaire@empa.ch", + "AFFILIATION": "Empa Duebendorf, Switzerland", + "DATE CREATED": time.ctime(time.time()), +} diff --git a/cases/config_tno.py b/cases/config_tno.py index 4591973..4d585ad 100644 --- a/cases/config_tno.py +++ b/cases/config_tno.py @@ -1,12 +1,13 @@ import os import time -from emiproc.grids import COSMOGrid, TNOGrid +from emiproc.grids import COSMOGrid, TNOGrid, ICONGrid # inventory inventory = 'TNO' -# model either "cosmo-art" or "cosmo-ghg" (affect the output units) +# model either "cosmo-art", "cosmo-ghg" or "icon" (affects the +# output units and handling of the output grid) model = 'cosmo-ghg' # path to input inventory @@ -52,11 +53,11 @@ # point and area sources # COSMO domain -cosmo_grid = COSMOGrid( - nx=900, - ny=600, - dx=0.01, - dy=0.01, +output_grid = COSMOGrid( + nx=450, #900 + ny=300, #600 + dx=0.02, #0.01 + dy=0.02, #0.01 xmin=-4.92, ymin=-3.18, pollon=-170.0, diff --git a/cases/test_swiss_cc.py b/cases/test_swiss_cc.py index 6efc96b..711538e 100644 --- a/cases/test_swiss_cc.py +++ b/cases/test_swiss_cc.py @@ -1,12 +1,13 @@ import os import time -from emiproc.grids import COSMOGrid, SwissGrid +from emiproc.grids import COSMOGrid, SwissGrid, ICONGrid # inventory inventory = 'swiss-cc' -# model either "cosmo-art" or "cosmo-ghg" (affects the output units) +# model either "cosmo-art", "cosmo-ghg" or "icon" (affects the +# output units and handling of the output grid) model = 'cosmo-ghg' # path to input inventory @@ -65,7 +66,7 @@ varname_format = '{species}_{category}' # COSMO domain -cosmo_grid = COSMOGrid( +output_grid = COSMOGrid( nx=90, ny=60, dx=0.1, diff --git a/cases/test_tno.py b/cases/test_tno.py index 9ada7f3..09a7ba6 100644 --- a/cases/test_tno.py +++ b/cases/test_tno.py @@ -1,12 +1,13 @@ import os import time -from emiproc.grids import COSMOGrid, TNOGrid +from emiproc.grids import COSMOGrid, TNOGrid, ICONGrid # inventory inventory = 'TNO' -# model either "cosmo-art" or "cosmo-ghg" (affect the output units) +# model either "cosmo-art", "cosmo-ghg" or "icon" (affects the +# output units and handling of the output grid) model = 'cosmo-ghg' # path to input inventory @@ -49,7 +50,7 @@ varname_format = '{species}_{category}_{source_type}' # COSMO domain -cosmo_grid = COSMOGrid( +output_grid = COSMOGrid( nx=90, ny=60, dx=0.1, diff --git a/emiproc/__init__.py b/emiproc/__init__.py index 3ccce41..2f4a893 100644 --- a/emiproc/__init__.py +++ b/emiproc/__init__.py @@ -24,7 +24,7 @@ def process_swiss(cfg, interpolation, country_mask, out, latname, lonname): total_flux = {} for var in cfg.species: var_name = cfg.in2out_species.get(var, var) - total_flux[var_name] = np.zeros((cfg.cosmo_grid.ny, cfg.cosmo_grid.nx)) + total_flux[var_name] = np.zeros((cfg.output_grid.ny, cfg.output_grid.nx)) for cat in cfg.categories: for var in cfg.species: @@ -44,18 +44,18 @@ def process_swiss(cfg, interpolation, country_mask, out, latname, lonname): print(filename) emi += util.read_emi_from_file(filename) # (lon,lat) - out_var = np.zeros((cfg.cosmo_grid.ny, cfg.cosmo_grid.nx)) + out_var = np.zeros((cfg.output_grid.ny, cfg.output_grid.nx)) for lon in range(np.shape(emi)[0]): for lat in range(np.shape(emi)[1]): for (x, y, r) in interpolation[lon, lat]: out_var[y, x] += emi[lon, lat] * r - cosmo_area = 1.0 / cfg.cosmo_grid.gridcell_areas() + output_area = 1.0 / cfg.output_grid.gridcell_areas() # convert units if cfg.model == 'cosmo-ghg': # COSMO-GHG: kg.year-1.cell-1 to kg.m-2.s-1 - out_var *= cosmo_area.T / util.SEC_PER_YR + out_var *= output_area.T / util.SEC_PER_YR unit = 'kg m-2 s-1' elif cfg.model == 'cosmo-art': # COSMO-ART: g.year-1.cell-1 to kg.h-1.cell-1 @@ -97,18 +97,18 @@ def process_rotgrid(cfg, interpolation, country_mask, out, latname, lonname): for var in cfg.species: print('Species', var, 'Category', cat) emi = cosmo_in[var]*(cfg.input_grid.gridcell_areas().T) - out_var = np.zeros((cfg.cosmo_grid.ny, cfg.cosmo_grid.nx)) + out_var = np.zeros((cfg.output_grid.ny, cfg.output_grid.nx)) for lon in range(np.shape(emi)[1]): for lat in range(np.shape(emi)[0]): for (x, y, r) in interpolation[lon, lat]: out_var[y, x] += emi[lat, lon] * r - cosmo_area = 1.0 / cfg.cosmo_grid.gridcell_areas() + output_area = 1.0 / cfg.output_grid.gridcell_areas() # convert units # COSMO-GHG: kg.s-1.cell-1 to kg.m-2.s-1 - out_var *= cosmo_area.T + out_var *= output_area.T unit = 'kg m-2 s-1' # only allow positive fluxes @@ -150,7 +150,7 @@ def process_tno(cfg, interpolation, country_mask, out, latname, lonname): selection_point = tno["source_type_index"][:] == 2 # Area of the COSMO grid cells - cosmo_area = 1.0 / cfg.cosmo_grid.gridcell_areas() + output_area = 1.0 / cfg.output_grid.gridcell_areas() for cat in cfg.categories: @@ -176,10 +176,10 @@ def process_tno(cfg, interpolation, country_mask, out, latname, lonname): for s in cfg.species: print("Species", s, "Category", cat) out_var_area = np.zeros( - (cfg.cosmo_grid.ny, cfg.cosmo_grid.nx) + (cfg.output_grid.ny, cfg.output_grid.nx) ) out_var_point = np.zeros( - (cfg.cosmo_grid.ny, cfg.cosmo_grid.nx) + (cfg.output_grid.ny, cfg.output_grid.nx) ) var = tno[s][:] @@ -195,7 +195,7 @@ def process_tno(cfg, interpolation, country_mask, out, latname, lonname): out_var_area[y, x] += var[i] * r if selection_cat_point[i]: try: - indx, indy = cfg.cosmo_grid.indices_of_point( + indx, indy = cfg.output_grid.indices_of_point( tno["longitude_source"][i], tno["latitude_source"][i], ) @@ -209,10 +209,10 @@ def process_tno(cfg, interpolation, country_mask, out, latname, lonname): print("Gridding took %.1f seconds" % (end - start)) # convert units - if cfg.model == 'cosmo-ghg': + if cfg.model == 'cosmo-ghg' or cfg.model == 'icon': # COSMO-GHG: kg.year-1.cell-1 to kg.m-2.s-1 - out_var_point *= cosmo_area.T / util.SEC_PER_YR - out_var_area *= cosmo_area.T / util.SEC_PER_YR + out_var_point *= output_area.T / util.SEC_PER_YR + out_var_area *= output_area.T / util.SEC_PER_YR unit = 'kg m-2 s-1' elif cfg.model == 'cosmo-art': @@ -231,10 +231,14 @@ def process_tno(cfg, interpolation, country_mask, out, latname, lonname): ): if sel.any(): out_var_name = util.get_out_varname(s, cat, cfg, - source_type=t) + source_type=t) print('Write as variable:', out_var_name) - util.write_variable(out, out_var, out_var_name, - latname, lonname, unit) + if cfg.model.startswith("cosmo"): + util.write_variable(out, out_var, out_var_name, + latname, lonname, unit) + elif cfg.model.startswith("icon"): + util.write_variable_ICON(out, out_var, + out_var_name, unit) @@ -247,44 +251,53 @@ def main(cfg): # Load or compute the country mask country_mask = util.get_country_mask( cfg.output_path, - cfg.cosmo_grid.name, - cfg.cosmo_grid, + cfg.output_grid.name, + cfg.output_grid, cfg.shpfile_resolution, cfg.nprocs, ) - # Load or compute the mapping between the inventory and COSMO grid + # Load or compute the mapping between the inventory and output grid interpolation = util.get_gridmapping( cfg.output_path, cfg.input_grid.name, - cfg.cosmo_grid, + cfg.output_grid, cfg.input_grid, cfg.nprocs, ) # Set names for longitude and latitude - if ( - cfg.cosmo_grid.pollon == 180 or cfg.cosmo_grid.pollon == 0 - ) and cfg.cosmo_grid.pollat == 90: - lonname = "lon" - latname = "lat" - print( - "Non-rotated grid: pollon = %f, pollat = %f" - % (cfg.cosmo_grid.pollon, cfg.cosmo_grid.pollat) - ) - else: - lonname = "rlon" - latname = "rlat" - print( - "Rotated grid: pollon = %f, pollat = %f" - % (cfg.cosmo_grid.pollon, cfg.cosmo_grid.pollat) - ) + if cfg.model.startswith("cosmo"): + if ( + cfg.output_grid.pollon == 180 or cfg.output_grid.pollon == 0 + ) and cfg.output_grid.pollat == 90: + lonname = "lon" + latname = "lat" + print( + "Non-rotated grid: pollon = %f, pollat = %f" + % (cfg.output_grid.pollon, cfg.output_grid.pollat) + ) + else: + lonname = "rlon" + latname = "rlat" + print( + "Rotated grid: pollon = %f, pollat = %f" + % (cfg.output_grid.pollon, cfg.output_grid.pollat) + ) + elif cfg.model.startswith("icon"): + lonname = None + latname = None + print("ICON grid") # Starts writing out the output file output_file = os.path.join(cfg.output_path, cfg.output_name) with Dataset(output_file, "w") as out: - util.prepare_output_file(cfg.cosmo_grid, cfg.nc_metadata, out) - util.add_country_mask(country_mask, out) + if cfg.model.startswith("cosmo"): + util.prepare_output_file(cfg.output_grid, cfg.nc_metadata, out) + util.add_country_mask(country_mask, out, "cosmo") + elif cfg.model.startswith("icon"): + util.prepare_ICON_output_file(cfg.output_grid, cfg.nc_metadata, out) + util.add_country_mask(country_mask, out, "icon") if cfg.inventory == 'TNO': process_tno(cfg, interpolation, country_mask, out, latname, diff --git a/emiproc/grids.py b/emiproc/grids.py index 86c1cd1..027b47a 100644 --- a/emiproc/grids.py +++ b/emiproc/grids.py @@ -6,7 +6,7 @@ import cartopy.crs as ccrs from netCDF4 import Dataset -from shapely.geometry import Polygon +from shapely.geometry import Polygon, Point class Grid: @@ -743,7 +743,7 @@ def intersected_cells(self, corners): # The inventory cell lies outside the cosmo grid return [] - molly = ccrs.Mollweide() + molly = ccrs.Mollweide(central_longitude = self.pollon) corners = molly.transform_points(self.projection,corners[:,0],corners[:,1]) inv_cell = Polygon(corners) @@ -761,3 +761,206 @@ def intersected_cells(self, corners): intersections.append((i, j, overlap.area / inv_cell.area)) return intersections + + +class ICONGrid(Grid): + """Class to manage an ICON-domain + This grid is defined as an unstuctured triangular grid (1D). + The cells are ordered in a deliberate way and indexed with ascending integer numbers. + The grid file contains variables like midpoint coordinates etc as a fct of the index. + """ + + def __init__(self, dataset_path, name="ICON"): + """Open the netcdf-dataset and read the relevant grid information. + + Parameters + ---------- + dataset_path : str + name : str, optional + """ + self.dataset_path = dataset_path + + with Dataset(dataset_path) as dataset: + self.clon_var = np.array(dataset["clon"][:]) * 180 / np.pi + self.clat_var = np.array(dataset["clat"][:]) * 180 / np.pi + self.cell_areas = np.array(dataset["cell_area"][:]) + self.vlat = np.array(dataset["vlat"][:]) * 180 / np.pi + self.vlon = np.array(dataset["vlon"][:]) * 180 / np.pi + self.vertex_of_cell = np.array(dataset["vertex_of_cell"][:]) + self.cell_of_vertex = np.array(dataset["cells_of_vertex"][:]) + + + self.ncell = len(self.clat_var) + + # Initiate a list of polygons, which is updated whenever the polygon of a cell + # is called for the first time + self.polygons = [None] * self.ncell + + # Consider the ICON-grid as a 1-dimensional grid where ny=1 + self.nx = self.ncell + self.ny = 1 + + super().__init__(name, ccrs.PlateCarree()) + + + def cell_corners(self, n, j): + """Return the corners of the cell with index n. + + Parameters + ---------- + n : int + j : int + + Returns + ------- + tuple(np.array(shape=(3,), dtype=float), + np.array(shape=(3,), dtype=float)) + Arrays containing the lon and lat coordinates of the corners + """ + + return self.vlon[self.vertex_of_cell[:,n]-1], self.vlat[self.vertex_of_cell[:,n]-1] + + + def indices_of_point(self, lon, lat, proj=ccrs.PlateCarree()): + """Return the indices of the grid cell that contains the point (lon, lat) + + Parameters + ---------- + lat : float + The latitude of the point source + lon : float + The longitude of the point source + + Returns + ------- + int + (icon_indn), + the index of the ICON grid cell containing the source. + + Raises + ------ + IndexError + If the point lies outside the grid. + """ + + indn = -1 + + pnt = Point(lon,lat) + + closest_vertex = ((self.vlon - lon)**2 + (self.vlat - lat)**2).argmin() + cell_range = self.cell_of_vertex[:,closest_vertex] - 1 + + for n in cell_range: + + if n == -1: + continue + + if self.polygons[n] is not None: + icon_cell = self.polygons[n] + else: + corners = np.array(list(zip(*self.cell_corners(n,0)))) + icon_cell = Polygon(corners) + self.polygons[n] = icon_cell + + if icon_cell.contains(pnt): + indn = n + break + + if indn == -1: + raise IndexError("Point lies outside the ICON Grid") + + return int(indn), 0 + + def gridcell_areas(self): + """Calculate 2D array of the areas (m^2) of a regular rectangular grid + on earth. + + Returns + ------- + np.array + 2D array containing the areas of the gridcells in m^2 + shape: (ncell) + """ + + return self.cell_areas + + def intersected_cells(self, corners): + """Given a inventory cell, return a list of ICON-cell-indices and + intersection fractions. + + The inventory cell is specified by it's corners. The output is a list + of tuples, specifying the indices and overlap as a fraction of the + inventory cell area. + + Parameters + ---------- + corners : np.array(shape=(4,2), dtype=float) + The corners of the inventory cell in the COSMO coordinate system + + Returns + ------- + list(tuple(int, float)) + A list containing triplets (x,y,r) + - n : index of ICON grid cell + - r : ratio of the area of the intersection compared to the total + area of the inventory cell. + r is in (0,1] (only nonzero intersections are reported) + """ + # Find around which ICON grid index the inventory cell lies. + + cell_xmin = min(k[0] for k in corners) + icon_xmax = max(self.vlon) + + if cell_xmin > icon_xmax: + # The inventory cell lies outside the cosmo grid + return [] + + cell_xmax = max(k[0] for k in corners) + icon_xmin = min(self.vlon) + + if cell_xmax < icon_xmin: + # The inventory cell lies outside the cosmo grid + return [] + + cell_ymin = min(k[1] for k in corners) + icon_ymax = max(self.vlat) + + if cell_ymin > icon_ymax: + # The inventory cell lies outside the cosmo grid + return [] + + cell_ymax = max(k[1] for k in corners) + icon_ymin = min(self.vlat) + + if cell_ymax < icon_ymin: + # The inventory cell lies outside the cosmo grid + return [] + + molly = ccrs.Mollweide() + corners = molly.transform_points(self.projection,corners[:,0],corners[:,1]) + inv_cell = Polygon(corners) + + + intersections = [] + # make sure we iterate only over valid gridpoint indices + for n in np.arange(self.ncell): + icon_cell_xmin = min(self.vlon[self.vertex_of_cell[:,n]-1]) + if cell_xmax < icon_cell_xmin: + continue + icon_cell_xmax = max(self.vlon[self.vertex_of_cell[:,n]-1]) + if cell_xmin > icon_cell_xmax: + continue + icon_cell_ymin = min(self.vlat[self.vertex_of_cell[:,n]-1]) + if cell_ymax < icon_cell_ymin: + continue + icon_cell_ymax = max(self.vlat[self.vertex_of_cell[:,n]-1]) + if cell_ymin > icon_cell_ymax: + continue + corners = np.array(list(zip(*self.cell_corners(n,0)))) + corners = molly.transform_points(self.projection,corners[:,0],corners[:,1]) + icon_cell = Polygon(corners) + if icon_cell.intersects(inv_cell): + overlap = icon_cell.intersection(inv_cell) + intersections.append((n, 0, overlap.area / inv_cell.area)) + + return intersections diff --git a/emiproc/utilities.py b/emiproc/utilities.py index c7dab1d..aa43319 100644 --- a/emiproc/utilities.py +++ b/emiproc/utilities.py @@ -71,6 +71,20 @@ def write_variable(ncfile, variable, var_name, latname, lonname, unit, ncfile[var_name][:] += variable +def write_variable_ICON(ncfile, variable, var_name, unit, overwrite=False): + """ + As function 'write_variable' but for the case of an ICON output grid + """ + if var_name not in ncfile.variables.keys(): + ncfile.createVariable(var_name, float, ("cell")) + + ncfile[var_name].units = unit + ncfile[var_name][:] = variable[0,:] + + elif overwrite: + ncfile[var_name][:] = variable[0,:] + else: + ncfile[var_name][:] += variable[0,:] def read_emi_from_file(path): @@ -110,7 +124,7 @@ def load_cfg(cfg_path): return cfg -def prepare_output_file(cosmo_grid, metadata, dataset): +def prepare_output_file(output_grid, metadata, dataset): """Add lat & lon dimensions and variables to the dataset, handle rotated pole Creates & writes dimensions and variables for longitude and latitude. @@ -119,8 +133,8 @@ def prepare_output_file(cosmo_grid, metadata, dataset): Parameters ---------- - cosmo_grid : COSMOGrid - Contains information about the cosmo grid + output_grid : COSMOGrid + Contains information about the output grid (only called for COSMOGrid) metadata : dict(str : str) Containing global file attributes. Used as argument to netCDF4.Dataset.setncatts. @@ -129,8 +143,8 @@ def prepare_output_file(cosmo_grid, metadata, dataset): """ # Create the dimensions and the rotated pole if ( - cosmo_grid.pollon == 180 or cosmo_grid.pollon == 0 - ) and cosmo_grid.pollat == 90: + output_grid.pollon == 180 or output_grid.pollon == 0 + ) and output_grid.pollat == 90: lonname = "lon" latname = "lat" else: @@ -138,30 +152,63 @@ def prepare_output_file(cosmo_grid, metadata, dataset): latname = "rlat" var_rotpol = dataset.createVariable("rotated_pole", str) var_rotpol.grid_mapping_name = "rotated_latitude_longitude" - var_rotpol.grid_north_pole_latitude = cosmo_grid.pollat - var_rotpol.grid_north_pole_longitude = cosmo_grid.pollon + var_rotpol.grid_north_pole_latitude = output_grid.pollat + var_rotpol.grid_north_pole_longitude = output_grid.pollon var_rotpol.north_pole_grid_longitude = 0.0 - dataset.createDimension(lonname, cosmo_grid.nx) - dataset.createDimension(latname, cosmo_grid.ny) + dataset.createDimension(lonname, output_grid.nx) + dataset.createDimension(latname, output_grid.ny) # Create the variables associated to the dimensions var_lon = dataset.createVariable(lonname, "float32", lonname) var_lon.axis = "X" var_lon.units = "degrees_east" var_lon.standard_name = "longitude" - var_lon[:] = cosmo_grid.lon_range() + var_lon[:] = output_grid.lon_range() var_lat = dataset.createVariable(latname, "float32", latname) var_lat.axis = "Y" var_lat.units = "degrees_north" var_lat.standard_name = "latitude" - var_lat[:] = cosmo_grid.lat_range() + var_lat[:] = output_grid.lat_range() dataset.setncatts(metadata) -def add_country_mask(country_mask, dataset): +def prepare_ICON_output_file(output_grid, metadata, dataset): + """Prepares the output file in the format of the unstructered ICON-grid. + + Copies all dimensions, variables and attributes from the input file + Adds the supplied metadata attributes. + + Parameters + ---------- + output_grid : ICONGrid + Contains information about the output grid (only called for ICONGrid) + metadata : dict(str : str) + Containing global file attributes. Used as argument to + netCDF4.Dataset.setncatts. + dataset : netCDF4.Dataset + Writable (empty) netCDF Dataset + """ + + with nc.Dataset(output_grid.dataset_path) as src: + # copy attributes + for name in src.ncattrs(): + dataset.setncattr(name, src.getncattr(name)) + # copy dimensions + for name, dimension in src.dimensions.items(): + dataset.createDimension( + name, (len(dimension))) + # copy all file data + for name, variable in src.variables.items(): + x = dataset.createVariable(name, variable.datatype, variable.dimensions) + dataset.variables[name][:] = src.variables[name][:] + + dataset.setncatts(metadata) + + +def add_country_mask(country_mask, dataset, model): """Create and write the country mask to the dataset. Parameters @@ -171,19 +218,26 @@ def add_country_mask(country_mask, dataset): dataset : netCDF4.Dataset Writable netCDF Dataset """ - if "rotated_pole" in dataset.variables: - var = dataset.createVariable("country_ids", "short", ("rlat", "rlon")) - var.grid_mapping = "rotated_pole" - else: - var = dataset.createVariable("country_ids", "short", ("lat", "lon")) + if model == "cosmo": + if "rotated_pole" in dataset.variables: + var = dataset.createVariable("country_ids", "short", ("rlat", "rlon")) + var.grid_mapping = "rotated_pole" + else: + var = dataset.createVariable("country_ids", "short", ("lat", "lon")) + elif model == "icon": + var = dataset.createVariable("country_ids", "short", ("cell")) var.long_name = "EMEP_country_code" # Transpose the country mask to conform with the storage of netcdf # python: (lon, lat), FORTRAN: (lat, lon) - var[:] = country_mask.T + + if model == "cosmo": + var[:] = country_mask.T + elif model == "icon": + var[:] = country_mask -def get_country_mask(output_path, suffix, cosmo_grid, resolution, nprocs): +def get_country_mask(output_path, suffix, output_grid, resolution, nprocs): """Returns the country-mask, either loaded from disk or computed. If there already exists a file at @@ -195,8 +249,8 @@ def get_country_mask(output_path, suffix, cosmo_grid, resolution, nprocs): output_path : str Path to the directory where the country-mask is stored suffix : str - cosmo_grid : grids.COSMOGrid - Contains all necessary information about the cosmo grid + output_grid : grids.COSMOGrid or grids.ICONGrid + Contains all necessary information about the output grid resolution : str The resolution for the used shapefile, used as argument for cartopy.io.shapereader.natural_earth() @@ -205,7 +259,7 @@ def get_country_mask(output_path, suffix, cosmo_grid, resolution, nprocs): Returns ------- - np.array(shape(cosmo_grid.nx, cosmo_grid.ny), dtype=int) + np.array(shape(output_grid.nx, output_grid.ny), dtype=int) """ cmask_path = os.path.join( output_path, f"country_mask_{resolution}_{suffix}.nc" @@ -229,10 +283,15 @@ def get_country_mask(output_path, suffix, cosmo_grid, resolution, nprocs): "AFFILIATION": "Empa Duebendorf, Switzerland", "DATE CREATED": time.ctime(time.time()), } - country_mask = compute_country_mask(cosmo_grid, resolution, nprocs) + country_mask = compute_country_mask(output_grid, resolution, nprocs) + print(np.nanmax(country_mask)) with nc.Dataset(cmask_path, "w") as dataset: - prepare_output_file(cosmo_grid, nc_metadata, dataset) - add_country_mask(country_mask, dataset) + if output_grid.__class__.__name__ == "COSMOGrid": + prepare_output_file(output_grid, nc_metadata, dataset) + add_country_mask(country_mask, dataset, "cosmo") + elif output_grid.__class__.__name__ == "ICONGrid": + prepare_ICON_output_file(output_grid, nc_metadata, dataset) + add_country_mask(country_mask, dataset, "icon") else: # Transpose country mask when reading in # netCDF/Fortran: (lat, lon), python: (lon, lat) @@ -241,7 +300,7 @@ def get_country_mask(output_path, suffix, cosmo_grid, resolution, nprocs): return country_mask -def compute_country_mask(cosmo_grid, resolution, nprocs): +def compute_country_mask(output_grid, resolution, nprocs): """Determine the country-code for each gridcell and return the grid. Each gridcell gets assigned to code of the country with the most @@ -255,8 +314,8 @@ def compute_country_mask(cosmo_grid, resolution, nprocs): Parameters ---------- - cosmo_grid : grids.COSMOGrid - Contains all necessary information about the cosmo grid + output_grid : grids.COSMOGrid or grids.ICONGrid + Contains all necessary information about the output grid resolution : str The resolution for the used shapefile, used as argument for cartopy.io.shapereader.natural_earth() @@ -265,7 +324,7 @@ def compute_country_mask(cosmo_grid, resolution, nprocs): Returns ------- - np.array(shape(cosmo_grid.nx, cosmo_grid.ny), dtype=int) + np.array(shape(output_grid.nx, output_grid.ny), dtype=int) """ print( f"Computing the country mask with {resolution} resolution.\n" @@ -282,18 +341,18 @@ def compute_country_mask(cosmo_grid, resolution, nprocs): projections = { # Projection of the shapefile: WGS84, which the PlateCarree defaults to "shapefile": ccrs.PlateCarree(), - "cosmo": cosmo_grid.get_projection(), + "cosmo": output_grid.get_projection(), } # arguments to assign_country_code() arguments = [ ( - [(i, j) for j in range(cosmo_grid.ny)], # indices of cells - cosmo_grid, + [(i, j) for j in range(output_grid.ny)], # indices of cells + output_grid, projections, countries, ) - for i in range(cosmo_grid.nx) + for i in range(output_grid.nx) ] if nprocs == 1: @@ -313,7 +372,7 @@ def compute_country_mask(cosmo_grid, resolution, nprocs): return country_mask -def assign_country_code(indices, cosmo_grid, projections, countries): +def assign_country_code(indices, output_grid, projections, countries): """Assign the country codes on the gridcells indicated by indices. Each gridcell gets assigned to code of the country with the most @@ -330,8 +389,8 @@ def assign_country_code(indices, cosmo_grid, projections, countries): indices : list(tuple(int, int)) A list of indices indicating for which gridcells to compute the country code. - cosmo_grid : grids.COSMOGrid - Contains all necessary information about the cosmo grid + output_grid : grids.COSMOGrid or grids.ICONGrid + Contains all necessary information about the output grid projections : dict(str, cartopy.crs.Projection) Dict containing two elements: 'shapefile': Projection used in the coordinates of the countries @@ -347,9 +406,10 @@ def assign_country_code(indices, cosmo_grid, projections, countries): # Have to recreate the COSMO projection since cartopy projections # don't work with deepcopy (which is used by multiprocessing to send # arguments of functions to other processes) - projections["cosmo"] = ccrs.RotatedPole( - pole_longitude=cosmo_grid.pollon, pole_latitude=cosmo_grid.pollat - ) + if output_grid.__class__.__name__ == "COSMOGrid": + projections["cosmo"] = ccrs.RotatedPole( + pole_longitude=output_grid.pollon, pole_latitude=output_grid.pollat + ) # Name of the country attribute holding the ISO3 country abbreviation iso3 = "ADM0_A3" @@ -362,11 +422,18 @@ def assign_country_code(indices, cosmo_grid, projections, countries): country_mask = np.empty(len(indices), dtype=int) for k, (i, j) in enumerate(indices): - cosmo_cell_x, cosmo_cell_y = cosmo_grid.cell_corners(i, j) - projected_corners = projections["shapefile"].transform_points( - projections["cosmo"], cosmo_cell_x, cosmo_cell_y - ) + output_cell_x, output_cell_y = output_grid.cell_corners(i, j) + + if output_grid.__class__.__name__ == "COSMOGrid": + projected_corners = projections["shapefile"].transform_points( + projections["cosmo"], output_cell_x, output_cell_y + ) + elif output_grid.__class__.__name__ == "ICONGrid": + projected_corners = np.array([ + output_cell_x, output_cell_y, np.zeros(output_cell_x.shape) + ]).T + projected_cell = Polygon(projected_corners) intersected_countries = [ @@ -460,7 +527,7 @@ def intersection_possible(points, c_minx, c_miny, c_maxx, c_maxy): return c_miny < p_maxy and c_maxy > p_miny -def compute_map_from_inventory_to_cosmo(cosmo_grid, inv_grid, nprocs): +def compute_map_from_inventory_to_cosmo(output_grid, inv_grid, nprocs): """Compute the mapping from inventory to cosmo grid. Loop over all inventory cells and determine which cosmo cells they overlap @@ -473,8 +540,8 @@ def compute_map_from_inventory_to_cosmo(cosmo_grid, inv_grid, nprocs): Parameters ---------- - cosmo_grid : grids.COSMOGrid - Contains all necessary information about the cosmo grid + output_grid : grids.COSMOGrid or grids.ICONGrid + Contains all necessary information about the output grid inv_grid : grids.InventoryGrid Contains all necessary information about the inventory grid nprocs : int @@ -497,7 +564,7 @@ def compute_map_from_inventory_to_cosmo(cosmo_grid, inv_grid, nprocs): inv_projection = inv_grid.get_projection() # Projection used to convert to the cosmo grid - cosmo_projection = cosmo_grid.get_projection() + output_projection = output_grid.get_projection() progress = ProgressIndicator(lon_size) @@ -509,12 +576,12 @@ def compute_map_from_inventory_to_cosmo(cosmo_grid, inv_grid, nprocs): inv_cell_corners_x, inv_cell_corners_y = inv_grid.cell_corners( i, j ) - cell_in_cosmo_projection = cosmo_projection.transform_points( + cell_in_output_projection = output_projection.transform_points( inv_projection, inv_cell_corners_x, inv_cell_corners_y ) - cells.append(cell_in_cosmo_projection) + cells.append(cell_in_output_projection) - mapping[i, : ] = [cosmo_grid.intersected_cells(c) for c in cells] + mapping[i, : ] = [output_grid.intersected_cells(c) for c in cells] else: with Pool(nprocs) as pool: for i in range(lon_size): @@ -524,12 +591,12 @@ def compute_map_from_inventory_to_cosmo(cosmo_grid, inv_grid, nprocs): inv_cell_corners_x, inv_cell_corners_y = inv_grid.cell_corners( i, j ) - cell_in_cosmo_projection = cosmo_projection.transform_points( + cell_in_output_projection = output_projection.transform_points( inv_projection, inv_cell_corners_x, inv_cell_corners_y ) - cells.append(cell_in_cosmo_projection) + cells.append(cell_in_output_projection) - mapping[i, :] = pool.map(cosmo_grid.intersected_cells, cells) + mapping[i, :] = pool.map(output_grid.intersected_cells, cells) end = time.time() @@ -538,7 +605,7 @@ def compute_map_from_inventory_to_cosmo(cosmo_grid, inv_grid, nprocs): return mapping -def get_gridmapping(output_path, suffix, cosmo_grid, inv_grid, nprocs): +def get_gridmapping(output_path, suffix, output_grid, inv_grid, nprocs): """Returns the interpolation between the inventory and COSMO grid. If there already exists a file at output_path/mapping_{suffix}.npy ask @@ -549,8 +616,8 @@ def get_gridmapping(output_path, suffix, cosmo_grid, inv_grid, nprocs): output_path : str Path to the directory where the country-mask is stored suffix : str - cosmo_grid : grids.COSMOGrid - Contains all necessary information about the cosmo grid + output_grid : grids.COSMOGrid or grids.ICONGrid + Contains all necessary information about the output grid inv_grid : grids.Grid Contains all necessary information about the inventory grid nprocs : int @@ -575,7 +642,7 @@ def get_gridmapping(output_path, suffix, cosmo_grid, inv_grid, nprocs): if compute_map: mapping = compute_map_from_inventory_to_cosmo( - cosmo_grid, inv_grid, nprocs + output_grid, inv_grid, nprocs ) np.save(mapping_path, mapping)