Feature/pcastell/improve update

CHANGELOG.md

@@ -9,8 +9,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 ### Added
+- add optional explicit bin ordering to aop.py yaml. this allows for calculating AOP's for a single bin, or subset of bins
+
 ### Changed
-### Remoed
+- use xesmf regridder to station sampling. this is more efficient that using xarray native interp
+
+### Removed
 ### Deprecated
 
 # [1.7.0] - yyyy-mm-dd

src/pyobs/aop.py

@@ -42,6 +42,12 @@
     - DU003
     - DU004
     - DU005
+  bin:
+    - 1
+    - 2
+    - 3
+    - 4
+    - 5
   shapefactor: 1.4
   rhod:
     - 2500
@@ -61,6 +67,12 @@
     - SS003
     - SS004
     - SS005
+  bin:
+    - 1
+    - 2
+    - 3
+    - 4
+    - 5
   shapefactor: 1
   rhod:
     - 2200
@@ -77,6 +89,9 @@
   tracers:
     - OCPHOBIC
     - OCPHILIC
+  bin:
+    - 1
+    - 2
   shapefactor: 1
   rhod:
     - 1800
@@ -89,6 +104,9 @@
   tracers:
     - BCPHOBIC
     - BCPHILIC
+  bin:
+    - 1
+    - 2
   shapefactor: 1
   rhod:
     - 1800
@@ -101,6 +119,9 @@
   tracers:
     - BRCPHOBIC
     - BRCPHILIC
+  bin:
+    - 1
+    - 2
   shapefactor: 1
   rhod:
     - 1800
@@ -112,6 +133,8 @@
   bandFile: ExtData/chemistry/AerosolOptics/v1.0.0/x/opticsBands_SU.v1_3.RRTMG.nc4
   tracers:
     - SO4
+  bin:
+    - 1
   shapefactor: 1
   rhod:
     - 1700
@@ -124,6 +147,10 @@
     - NO3AN1
     - NO3AN2
     - NO3AN3
+  bin:
+    - 1
+    - 2
+    - 3
   shapefactor: 1
   rhod:
     - 1725
