Merge pull request #1216 from cbworden/realz_data

Realz data

Author: cbworden

azure-pipelines.yml

@@ -28,6 +28,12 @@ jobs:
           #MacOS_11_Python38:
           #imageName: 'macos-11'
           #python.version: '3.8'
+      MacOS_12_Python39:
+        imageName: 'macos-latest'
+        python.version: '3.9'
+      MacOS_12_Python310:
+        imageName: 'macos-latest'
+        python.version: '3.10'
       MacOS_11_Python39:
         imageName: 'macos-11'
         python.version: '3.9'
@@ -68,15 +74,18 @@ jobs:
           echo "MacOS 10.14"
         else
           export CPATH=/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include
-          echo "MacOS 10.15 or MacOS 11"
+          echo "MacOS 10.15 or MacOS 11 or MacOS 12"
         fi
       fi
       echo $(Agent.NAME)
       echo $(python.version)
+      echo "Conda is " $CONDA " user is " $USER
       bash install.sh -p $(python.version)
     displayName: Create environment
 
-  - bash: conda init bash                            
+  - bash: |
+      echo "##vso[task.prependpath]$CONDA/bin"
+      conda init bash
     displayName: Init conda for bash                                         
 
   - bash: |

install.sh

@@ -157,6 +157,7 @@ package_list=(
       "python=$py_ver"
       "$CC_PKG"
       "cartopy>=0.20"
+      "clang>=12.0.1"
       "configobj>=5.0.6"
       "cython>=0.29.23"
       "defusedxml>=0.7.1"
@@ -178,11 +179,11 @@ package_list=(
       "ps2ff>=1.5.2"
       "psutil>=5.6.7"
       "pyproj>=2.6.1"
-      "pytest>=6.2.4"
-      "pytest-cov>=2.12.1"
+      "pytest==6.2.4"
+      "pytest-azurepipelines==0.8.0"
+      "pytest-cov==2.12.1"
+      "pytest-faulthandler==2.0.1"
       "python-daemon>=2.3.0"
-      "pytest-faulthandler>=2.0.1"
-      "pytest-azurepipelines>=0.8.0"
       "pyzmq"
       "rasterio>=1.2.5"
       "scikit-image>=0.16.2"

setup.py

@@ -1,7 +1,8 @@
 import os
 from distutils.core import setup
 import os.path
-import versioneer
+
+# import versioneer
 from distutils.extension import Extension
 from Cython.Distutils import build_ext
 from Cython.Build import cythonize
@@ -13,7 +14,7 @@
 if shutil.which("clang") is None:
     os.environ["CC"] = "gcc"
 else:
-    os.environ["CC"] = "clang"
+    os.environ["CC"] = "/usr/bin/clang"
 
 sourcefiles = ["shakemap/c/pcontour.pyx", "shakemap/c/contour.c"]
 
@@ -37,12 +38,12 @@
     ),
 ]
 
-cmdclass = versioneer.get_cmdclass()
-cmdclass["build_ext"] = build_ext
+# cmdclass = versioneer.get_cmdclass()
+# cmdclass["build_ext"] = build_ext
 
 setup(
     name="shakemap",
-    version=versioneer.get_version(),
+    #    version=versioneer.get_version(),
     description="USGS Near-Real-Time Ground Motion Mapping",
     author="Bruce Worden, Mike Hearne, Eric Thompson",
     author_email="[email protected],[email protected],[email protected]",
@@ -93,6 +94,7 @@
         "bin/sm_sync",
         "bin/fix_netcdf",
     ],
-    cmdclass=cmdclass,
+    #    cmdclass=cmdclass,
+    cmdclass={"build_ext": build_ext},
     ext_modules=cythonize(ext_modules),
 )

shakemap/coremods/model.py

@@ -1345,7 +1345,7 @@ def _computeBias(self):
             self.ccf.getCorrelation(t1_22, t2_22, matrix22)
             sta_phi_flat = sta_phi.flatten()
             make_sigma_matrix(matrix22, sta_phi_flat, sta_phi_flat)
