Update Develop

CITATION.cff

@@ -6,8 +6,12 @@ authors:
     orcid: https://orcid.org/0000-0002-7127-6907
   - family-names: Levis
     given-names: Aviad
+    orcid: https://orcid.org/0000-0001-7307-632X
+  - family-names: Aides
+    given-names: Amit
   - family-names: Forster
     given-names: Linda
+    orcid: https://orcid.org/0000-0002-9738-9571
   - family-names: Holodovsky
     given-names: Vadim
 title: "Atmospheric Tomography with 3D Radiative Transfer"

README.md

@@ -1,17 +1,32 @@
 # Atmospheric Tomography with 3D Radiative Transfer (AT3D)
 
+AT3D performs 3D reconstruction of cloud/aerosol microphysical properties from multi-angle, multi-spectral solar reflected radiation using a non-linear optimization procedure [[1],[2],[3]].
+The core radiative transfer routines are sourced from the Fortran SHDOM (Spherical Harmonic Discrete Ordinate Method for 3D Atmospheric Radiative Transfer) code by Frank K. Evans [[4],[5]].
+
+The python package was created by [Aviad Levis](https://www.aviadlevis.info), Amit Aides (Technion - Israel Institute of Technology) and [Jesse Loveridge](https://atmos.illinois.edu/directory/profile/jesserl2) (University of Illinois). Code contributions have been made so far by Linda Forster and Vadim Holodovsky.
+
+## Usage
+
+The AT3D software is built around SHDOM, which is freely distributed online [[4]]. Please contact Frank Evans if you have concerns about the licensing of his code as it appears in this package. This package (AT3D) is distributed under the GNU General Public License (see the `LICENSE` file).
+If you want to acknowledge the use of this repository in a publication (e.g. scientific journal article), then please cite the appropriate release, or the most recent release, which is available at the following DOI. See the `CITATION.cff` file for how to reference this repository.
+
 [![DOI](https://zenodo.org/badge/342386439.svg)](https://zenodo.org/badge/latestdoi/342386439)
 
-AT3D performs 3D reconstruction of cloud/aerosol microphysical properties from multi-angle, multi-spectral solar reflected radiation using a non-linear optimization procedure [[1],[2],[3]].
-The core radiative transfer routines are sourced from the Fortran SHDOM (Spherical Harmonic Discrete Ordinate Method for 3D Atmospheric Radiative Transfer) code by Frank K. Evans [[4]].
+If you want to acknowledge the scientific origin of a particular feature of this software in a publication, then please cite the appropriate journal or conference articles in which the feature originates [[1],[2],[3]]. This work relies on the generosity of Frank Evans in making his code publicly available, for which we are very grateful. Please acknowledge his work appropriately. In particular, use of the SHDOM solver as a part of the AT3D software in a scientific work should cite [[5]].
 
-The python package was created by [Aviad Levis](https://www.aviadlevis.info), Amit Aides (Technion - Israel Institute of Technology) and [Jesse Loveridge](https://atmos.illinois.edu/directory/profile/jesserl2) (University of Illinois). Code contributions were made by Linda Forster and Vadim Holodovsky.
+Any publications using the synthetic les clouds in the ./data/synthetic_cloud_fields/jpl_les directory which is included in the distribution must cite the following work [[7]].
 
+## Contact
+
+If you find this package useful and/or would like to contribute code please let us know: aviad.levis@gmail.com; jesserl2@illinois.edu. 
 
 [1]: http://openaccess.thecvf.com/content_iccv_2015/html/Levis_Airborne_Three-Dimensional_Cloud_ICCV_2015_paper.html
 [2]: http://openaccess.thecvf.com/content_cvpr_2017/html/Levis_Multiple-Scattering_Microphysics_Tomography_CVPR_2017_paper.html
-[3]: https://www.mdpi.com/2072-4292/12/17/2831
+[3]: https://doi.org/10.3390/rs12172831
 [4]: http://coloradolinux.com/~evans/shdom.html
+[5]: https://doi.org/10.1175/1520-0469(1998)055<0429:TSHDOM>2.0.CO;2
+[6]: https://doi.org/10.1175/1520-0477(1998)079<0831:OPOAAC>2.0.CO;2
+[7]: https://doi.org/10.1175/JAS-D-13-0306.1
 
 &nbsp;
 
@@ -28,7 +43,7 @@ The key features of polarized SHDOM are included
 Note that each RTE solution is serial (**unlike SHDOM**) but independent wavelengths and pixel radiance calculations are parallelized using either MPI or a multi-threading shared memory framework.
 Other key features that are implemented are:
   * Several sensor configurations (e.g. Perspective, Orthographic) and arbitrary observation geometries.
-  * Mie & Rayleigh scattering optical property calculations. Optical properties of other species (e.g. non-spherical ice or aerosol or absorbing gases) can be included but must be calculated externally.
+  * Mie & Rayleigh scattering optical property calculations including [OPAC aerosols](6). Optical properties of other species (e.g. non-spherical ice or aerosol or absorbing gases) can be included but must be calculated externally.
   * Microphysical/optical properties can be generated or be read from netCDF or the SHDOM/[I3RC](https://i3rc.gsfc.nasa.gov/) file format.
 
 #### Inverse (remote-sensing):
@@ -41,7 +56,7 @@ Future improvement include:
 * Parallelize RTE solution with MPI.
 * Include retrieval of surface BRDF.
 
-To contribute to the development effort, contact us! see `Usage and Contact` section below.
+To contribute to the development effort, contact us! see `Contact` section above.
 
 &nbsp;
 
@@ -52,7 +67,9 @@ To contribute to the development effort, contact us! see `Usage and Contact` sec
 &nbsp;
 
 ## Installation
-Compilation of this package requires Fortran & C compilers (e.g. GCC 9.3.0_1) to be available and correctly linked. Installation has been tested on Mac and Linux using using [anaconda](https://www.anaconda.com/) package management.
+Compilation of this package requires Fortran & C compilers (e.g. GCC 9.3.0_1) to be available and correctly linked. Installation has been tested on Mac and Linux using using [anaconda](https://www.anaconda.com/) package management. 
+
+The treatment of legacy Fortran code has changed from GCC 9.X to 10.X+ so currently there is a flag in the setup.py script which needs to be commented if trying to install using GCC 9.X or earlier versions. The flag is `extra_f77_compile_args=["-fallow-argument-mismatch"]`. There is additional discussion of this point written as comments in setup.py but please feel raise an issue or discussion if you run into issues with the compiler version. The default version of the package should compile with GCC v11.3.
 
 Clone the repository into your local machine
 ```
@@ -110,10 +127,3 @@ For basic usage follow the following jupyter notebook tutorials under the notebo
 For generating rendering and optimization scripts see the list below.
 The scripts folder contains another readme file with examples of how to run each script.
 TODO
-
-
-&nbsp;
-
-## Usage and Contact
-If you find this package useful and/or would like to contribute code please let us know: aviad.levis@gmail.com; jesserl2@illinois.edu.  
-If you use this package in an academic publication please acknowledge the appropriate publications (see LICENSE file).

at3d/checks.py

@@ -290,9 +290,9 @@ def check_sensor(dataset):
                 " Valid values are 'I', 'Q', 'U', 'V'.".format(i)
                 )
     check_hasdim(dataset, stokes='stokes_index')
-    if dataset.stokes.dtype != np.bool:
+    if dataset.stokes.dtype != bool:
         raise TypeError("'stokes' variable in sensor dataset should be of boolean type.")
-    if dataset.use_subpixel_rays.dtype != np.bool:
+    if dataset.use_subpixel_rays.dtype != bool:
         raise TypeError("'use_subpixel_rays' variable in sensor dataset should be of boolean type.")
     if dataset.use_subpixel_rays.size != 1:
         raise ValueError("'use_subpixel_rays' variable should have a single value shape=(1,).")

at3d/containers.py

@@ -179,7 +179,7 @@ def get_measurements(self, solvers, n_jobs=1, mpi_comm=None, maxiter=100, verbos
                 "`solvers` should be of type '{}' not '{}'".format(
                     SolversDict, type(solvers))
                 )
-        if not isinstance(destructive, np.bool):
+        if not isinstance(destructive, bool):
             raise TypeError("`destructive` should be a boolean.")
 
         rte_sensors, sensor_mappings = self.sort_sensors(solvers)

at3d/gradient.py

@@ -84,7 +84,7 @@ def __init__(self, measurements, solvers, forward_sensors,
             uncertainty_kwargs['add_noise'] = False
             warnings.warn("'add_noise' flag was unspecified. It has been set to False.")
         else:
-            if not isinstance(uncertainty_kwargs['add_noise'], np.bool):
+            if not isinstance(uncertainty_kwargs['add_noise'], bool):
                 raise TypeError("uncertainty_kwargs['add_noise'] should be of boolean type.")
         if 'cost_function' not in gradient_kwargs:
             raise ValueError("'cost_function' must be specified in `gradient_kwargs`. Supported values "

at3d/solver.py

@@ -289,11 +289,11 @@ def solve(self, maxiter, init_solution=True, setup_grid=True, verbose=True,
         verbose: boolean
             True will output solution iteration information into stdout.
         """
-        if not isinstance(verbose, np.bool):
+        if not isinstance(verbose, bool):
             raise TypeError("`verbose` should be a boolean.")
-        if not isinstance(init_solution, np.bool):
+        if not isinstance(init_solution, bool):
             raise TypeError("`init_solution` should be a boolean.")
-        if not isinstance(setup_grid, np.bool):
+        if not isinstance(setup_grid, bool):
             raise TypeError("`setup_grid` should be a boolean.")
 
         # Part 1: Initialize solution (from a 1D layered model)
@@ -1832,9 +1832,9 @@ def _setup_numerical_params(self, numerical_params):
         self._transcut = numerical_params.transcut.data
         self._transmin = numerical_params.transmin.data
 
-        if self._deltam.dtype != np.bool:
+        if self._deltam.dtype != bool:
             raise TypeError("numerical_params.deltam should be of boolean type.")
-        if self._highorderrad.dtype != np.bool:
+        if self._highorderrad.dtype != bool:
             raise TypeError("numerical_params.high_order_radiance should be of boolean type.")
 
         if self._nphi < self._nmu:
@@ -1979,7 +1979,8 @@ def ibits(val, bit, ret_val):
             delyp=self._pa.dely,
             zlevels=self._pa.zlevels
         )
-
+        self._xgrid = self._xgrid[:self._nx1]
+        self._ygrid = self._ygrid[:self._ny1]
         self._setup_memory()
 
         self._npts, self._ncells, self._gridpos, self._gridptr, self._neighptr, \

at3d/util.py

@@ -523,8 +523,8 @@ def load_forward_model(file_name):
         for i, image in sensor.groups.items():
             sensor_dataset = xr.open_dataset(xr.backends.NetCDF4DataStore(dataset[
                 'sensors/'+str(key)+'/'+str(i)]))
-            sensor_dataset['stokes'] = (['stokes_index'],sensor_dataset['stokes'].data.astype(np.bool))
-            sensor_dataset['use_subpixel_rays'] = sensor_dataset['use_subpixel_rays'].data.astype(np.bool)
+            sensor_dataset['stokes'] = (['stokes_index'],sensor_dataset['stokes'].data.astype(bool))
+            sensor_dataset['use_subpixel_rays'] = sensor_dataset['use_subpixel_rays'].data.astype(bool)
             sensor_dict.add_sensor(key, sensor_dataset)
         #     sensor_list.append(xr.open_dataset(xr.backends.NetCDF4DataStore(dataset[
         #         'sensors/'+str(key)+'/'+str(i)])))
@@ -533,8 +533,8 @@ def load_forward_model(file_name):
     for key, solver in solvers.items():
         numerical_params = xr.open_dataset(xr.backends.NetCDF4DataStore(dataset[
                         'solvers/'+str(key)+'/numerical_parameters']))
-        numerical_params['deltam'] = numerical_params.deltam.data.astype(np.bool)
-        numerical_params['high_order_radiance'] = numerical_params.high_order_radiance.data.astype(np.bool)
+        numerical_params['deltam'] = numerical_params.deltam.data.astype(bool)
+        numerical_params['high_order_radiance'] = numerical_params.high_order_radiance.data.astype(bool)
         num_stokes = numerical_params.num_stokes.data
         surface = xr.open_dataset(xr.backends.NetCDF4DataStore(dataset['solvers/'+str(key)+'/surface']))
         source = xr.open_dataset(xr.backends.NetCDF4DataStore(dataset['solvers/'+str(key)+'/source']))

tests/test_derivatives.py

@@ -4,10 +4,15 @@
 import xarray as xr
 import at3d
 from scipy import stats
-
+import sys
 import warnings
 warnings.filterwarnings('ignore')
 
+import builtins as __builtin__
+def print(*args, **kwargs):
+    if '-vv' in sys.argv:
+        return __builtin__.print(*args, **kwargs)
+
 class PlanckDerivative(TestCase):
 
     @classmethod

tests/test_direct_beam.py

@@ -5,11 +5,16 @@
 import numpy as np
 import xarray as xr
 import at3d
-from scipy import stats
+import sys
 
 import warnings
 warnings.filterwarnings('ignore')
 
+import builtins as __builtin__
+def print(*args, **kwargs):
+    if '-vv' in sys.argv:
+        return __builtin__.print(*args, **kwargs)
+
 def cloud_direct_beam(mie_mono_table,ext,veff,reff,ssalb,solarmu,surfacealb, ground_temperature, step=0.0, index=(1,1,1),
           nmu=16, split=0.03, load_solution=None, resolution=1, boundary='open',deltam=True):
     rte_grid = at3d.grid.make_grid(0.05/resolution, 8*resolution, 0.05/resolution, 3*resolution,

tests/test_shdom.py

@@ -4,9 +4,16 @@
 import numpy as np
 import xarray as xr
 import at3d
+import sys
 
 warnings.filterwarnings('ignore')
 
+import builtins as __builtin__
+def print(*args, **kwargs):
+    if '-vv' in sys.argv:
+        return __builtin__.print(*args, **kwargs)
+
+
 def parse_shdom_output(filename, comment='*'):
     """
     Parse output from original SHDOM files
@@ -471,7 +478,7 @@ def test_no_subpixel_args_reference(self):
         print(np.abs(diff).max())
         print(np.argmax(np.abs(diff)), diff.shape)
         print(np.sqrt(np.mean(diff**2)))
-        self.assertTrue(test.equals(self.Sensordict['MISR']['sensor_list'][0]))
+        xr.testing.assert_allclose(test, self.Sensordict['MISR']['sensor_list'][0], rtol=1e-5)
 
 
 class Parallelization_SubpixelRays(TestCase):
@@ -584,7 +591,7 @@ def test_subpixel_args_reference(self):
         print(np.abs(diff).max())
         print(np.argmax(np.abs(diff)), diff.shape)
         print(np.sqrt(np.mean(diff**2)))
-        self.assertTrue(test.equals(self.Sensordict['MISR']['sensor_list'][0]))
+        xr.testing.assert_allclose(test, self.Sensordict['MISR']['sensor_list'][0], rtol=1e-5)
 
 
 class Verify_Lambertian_Surfaces(TestCase):
@@ -969,7 +976,6 @@ def setUpClass(cls):
         cls.rad = rad
 
     def test_radiance(self):
-
         #print(self.rad,self.rad2, self.integrated_rays.I.data)
         print(100*np.max(np.abs(self.rad - self.integrated_rays.I.data)/self.rad))
         print(np.max(np.abs(self.rad - self.integrated_rays.I.data)))
@@ -1044,7 +1050,6 @@ def setUpClass(cls):
         cls.rad = rad
 
     def test_radiance(self):
-
         #print(self.rad,self.rad2, self.integrated_rays.I.data)
         print(100*np.max(np.abs(self.rad - self.integrated_rays.I.data)/self.rad))
         print(np.max(np.abs(self.rad - self.integrated_rays.I.data)))