@@ -322,8 +349,11 @@ def getAOPrt(self,Species=None,wavelength=None,vector=False,
             Tracers = self.mieTable[s]['tracers']
             mie = self.mieTable[s]['mie']
 
-            bin = 1
-            for q in Tracers:
+            bins = self.mieTable[s].get('bin',[])
+            if not bins:
+                bins = list(range(1,len(Tracers)+1))
+
+            for q,bin in zip(Tracers,bins):
 
                 if self.verbose:
                     print('   -',q)
@@ -357,8 +387,6 @@ def getAOPrt(self,Species=None,wavelength=None,vector=False,
                     g   += sca_ * g_
 
 
-                bin += 1
-
         # Final normalization of SSA and g
         # protect against divide by zero
         # this can happen if you ask for the AOP of an individual species
@@ -497,8 +525,11 @@ def getAOPext(self,Species=None,wavelength=None,fixrh=None,doaback=True):
             Tracers = self.mieTable[s]['tracers']
             mie = self.mieTable[s]['mie']
 
-            bin = 1
-            for q in Tracers:
+            bins = self.mieTable[s].get('bin',[])
+            if not bins:
+                bins = list(range(1,len(Tracers)+1))
+
+            for q,bin in zip(Tracers,bins):
 
                 if self.verbose:
                     print('   -',q)
@@ -521,7 +552,6 @@ def getAOPext(self,Species=None,wavelength=None,fixrh=None,doaback=True):
                 depol1 += (pback11_-pback22_) * sca_
                 depol2 += (pback11_+pback22_) * sca_
 
-                bin += 1
 
         if doaback:
             # Compute Molecular Scattering and Total Attenuated Backscatter Coefficient
@@ -748,8 +778,14 @@ def getPM(self,Species=None,pmsize=None,fixrh=None,aerodynamic=False,vacuum_aero
             Tracers = self.mieTable[s]['tracers']
             mie = self.mieTable[s]['mie']
 
-            bin = 1
-            for q in Tracers:
+            bins = self.mieTable[s].get('bin',[])
+            if not bins:
+                bins = list(range(1,len(Tracers)+1))
+
+            # Dry aerosol density in kg m-3
+            # rhod is not in all of the standard optics files, and is for now read from the yaml config
+            rhod = self.mieTable[s]['rhod']
+            for q,bin,rhod_ in zip(Tracers,bins,rhod):
 
                 if self.verbose:
                     print('   -',q)
@@ -758,10 +794,7 @@ def getPM(self,Species=None,pmsize=None,fixrh=None,aerodynamic=False,vacuum_aero
                 # Aerosol mass concentration in kg/m3                
                 q_conc = (a['AIRDENS'] * a[q]).values
 
-                # Dry aerosol density in kg m-3
-                # rhod is not in all of the standard optics files, and is for now read from the yaml config 
                 # rhod_ = mie.getAOP('rhod',  bin, rh, wavelength=wavelength).values
-                rhod_ = self.mieTable[s]['rhod'][bin-1] 
 
                 # Lower and upper bound of the bin's radius converted from meters to microns
                 rLow_ = mie.getBinInfo('rLow', bin)*1000000 
@@ -816,8 +849,6 @@ def getPM(self,Species=None,pmsize=None,fixrh=None,aerodynamic=False,vacuum_aero
                 waterfactor = 1 - (1./growthfactor)
                 wm += pm_ * waterfactor
 
-                bin += 1
-
 
         # Get the fraction of the total PM that is water
         # ------------------------------------------------

src/pyobs/sampler.py

@@ -8,13 +8,16 @@
 
 import numpy  as np
 import xarray as xr
+import xesmf as xe
 import pandas as pd
 
 from datetime import datetime, timedelta
 from dateutil.parser import parse as isoparser
 from glob import glob
 
 from . import xrctl as xc
+import fsspec
+import dask.distributed
 
 os.environ['HDF5_USE_FILE_LOCKING']='FALSE'
 
@@ -59,7 +62,12 @@ def __init__(self, stations, lons, lats,
         # -------------------------------------------------
         if isinstance(dataset,xr.Dataset):
             self.ds = dataset # we are good to go...
-
+        elif isinstance(dataset,dask.distributed.Future):
+            self.ds = xr.open_dataset(dataset, engine="zarr",chunks=chunks, backend_kwargs={"consolidated": False})
+        #  if dataset is a parquet referece file path
+        elif dataset[-4:] == 'parq':
+            fs = fsspec.filesystem("reference", fo=dataset, remote_protocol='file', lazy=True)
+            self.ds = xr.open_dataset(fs.get_mapper(), engine="zarr",chunks=chunks, backend_kwargs={"consolidated": False})        
         # If dataset is a list of files...
         # OR GrADS-style ctl
         # OR a glob type of template
@@ -79,11 +87,12 @@ def __init__(self, stations, lons, lats,
         if times is not None:
             self.times = xr.DataArray(times,dims='time')
         else:
-            self.times = self.ds.time
+            self.times = times
 
-        # TO DO: when using xESMF for regridding, pre-compute transforms here
+        # Use xESMF for regridding, pre-compute transforms here
         # -------------------------------------------------------------------
-
+        ds_loc = xr.Dataset({"lon": self.lons, "lat": self.lats})
+        self.regridder = xe.Regridder(self.ds, ds_loc, "bilinear", locstream_out=True)
 
     #--
     def sample(self,Variables=None,method='linear'):
@@ -101,7 +110,11 @@ def sample(self,Variables=None,method='linear'):
 
         for vn in Variables:
             if self.verb: print('[ ] sampling ',vn)
-            sampled[vn] = self.ds[vn].interp(time=self.times,lon=self.lons,lat=self.lats,method=method)
+            sampled[vn] = self.regridder(self.ds[vn])
+#            sampled[vn] = self.ds[vn].interp(time=self.times,lon=self.lons,lat=self.lats,method=method)
+
+            if self.times is not None:
+                sampled[vn] = sampled[vn].interp(time=self.times,method=method)
 
         return xr.Dataset(sampled).assign_coords({'station': self.stations})
 
@@ -136,7 +149,10 @@ def __init__(self, times, lons, lats, dataset, parallel=True,chunks='auto',verbo
         # -------------------------------------------------
         if isinstance(dataset,xr.Dataset):
             self.ds = dataset # we are good to go...
-
+        #  if dataset is a parquet referece file
+        elif dataset[-4:] == 'parq':
+            fs = fsspec.filesystem("reference", fo=dataset, remote_protocol='file', lazy=True)
+            self.ds = xr.open_dataset(fs.get_mapper(), engine="zarr", chunks=chunks, backend_kwargs={"consolidated": False})
         # If dataset is a list of files...
         # OR GrADS-style ctl 
         # OR a glob type of template