-            np.fill_diagonal(matrix22, np.diag(matrix22) + sta_sig_extra**2)
+            np.fill_diagonal(matrix22, np.diag(matrix22) + sta_sig_extra ** 2)
             cov_WD_WD_inv = np.linalg.pinv(matrix22)
             #
             # Hold on to some things we'll need later
@@ -1508,6 +1508,13 @@ def _computeMVN(self, imtstr):
             self.outsd[imtstr] = pout_sd[0]
             self.outphi[imtstr] = self.psd[imtstr]
             self.outtau[imtstr] = self.tsd[imtstr]
+            # Special stuff for the MMI priors. When we make this apply to other
+            # IMTs we'll need to mess around with this block
+            if imtstr == "MMI":
+                self.MMI_add_uncertainty = np.array([])
+                self.MMI_Sigma_HH_YD = np.array([])
+                self.MMI_C = np.array([])
+                self.MMI_sta_per_ix = np.array([])
             return
 
         pout_sd2_phi = np.power(self.psd[imtstr], 2.0)
@@ -1554,6 +1561,14 @@ def _computeMVN(self, imtstr):
         sdsta_phi = sta_phi.flatten()
         matrix12_phi = np.empty(t2_12.shape, dtype=np.double)
         rcmatrix_phi = np.empty(t2_12.shape, dtype=np.double)
+        # Allocate the full C matrix only for the desired IMT (MMI). This
+        # is for generating realizations. Someday we will want to make this
+        # apply to other IMTs and become an optional flag for this module
+        # where we save the priors optionally
+        # Note: We would have to set up C_complete = {} at the top of the
+        # module, and then here, it would be C_complete[imtstr] = ...
+        if imtstr == "MMI":
+            C_complete = np.empty_like(T_Y0)
         for iy in range(self.smny):
             ss = iy * self.smnx
             se = (iy + 1) * self.smnx
@@ -1631,6 +1646,11 @@ def _computeMVN(self, imtstr):
             sdgrid_tau[iy, :] = np.sum(
                 np.multiply(C_tmp1, C, out=C_tmp2), out=s_tmp3, axis=1
             )
+            # This is where we would have a flag if we were saving the priors
+            # as an option
+            if imtstr == "MMI":
+                # Save C to complete full size C
+                C_complete[ss:se, :] = C
 
             mtime += time.time() - time4
         #
@@ -1655,6 +1675,14 @@ def _computeMVN(self, imtstr):
         self.outphi[imtstr] = self.psd[imtstr]
         self.outtau[imtstr] = np.sqrt(sdgrid_tau, out=sdgrid_tau)
 
+        # Special stuff for the MMI priors. When we make this apply to other
+        # IMTs we'll need to mess around with this block
+        if imtstr == "MMI":
+            self.MMI_add_uncertainty = self.sta_sig_extra[imtstr]
+            self.MMI_Sigma_HH_YD = self.cov_HH_yD[imtstr]
+            self.MMI_C = C_complete
+            self.MMI_sta_per_ix = sta_per_ix
+
         self.logger.debug(f"\ttime for {imtstr} distance={ddtime:f}")
         self.logger.debug(f"\ttime for {imtstr} correlation={ctime:f}")
         self.logger.debug(f"\ttime for {imtstr} sigma={stime:f}")
@@ -2108,7 +2136,7 @@ def _fillStationJSON(self):
                 else:
                     mytau = sdf[key + "_tau"][six]
                 myphi = sdf[key + "_phi"][six]
-                mysigma = np.sqrt(mytau**2 + myphi**2)
+                mysigma = np.sqrt(mytau ** 2 + myphi ** 2)
                 mysigma_rock = sdf[key + "_sigma_rock"][six]
                 mysigma_soil = sdf[key + "_sigma_soil"][six]
                 imt_name = key.lower().replace("_pred", "")
@@ -2305,6 +2333,15 @@ def _storeGriddedData(self, oc):
                 self.outphi[key],
                 self.outtau[key],
             )
