Merge pull request #154 from a2e-mmc/dev

Merge dev branch

README.md

@@ -82,8 +82,29 @@ df = read_dir(dpath, file_filter='*_w*', reader=profiler)
 
 ## Installation
 
-To install, run `pip install -e mmctools` after cloning the repository (or `pip install -e .` from inside a2e-mmc/mmctools).
-
+The recommended approach is to first create a new conda environment:
+```
+conda create -n mmc python=3.7
+conda activate mmc
+conda install -y -c conda-forge jupyterlab matplotlib scipy xarray dask pyarrow gdal rasterio elevation pyyaml netcdf4 wrf-python cdsapi cfgrib
+```
+Note: All packages after `xarray` are optional:
+- `dask` makes netcdf data processing more efficient
+- `pyarrow` is a dependency for the "feather" data format, an *extremely* efficient
+  way to save dataframe data (in terms file I/O time and file size)
+- `gdal`, `rasterio`, and `elevation` are required for processing terrain data
+- `netcdf4` and `wrf-python` are for the NCAR-provided WRF utilities, which are
+  useful for interpolating and slicing data
+- `cdsapi` is needed for `wrf.preprocessing` to retrieve Copernicus ERA5
+  reanalysis data
+- `cfgrib` enables xarray to load grib files
+
+
+Then create an "editable" installation of the mmctools repository:
+```
+cd /path/to/a2e-mmc/mmctools
+pip install -e .
+```
 
 ## Code Development Principles
 

environment.yaml

@@ -1,16 +1,17 @@
 channels:
   - conda-forge
 dependencies:
-  - python=3.7
-  - matplotlib=3.1.3
-  - numpy=1.18.1
-  - scipy=1.4.1
-  - pandas=1.0.1
-  - xarray=0.15.0
-  - dask=2.11.0
-  - netcdf4=1.5.3
-  - wrf-python=1.3.2
-  - utm=0.5.0
-  - f90nml=1.1.2
-  - elevation=1.0.6
-  - rasterio=1.0.25
+  - python>=3.10.6
+  - matplotlib>=3.5.3
+  - numpy>=1.23.2
+  - scipy>=1.9.0
+  - pandas>=1.4.3
+  - xarray>=2022.6.0
+  - dask>=2022.8.1
+  - netcdf4>=1.6.0
+  - wrf-python>=1.3.4
+  - utm>=0.7.0
+  - f90nml>=1.4.3
+  - elevation>=1.1.3
+  - rasterio>=1.3.2
+  - richdem>=2.3.0

example_plotting.ipynb

@@ -24,17 +24,24 @@
     "import pandas as pd"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note: the following cell should not be needed if a `pip install [-e]` was performed"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
-    "#Make sure a2e-mmc repositories are in the pythonpath\n",
-    "a2epath = '/home/equon/a2e-mmc'\n",
-    "import sys\n",
-    "if not a2epath in sys.path:\n",
-    "    sys.path.append(a2epath)"
+    "# #Make sure a2e-mmc repositories are in the pythonpath\n",
+    "# a2epath = '/home/equon/a2e-mmc'\n",
+    "# import sys\n",
+    "# if not a2epath in sys.path:\n",
+    "#     sys.path.append(a2epath)"
    ]
   },
   {
@@ -85,7 +92,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "datadir = '/Users/equon/a2e-mmc/assessment/datasets/SWiFT/data'\n",
+    "datadir = '/home/equon/a2e-mmc/assessment/datasets/SWiFT/data'\n",
     "TTUdata = 'TTU_tilt_corrected_20131108-09.csv'"
    ]
   },
@@ -119,9 +126,8 @@
        "      <th>u</th>\n",
        "      <th>v</th>\n",
        "      <th>w</th>\n",
-       "      <th>Ts</th>\n",
-       "      <th>T</th>\n",
-       "      <th>RH</th>\n",
+       "      <th>t</th>\n",
+       "      <th>ts</th>\n",
        "      <th>p</th>\n",
        "    </tr>\n",
        "    <tr>\n",
@@ -133,7 +139,6 @@
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
-       "      <th></th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
@@ -143,63 +148,58 @@
        "      <td>-0.138929</td>\n",
        "      <td>2.637817</td>\n",
        "      <td>0.074016</td>\n",
-       "      <td>289.410000</td>\n",
        "      <td>284.794</td>\n",
-       "      <td>26.186</td>\n",
+       "      <td>289.410000</td>\n",
        "      <td>908.547754</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2.4</th>\n",
        "      <td>-0.601111</td>\n",
        "      <td>2.783204</td>\n",
        "      <td>0.487330</td>\n",
-       "      <td>290.979994</td>\n",
        "      <td>284.932</td>\n",
-       "      <td>25.810</td>\n",
+       "      <td>290.979994</td>\n",
        "      <td>908.723508</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>4.0</th>\n",
        "      <td>0.416792</td>\n",
        "      <td>4.043940</td>\n",
        "      <td>0.295800</td>\n",
-       "      <td>287.520000</td>\n",
        "      <td>285.166</td>\n",
-       "      <td>25.380</td>\n",
+       "      <td>287.520000</td>\n",
        "      <td>908.215548</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>10.1</th>\n",
        "      <td>-0.276479</td>\n",
        "      <td>5.227110</td>\n",
        "      <td>-0.418065</td>\n",
-       "      <td>287.250000</td>\n",
        "      <td>285.298</td>\n",
-       "      <td>25.264</td>\n",
+       "      <td>287.250000</td>\n",
        "      <td>907.611414</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>16.8</th>\n",
        "      <td>0.034364</td>\n",
        "      <td>5.908367</td>\n",
        "      <td>-0.173836</td>\n",
-       "      <td>287.610000</td>\n",
        "      <td>285.414</td>\n",
