Merge pull request #98 from brian-rose/fix_interpolation

Fix interpolation

climlab/__init__.py

@@ -7,7 +7,7 @@
 has been documented for a comprehensive understanding and traceability.
 '''
 
-__version__ = '0.7.2.dev1'
+__version__ = '0.7.2.dev2'
 
 # this should ensure that we can still import constants.py as climlab.constants
 from .utils import constants, thermo, legendre

climlab/radiation/radiation.py

@@ -143,16 +143,26 @@ def default_absorbers(Tatm,
         ozonedata = xr.open_dataset(ozonefilepath)
         ##  zonal and time average
         ozone_zon = ozonedata.OZONE.mean(dim=('time','lon')).transpose('lat','lev')
-        weight = np.cos(np.deg2rad(ozonedata.lat))
-        ozone_global = (ozone_zon * weight).mean(dim='lat') / weight.mean(dim='lat')
+        if ('lat' in xTatm.dims):
+            O3source = ozone_zon
+        else:
+            weight = np.cos(np.deg2rad(ozonedata.lat))
+            ozone_global = (ozone_zon * weight).mean(dim='lat') / weight.mean(dim='lat')
+            O3source = ozone_global
         try:
-            if ('lat' in xTatm.dims):
-                O3 = ozone_zon.interp_like(xTatm)
-            else:
-                O3 = ozone_global.interp_like(xTatm)
+            O3 = O3source.interp_like(xTatm)
+            # There will be NaNs for gridpoints outside the ozone file domain
+            assert not np.any(np.isnan(O3))
         except:
-            print('Interpolation of ozone data failed.')
-            print('Setting O3 to zero instead.')
+            warnings.warn('Some grid points are beyond the bounds of the ozone file. Ozone values will be extrapolated.')
+            try:
+                # passing fill_value=None to the underlying scipy interpolator
+                # will result in extrapolation instead of NaNs
+                O3 = O3source.interp_like(xTatm, kwargs={'fill_value':None})
+                assert not np.any(np.isnan(O3))
+            except:
+                warnings.warn('Interpolation of ozone data failed. Setting O3 to zero instead.')
+                O3 = 0. * xTatm
     absorber_vmr['O3'] = O3.values
     return absorber_vmr
 

climlab/tests/test_rrtm.py

@@ -179,3 +179,13 @@ def test_fixed_insolation():
     ins_array = insolation.values
     rad = climlab.radiation.RRTMG(name='Radiation', state=state, insolation=ins_array)
     rad.step_forward()
+
+@pytest.mark.compiled
+@pytest.mark.fast
+def test_large_grid():
+    num_lev = 50; num_lat=90
+    state = climlab.column_state(num_lev=num_lev, num_lat=num_lat, water_depth=10.)
+    rad1 = climlab.radiation.CAM3(state=state)
+    rad1.step_forward()
+    rad2 = climlab.radiation.RRTMG(state=state)
+    rad2.step_forward()

conda-recipe/meta.yaml

@@ -1,4 +1,4 @@
-{% set version = "0.7.2.dev1" %}
+{% set version = "0.7.2.dev2" %}
 
 package:
   name: climlab

setup.py

@@ -1,7 +1,7 @@
 import os, sys
 import textwrap
 
-VERSION = '0.7.2.dev1'
+VERSION = '0.7.2.dev2'
 
 # BEFORE importing setuptools, remove MANIFEST. Otherwise it may not be
 # properly updated when the contents of directories change (true for distutils,