HiTRes Inversion Options

CHANGELOG.md

@@ -34,6 +34,7 @@ All calls to species_info.json and site_info.json are now made through the openg
 - hbmcmc code now has a function (and relevant additions to output) that allows the inversion to be rerun (i.e. reproduced)
   using only the output as inputs and ACRG repository.
 - Country mask code has now been updated to allow ocean territories (Economic Exclusion Zone, EEZ) to be included.
+- hbmcmc code now has options to use HiTRes processes for the set up
 - Add function for embedding regional fields into larger fields, in acrg.name.emissions_helper_func
 
 ### Changed

acrg/BC/timevarying_BC.py

@@ -32,15 +32,23 @@
 -------
 
 '''
-if 'os'    not in dir(): import os
-if 'np'    not in dir(): import numpy as np
-if 'xr'    not in dir(): import xarray as xr
-if 'paths' not in dir(): from acrg_config.paths import paths
+import os
+import json
+import numpy as np
+import xarray as xr
+from scipy import interpolate
+from collections import OrderedDict
+
+from acrg.name.name import open_ds
+from acrg.config.paths import Paths
+from acrg.convert import concentration
+from acrg.countrymask import domain_volume
+
 
 class BoundaryConditions:
     def __init__(self, vmr_var, gph_height_var, dataset=None, filename=None,
-                 file_dir=None, time_coord='time', height_coord='height', species=None,
-                 domain='EUROPE', start_date=None):
+                 file_dir=None, time_coord='time', species=None,
+                 domain='EUROPE', start_date=None, adjust=None):
         '''
         vmr_var : str
             The VMR variable name in the nesw dataset