-       "      <td>24.934</td>\n",
+       "      <td>287.610000</td>\n",
        "      <td>907.307654</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "text/plain": [
-       "                          u         v         w          Ts        T      RH  \\\n",
-       "datetime   height                                                              \n",
-       "2013-11-08 0.9    -0.138929  2.637817  0.074016  289.410000  284.794  26.186   \n",
-       "           2.4    -0.601111  2.783204  0.487330  290.979994  284.932  25.810   \n",
-       "           4.0     0.416792  4.043940  0.295800  287.520000  285.166  25.380   \n",
-       "           10.1   -0.276479  5.227110 -0.418065  287.250000  285.298  25.264   \n",
-       "           16.8    0.034364  5.908367 -0.173836  287.610000  285.414  24.934   \n",
+       "                          u         v         w        t          ts  \\\n",
+       "datetime   height                                                      \n",
+       "2013-11-08 0.9    -0.138929  2.637817  0.074016  284.794  289.410000   \n",
+       "           2.4    -0.601111  2.783204  0.487330  284.932  290.979994   \n",
+       "           4.0     0.416792  4.043940  0.295800  285.166  287.520000   \n",
+       "           10.1   -0.276479  5.227110 -0.418065  285.298  287.250000   \n",
+       "           16.8    0.034364  5.908367 -0.173836  285.414  287.610000   \n",
        "\n",
        "                            p  \n",
        "datetime   height              \n",
@@ -235,7 +235,7 @@
    "source": [
     "# Calculate wind speed and direction\n",
     "df['wspd'], df['wdir'] = calc_wind(df)\n",
-    "df['theta'] = theta(df['T'],df['p'])"
+    "df['theta'] = theta(df['t'],df['p'])"
    ]
   },
   {
@@ -483,9 +483,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.1"
+   "version": "3.7.6"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }

mmctools/__init__.py

@@ -0,0 +1,24 @@
+
+#import windtools.windtools.plotting as plotting
+#
+# <module 'windtools.windtools.plotting' from '/home/equon/windtools/windtools/plotting.py'>
+# -- the wrong module was imported!
+
+#from .windtools.windtools import plotting
+#print(plotting)
+#
+# <module 'mmctools.windtools.windtools.plotting' from '/home/equon/a2e-mmc/mmctools/mmctools/windtools/windtools/plotting.py'>
+# ImportError: cannot import name 'plot_*' from 'mmctools.plotting' (unknown location)
+
+#from .windtools.windtools import *
+#print(plotting)
+#
+# NameError: name 'plotting' is not defined
+
+
+# for backwards compatibility, enable `from mmctools.plotting import plot_something`
+# enable `from mmctools.foo.bar import baz # to import baz from windtools.foo.bar`
+import os
+mmctools_module = os.path.split(__file__)[0]
+__path__.append(os.path.join(mmctools_module,'windtools','windtools'))
+

mmctools/coupling/sowfa.py

@@ -113,30 +113,50 @@ def write_BCs(self,
         =====
         fname : str
             Filename
-        fieldname : str
-            Name of the field to be written out
+        fieldname : str or list-like
+            Name of the scalar field (or a list of names of vector field
+            components) to be written out; 0 may be substituted to
+            indicate an array of zeroes
         fact : float
             Scale factor for the field, e.g., to scale heat flux to follow
             OpenFOAM sign convention that boundary fluxes are positive if
             directed outward
         """
 
-        # extract time and height array
+        # extract time array
         ts = self.df.t_index.values
         nt = ts.size
-
-        # assert field exists and is complete
-        assert(fieldname in self.df.columns), 'Field '+fieldname+' not in df'
-        assert(~pd.isna(self.df[fieldname]).any()), 'Field '+fieldname+' is not complete (contains NaNs)'
-
-        # scale field with factor,
-        # e.g., scale heat flux with fact=-1 to follow OpenFOAM sign convention
-        fieldvalues = fact * self.df[fieldname].values
-
-        with open(os.path.join(self.dpath,fname),'w') as fid:
+
+        # check if scalar or vector
+        if isinstance(fieldname, (list,tuple)):
+            assert len(fieldname) == 3, 'expected 3 vector components'
+            fieldnames = fieldname
+            fmt = ['    (%g', '(%.12g', '%.12g', '%.12g))',]
+        else:
+            fieldnames = [fieldname]
             fmt = ['    (%g', '%.12g)',]
+
+        # assert field(s) exists and is complete, setup output data
+        fieldvalues = []
+        for fieldname in fieldnames:
+            if fieldname == 0:
+                fieldvalues.append(np.zeros_like(ts))
+            else:
+                assert(fieldname in self.df.columns), \
+                        'Field '+fieldname+' not in df'
+                assert(~pd.isna(self.df[fieldname]).any()), \
+                        'Field '+fieldname+' is not complete (contains NaNs)'
+                fieldvalues.append(self.df[fieldname].values)
+
+                # scale field with factor,
+                # e.g., scale heat flux with fact=-1 to follow OpenFOAM sign convention
+                fieldvalues[-1] = fact * fieldvalues[-1]
+
+        fieldvalues = np.array(fieldvalues).T  # result: fieldvalues.shape==(nt, ndim)
+
+        with open(os.path.join(self.dpath,fname),'w') as fid:
             np.savetxt(fid,np.concatenate((ts.reshape((nt,1)),
-                                          fieldvalues.reshape((nt,1))
+                                          fieldvalues
                                           ),axis=1),fmt=fmt)
 
         return
@@ -172,7 +192,7 @@ def write_ICs(self,
             if not field in df.columns:
                 df.loc[:,field] = 0.0
 
-        # extract time and height array
+        # extract height array
         zs = df.height.values
         nz = zs.size