+            if key == "MMI":
+                sub_groups = ["imts", component, key]
+                oc.setArray(sub_groups, "add_uncertainty", self.MMI_add_uncertainty)
+                oc.setArray(sub_groups, "Sigma_HH_YD", self.MMI_Sigma_HH_YD)
+                oc.setArray(sub_groups, "C", self.MMI_C)
+                oc.setArray(sub_groups, "sta_per_ix", self.MMI_sta_per_ix)
+                oc.setArray(sub_groups, "sta_phi", self.sta_phi["MMI"])
+                oc.setArray(sub_groups, "sta_lons_rad", self.sta_lons_rad["MMI"])
+                oc.setArray(sub_groups, "sta_lats_rad", self.sta_lats_rad["MMI"])
         #
         # Directivity
         #
@@ -2561,7 +2598,7 @@ def _adjustResolution(self):
             target_res = (
                 -(latspan + lonspan)
                 - np.sqrt(
-                    latspan**2 + lonspan**2 + 2 * latspan * lonspan * (2 * nmax - 1)
+                    latspan ** 2 + lonspan ** 2 + 2 * latspan * lonspan * (2 * nmax - 1)
                 )
             ) / (2 * (1 - nmax))
 

utils/hdf2dict.py

@@ -41,6 +41,23 @@ def get_result_dict(filename):
         'ymax': the maximum latitude of the grid
         'ymin': the minimum latitude of the grid
     The grids themselves will be of dimension (ny x nx) (rows x cols).
+    In addition to the above grids, the following arrays will be returned:
+        'Sigma_HH_YD': A numpy array, see below for possible values.
+        'Sigma_HH_YD_metadata': This will always be 'None'.
+        'C': A numpy array, see below for possible values.
+        'C_metadata': This will always be 'None'.
+        'add_uncertainty': A numpy array, see below for possible values.
+        'add_uncertainty_metadata': This will always be 'None'.
+        'sta_per_ix': A numpy array, see below for possible values.
+        'sta_per_ix_metadata': This will always be 'None'.
+    These arrays will have the following contents:
+        'None': If this function is run on an older version of shake_results.hdf
+                that was produced by model.py before it generated the output
+                needed to produce these arrays.
+        '[]': I.e., empty arrays. When the shakemap contained no input station data
+              sufficient to generate the arrays.
+        ...: Arrays of shapes appropriate to the input and output data dimensions.
+             It is incumbent upon the user to make proper use of these arrays.
     """
     ddict = {}
     hdfobj = h5py.File(filename, "r+")
@@ -52,6 +69,21 @@ def get_result_dict(filename):
         for key, value in dset.attrs.items():
             metadata[key] = value
         ddict[name + "_metadata"] = metadata
+    for name in (
+        "add_uncertainty",
+        "Sigma_HH_YD",
+        "C",
+        "sta_per_ix",
+        "sta_phi",
+        "sta_lons_rad",
+        "sta_lats_rad",
+    ):
+        ddict[name + "_metadata"] = None
+        try:
+            dset = group[name]
+            ddict[name] = dset[()].astype(np.float64)
+        except Exception:
+            ddict[name] = None
 
     hdfobj.close()
     return ddict
@@ -63,6 +95,24 @@ def get_result_dict(filename):
 
     print("ddict keys: ", ddict.keys())
     print()
+    for key in ddict.keys():
+        if ddict[key] is None:
+            print(f"ddict {key} is None")
+    for name in ("mean", "std", "phi", "tau"):
+        print(f"{name} shape is {(ddict[name].shape)}")
+    for key in (
+        "add_uncertainty",
+        "Sigma_HH_YD",
+        "C",
+        "sta_per_ix",
+        "sta_phi",
+        "sta_lons_rad",
+        "sta_lats_rad",
+    ):
+        if ddict[key] is None:
+            print(f"{key} shape is None")
+        else:
+            print(f"{key} shape is {(ddict[key].shape)}")
     print("Metadata:")
     print(ddict["mean_metadata"].items())
     print()