@@ -73,13 +81,12 @@ def __init__(self, vmr_var, gph_height_var, dataset=None, filename=None,
             if filename is None:
                 print('Please provide either a dataset or filename containing vmr')
             elif os.path.isfile(filename):
-                from acrg_name.name import open_ds
                 dataset = open_ds(filename)
             else:
                 print(f'Cannot find file: {filename}')
 
-        self.data_path    = paths.data
-        self.acrg_path    = paths.acrg
+        self.data_path    = Paths.data
+        self.acrg_path    = Paths.acrg
 
         self.vmr_var      = vmr_var
         self.height_var   = gph_height_var
@@ -91,7 +98,12 @@ def __init__(self, vmr_var, gph_height_var, dataset=None, filename=None,
         # rename latitude and longitude to short version
         if 'latitude'  in dataset: dataset.rename({'latitude'  : 'lat'})
         if 'longitude' in dataset: dataset.rename({'longitude' : 'lon'})
-        self.dataset      = dataset
+
+        if adjust is not None:
+            print(f'Adjusting vmr by {adjust}')
+            dataset[vmr_var] = dataset[vmr_var] + adjust
+
+        self.dataset = dataset
 
     def make_bc_file(self, fp_directory=None, fp_height_coord='height', reverse=None,
                      convert_units=True, datasource=None, out_path=None, glob_attrs={},
@@ -145,7 +157,7 @@ def make_bc_file(self, fp_directory=None, fp_height_coord='height', reverse=None
         if verbose: print('\nSaving to netcdf\n----------------')
         # save to a netcdf file with a standardised filename
         self.to_netcdf(datasource=datasource, out_path=out_path, glob_attrs=glob_attrs,
-                       copy_glob_attrs=copy_glob_attrs, verbose=verbose)
+                       copy_glob_attrs=copy_glob_attrs, verbose=verbose, overwrite=False)
 
 
     def get_unit_conversion(self, verbose=True):
@@ -190,8 +202,6 @@ def cut_to_edge(self, fp_directory=None, fp_height_coord='height', convert_units
                 for north, south, east, and west
                 Can be accessed using BoundaryConditions_obj.edges
         '''
-        from acrg_countrymask import domain_volume
-
         fp_directory = os.path.join(self.data_path, 'LPDM', 'fp_NAME') if fp_directory is None else fp_directory
         self.fp_lat, self.fp_lon, self.fp_height = domain_volume(self.domain, fp_directory=fp_directory)
 
@@ -205,16 +215,16 @@ def cut_to_edge(self, fp_directory=None, fp_height_coord='height', convert_units
             lat_n -= 1
 
         lat_s = (np.abs(self.dataset.coords['lat'].values - min(self.fp_lat))).argmin()
-        if self.dataset.coords['lat'].values[lat_s] > min(self.fp_lat) and lat_s != (len(dataset.coords['lat'].values)-1):
+        if self.dataset.coords['lat'].values[lat_s] > min(self.fp_lat) and lat_s != (len(self.dataset.coords['lat'].values)-1):
             lat_s += 1
 
         lon_e = (np.abs(self.dataset.coords['lon'].values - max(self.fp_lon))).argmin()
-        if self.dataset.coords['lon'].values[lon_e] < max(self.fp_lon) and lon_e != (len(dataset.coords['lon'].values)-1):
-            lat_s += 1
+        if self.dataset.coords['lon'].values[lon_e] < max(self.fp_lon) and lon_e != (len(self.dataset.coords['lon'].values)-1):
+            lon_e += 1
 
         lon_w = (np.abs(self.dataset.coords['lon'].values - min(self.fp_lon))).argmin()
         if self.dataset.coords['lon'].values[lon_w] > min(self.fp_lon) and lon_w != 0:
-            lat_w -= 1
+            lon_w -= 1
 
         # Cut to these
         north = self.dataset.sel(lat  = self.dataset.coords['lat'][lat_n],
@@ -288,8 +298,6 @@ def interp_height_single(self, direction, height_coord='height', reverse=None, n
                 with a 3D VMR array interpolated to the NAME heights
                 Can be accessed using BoundaryConditions_obj.edges
         '''
-        if 'interpolate' not in dir(): from scipy import interpolate
-
         # get the axis along which to interpolate
         # either 'longitude' (or 'lon', N or S) or 'latitude' (or 'lat', E or W)
         latorlon = 'lon' if 'lon' in self.edges[direction] else 'lat'
@@ -358,8 +366,6 @@ def interp_latlon_single(self, direction=None, height_coord='height', reverse=No
                 with a 3D VMR array interpolated to the NAME latitudes and longitudes
                 Can be accessed using BoundaryConditions_obj.edges
         '''
-        if 'interpolate' not in dir(): from scipy import interpolate
-
         # get the axis along which to interpolate
         # either 'longitude' (or 'lon', N or S) or 'latitude' (or 'lat', E or W)
         latorlon = 'lon'       if 'lon'       in self.edges[direction] else \
@@ -418,7 +424,7 @@ def bc_filename(self, from_climatology=False, verbose=True):
         if verbose: print(f'Output filename : {self.out_filename}')
 
     def to_netcdf(self, datasource=None, out_path=None, glob_attrs={}, from_climatology=False,
-                  copy_glob_attrs=False, verbose=True):
+                  copy_glob_attrs=False, verbose=True, overwrite=False):
         '''
         Args
             datasource (str, optional)
@@ -449,7 +455,7 @@ def to_netcdf(self, datasource=None, out_path=None, glob_attrs={}, from_climatol
 
             Includes DataArrays for the n, s, e, and w boundaries
         '''        
-        edges    = {dd : self.edges[dd].rename({self.vmr_var: f'vmr_{dd}'})
+        edges    = {dd : self.edges[dd].rename({self.vmr_var: f'vmr_{dd[0]}'})
                     for dd in ['north', 'south', 'east', 'west']}
 
         # merge the interpolated n, s, e, & w Datasets
@@ -479,8 +485,22 @@ def to_netcdf(self, datasource=None, out_path=None, glob_attrs={}, from_climatol
                        os.path.join(self.data_path, 'LPDM', 'bc', self.domain)
 
         self.bc_filename(from_climatology=from_climatology, verbose=verbose)
-        if verbose: print(f'Saving boundary conditions to : {os.path.join(out_path, self.out_filename)}')
-        BC_edges.to_netcdf(path = os.path.join(out_path, self.out_filename), mode = 'w')
+        if overwrite and verbose:
+            save = True
+            if os.path.isfile(os.path.join(out_path, self.out_filename)):
+                print(f'Boundary condition file {os.path.join(out_path, self.out_filename)} already exists and is being overwritten.')
+        elif os.path.isfile(os.path.join(out_path, self.out_filename)) and not overwrite:
+            print(f'Boundary condition file {os.path.join(out_path, self.out_filename)} already exists.')
+            answer = input("You are about to overwrite an existing file, do you want to continue? Y/N ")
+            if answer.upper() == 'N':
+                save = False
+            elif answer.upper() == 'Y':
+                save = True
+        else:
+            save = True
+        if save == True:
+            if verbose: print(f'Saving boundary conditions to : {os.path.join(out_path, self.out_filename)}')
+            BC_edges.to_netcdf(path = os.path.join(out_path, self.out_filename), mode = 'w')
 
 
 

acrg/convert.py

@@ -8,6 +8,8 @@
 
 acrg_path = Paths.acrg
 
+import numpy as np
+
 def molar_mass(species):
     '''
     This function extracts the molar mass of a species from the acrg_species_info.json file.
@@ -69,5 +71,4 @@ def convert_lons_0360(lons):
 
     div = lons // 360
 
-    return lons - div*360
-
+    return lons - div*360

acrg/hbmcmc/hbmcmc.py

@@ -36,6 +36,8 @@
 def fixedbasisMCMC(species, sites, domain, meas_period, start_date, 
                    end_date, outputpath, outputname,
                    met_model = None,
+                   species_footprint = None,
+                   HiTRes = False,
                    xprior={"pdf":"lognormal", "mu":1, "sd":1},
                    bcprior={"pdf":"lognormal", "mu":0.004, "sd":0.02},
                    sigprior={"pdf":"uniform", "lower":0.5, "upper":3},
@@ -71,6 +73,11 @@ def fixedbasisMCMC(species, sites, domain, meas_period, start_date,
             End time of inversion "YYYY-mm-dd"
         outputname (str):
             Unique identifier for output/run name.
+        species_footprint (str, optional):
+            species of the footprint to be imported, if different to the species
+            of interest (e.g. import co2 footprints for HiTRes studies)
+        HiTRes (bool)
+            True if using HiTRes footprints
         outputpath (str):
             Path to where output should be saved.
         xprior (dict):
@@ -174,15 +181,23 @@ def fixedbasisMCMC(species, sites, domain, meas_period, start_date,
     TO DO:
         Add a wishlist...
     """    
+    keep_missing = True if HiTRes else False
+    if verbose and species_footprint is not None:
+        print(f'species_footprint: {species_footprint}')
+
     data = getobs.get_obs(sites, species, start_date = start_date, end_date = end_date, 
                          average = meas_period, data_directory=obs_directory,
-                          keep_missing=False,inlet=inlet, instrument=instrument, max_level=max_level)
-    fp_all = name.footprints_data_merge(data, domain=domain, met_model = met_model, calc_bc=True, 
-                                        height=fpheight, 
+                          keep_missing=keep_missing,inlet=inlet, instrument=instrument,
+                          max_level=max_level)
+    fp_all = name.footprints_data_merge(data, domain=domain, met_model = met_model, calc_bc=True,
+                                        HiTRes = HiTRes,
+                                        height = fpheight,
+                                        calc_timeseries = False,
                                         fp_directory = fp_directory,
                                         bc_directory = bc_directory,
                                         flux_directory = flux_directory,
-                                        emissions_name=emissions_name)
+                                        emissions_name = emissions_name,
+                                        species_footprint = species_footprint)
 
     for site in sites:
         for j in range(len(data[site])):
@@ -221,10 +236,15 @@ def fixedbasisMCMC(species, sites, domain, meas_period, start_date,
             basis_directory = tempdir
     else:
         basis_directory = basis_directory
-
-    fp_data = name.fp_sensitivity(fp_all, domain=domain, basis_case=fp_basis_case,basis_directory=basis_directory)
-    fp_data = name.bc_sensitivity(fp_data, domain=domain,basis_case=bc_basis_case)
 
+    fp_data = name.fp_sensitivity(fp_all, domain=domain, basis_case=fp_basis_case, basis_directory=basis_directory,
+                                  calc_timeseries = True)
+    fp_data = name.bc_sensitivity(fp_data, domain=domain, basis_case=bc_basis_case)
+
+    if HiTRes:
+        for site in sites:
+            fp_data[site] = fp_data[site].dropna(dim='time')
+
     #apply named filters to the data
     fp_data = name.filtering(fp_data, filters)
 

acrg/hbmcmc/inversion_pymc3.py

@@ -331,10 +331,6 @@ def inferpymc3_postprocessouts(xouts,bcouts, sigouts, convergence,
                 with {source_name: emissions_file_identifier} (e.g. {'anth':'co2-ff-mth'}). This way
                 multiple sources can be read in simultaneously if they are added as separate entries to
                 the emissions_name dictionary.
-                If using HiTRes footprints, both the high and low frequency emissions files must be specified
-                in a second dictionary like so: {'anth': {'high_freq':'co2-ff-2hr', 'low_freq':'co2-ff-mth'}}.
-                It is not a problem to have a mixture of sources, with some that use HiTRes footprints and some
-                that don't.
             basis_directory (str, optional):
                 Directory containing basis function file
             country_file (str, optional):
@@ -407,10 +403,8 @@ def inferpymc3_postprocessouts(xouts,bcouts, sigouts, convergence,
             # it varies over the inversion period
             emissions_flux = np.expand_dims(rerun_file.fluxapriori.values,2)
         else:
-            if emissions_name == None:
-                emds = name.flux(domain, species, start = start_date, end = end_date, flux_directory=flux_directory)
-            else:
-                emds = name.flux(domain, list(emissions_name.values())[0], start = start_date, end = end_date, flux_directory=flux_directory)
+            spec = species if emissions_name == None else list(emissions_name.values())[0]
+            emds = name.flux(domain=domain, species=spec, start=start_date, end=end_date, flux_directory=flux_directory)
             emissions_flux = emds.flux.values
         flux = scalemap*emissions_flux[:,:,0]