raster compat w/RAS2025, ditch multiprocessing, more RAS reading

Author: Rusty Holleman

stompy/model/openfoam/depth_average.py

@@ -7,7 +7,7 @@
 @author: chrhol4587
 """
 import numpy as np
-import os
+import os, time
 import hashlib, pickle
 from multiprocessing import Pool
 import subprocess
@@ -364,26 +364,30 @@ def precalc_raster_info_bbox(self,fld):
                     raster_weights[pix,cell] = A[cell]/Apixel
         self.raster_precalc = raster_weights
 
-        
     def precalc_raster_info_clipping(self, fld, force=False):
         # weights is a matrix with columns corresponding to openfoam cells
         # and rows corresponding to output pixels in row-major order
 
         raster_weights = None
         tasks = [ (self.proc_dir(proc),fld) 
                  for proc in range(self.n_procs) ]
-              
-        with Pool(self.n_tasks) as pool: 
-            
-            results = pool.starmap(precalc_raster_weights_proc,tasks)
+
+        # fallback to serial...       
+        if 0:
+            with Pool(self.n_tasks) as pool: 
+                
+                results = pool.starmap(precalc_raster_weights_proc,tasks)
+        else:
+            results = [precalc_raster_weights_proc(*t)
+                       for t in tasks]
 
-            for proc,proc_raster_weights in enumerate(results):
-                print(f"Assembling from processor {proc}")
-                # Want to stack these left-to-right
-                if raster_weights is None:
-                    raster_weights = proc_raster_weights
-                else:
-                    raster_weights = sparse.hstack( (raster_weights,proc_raster_weights) )
+        for proc,proc_raster_weights in enumerate(results):
+            print(f"Assembling from processor {proc}")
+            # Want to stack these left-to-right
+            if raster_weights is None:
+                raster_weights = proc_raster_weights
+            else:
+                raster_weights = sparse.hstack( (raster_weights,proc_raster_weights) )
         self.raster_precalc = raster_weights
 
     def to_raster(self, variable, timename):
@@ -476,7 +480,10 @@ def face_centers(self,proc=None):
 
     def compute_cell_centers(self,center_fn):
         case = os.path.dirname(center_fn)
-        completed = subprocess.run([self.writeMeshObj,"-constant","-time","none"], cwd=case)
+        try:
+            completed = subprocess.run([self.writeMeshObj,"-constant","-time","none"], cwd=case)
+        except FileNotFoundError:
+            return False
         if completed.returncode!=0:
             return False
         if not os.path.exists(center_fn):

stompy/model/ras/result_reader.py

@@ -6,9 +6,11 @@
 """
 
 from ...grid import unstructured_grid
-from ... import utils
+from ... import utils, memoize
 import h5py
 import numpy as np
+import datetime
+import glob, os
 
 # Would like to replicate what RAS does
 # technically that requires the terrain
@@ -61,6 +63,54 @@ def face_interp_ben(wse_j1, wse_j2, z_j1, z_j2, face_min_z):
           face_wse = min(face_wse, max_wse)
     return face_wse
 
+class RasProject:
+    def __init__(self,prj_fn):
+        self.prj_fn = prj_fn
+        self.prj_dir = os.path.dirname(prj_fn)
+        self.prj_name = os.path.basename(prj_fn).replace(".prj","")
+        
+    def getPlanFile(self, plan_name):
+        '''
+        reads the plan file from the plan HDF5 file and finds correct corresponding plan file to the short ID
+        '''
+        plan_files = self.getPlanFiles()
+        h=None
+        for plf in plan_files:
+            if plf.endswith('.tmp.hdf'):
+                continue
+            try:
+                h = h5py.File(os.path.join(self.prj_dir, plf), 'r')
+                try:
+                    short_title = h['Results/Unsteady'].attrs['Short ID'].decode('utf-8')
+                finally:
+                    h.close()
+                if short_title == plan_name:
+                    print(f'Found plan file {plf} for plan: {plan_name}')
+                    return plf
+            except Exception as e:
+                print('Error reading plan file:', plf)
+                # print(e)
+                continue
+
+        print('Unable to find plan file for plan:', plan_name)
+        #print('Please specficy plan file using variable plan_file=... in inputs.')
+        return None
+    
+    def getPlanFiles(self):
+        '''
+        gets all the plan files in the project directory
+        '''
+        plan_file_form = f'{self.prj_name}.p*.hdf' # "*.p[!*.]*"
+        plan_files = glob.glob(os.path.join(self.prj_dir, plan_file_form))
+        return plan_files
+
+    @memoize.imemoize()
+    def get_results(self,plan_name,area_name=None):
+        plan_file = self.getPlanFile(plan_name)
+        if plan_file is None:
+            return None
+        return RasReader(plan_file, area_name=area_name)
+
 class RasReader:
     # Where to find results within h5 file
     unsteady_base_ras6=('/Results/Unsteady/Output/Output Blocks/'
@@ -73,6 +123,10 @@ def __init__(self,results_fn, area_name=None):
 
         self.load_h5()
         self.load_grid()
+        
+    def short_title(self):
+        return self.h5['Results/Unsteady'].attrs['Short ID'].decode('utf-8')
+        
     def load_h5(self):
         self.h5 = h5py.File(self.results_fn,'r')
         # Is this RAS6 or RAS2025?
@@ -81,6 +135,10 @@ def load_h5(self):
             self.version="RAS6"
         except KeyError:
             self.version="RAS2025"
+            
+    def close(self):
+        self.h5.close()
+        self.h5 = None
 
     def terrain_file(self):
         key = f"/Geometry/2D Flow Areas/{self.area_name}"
@@ -113,8 +171,19 @@ def area_base(self):
         else:
             raise Exception("Bad version: "+self.version)
 
-    def time_days(self):
+    def time_relative_days(self):
         return self.h5[self.unsteady_base+'/Time']
+
+    def time_start(self):
+        t = self.h5['Plan Data/Plan Information'].attrs['Simulation Start Time'].decode('ascii')
+        # '22Aug2022 02:00:00'
+        t_datetime = datetime.datetime.strptime(t,"%d%b%Y %H:%M:%S")
+        return np.datetime64(t_datetime)
+
+    def times(self):
+        t0 = self.time_start()
+        offset_days = self.time_relative_days()[:]
+        return t0 + np.array(offset_days*86400,np.int64)*np.timedelta64(1,'s')
 
     def cell_wse(self,time_step,trim_virtual=True):
         key=self.area_base+'Water Surface'

stompy/spatial/rasterize_stl.py

@@ -63,10 +63,17 @@ def stl_to_field(fp_stl,R,translate=[0,0,0],dx=None,dy=None,max_dim=None,pad=[0,
     else:
         # with rounding there's some one-off in here, but I'm not worrying
         # about that right now.
-        dx=max( xxyy[1]-xxyy[0], xxyy[3]-xxyy[2] ) / float(max_dim)
+        # try to be consistent with SimpleGrid.delta()
+        # dx = (xmax-xmin)/(self.F.shape[1]-1.0)
+        # i.e. xmin and xmax reference pixel centers
+        # 
+        dx=max( xxyy[1]-xxyy[0], xxyy[3]-xxyy[2] ) / (float(max_dim)-1.0)
         dy=dx
-        
-    Z = np.full( (int( (xxyy[3]-xxyy[2])/dx),int( (xxyy[1]-xxyy[0])/dx)), np.float64(nodata))
+
+    # 1.0 to account for pixel centers, and 0.5 so int() will round to nearest 
+    Z = np.full( (int(1.5 + (xxyy[3]-xxyy[2])/dx),
+                  int(1.5 + (xxyy[1]-xxyy[0])/dx)),
+                 np.float64(nodata))
     print("Output shape will be ",Z.shape)
 
     fld = field.SimpleGrid(extents=xxyy,F=Z)
@@ -161,10 +168,13 @@ def main():
         kwargs['dy']=args.res
     else:
         kwargs['max_dim']=args.size
-    fld = stl_to_field(fp_stl, R, translate=translate, nodata=args.nodata,
+    fld = stl_to_field(fp_stl, R, translate=translate, 
                        pad=args.pad, **kwargs)
     fld.assign_projection(args.proj)
-    fld.write_gdal(args.tif_file,overwrite=args.force)
+    fld.F = fld.F.astype(np.float32) # RAS2025 doesn't like doubles
+    # RAS 6.6 defaults to DEFLATE. Stick with that for better compatibility
+    fld.write_gdal(args.tif_file,overwrite=args.force, nodata=args.nodata,
+                   options=["COMPRESS=DEFLATE"])
 
     if args.plot:    
         import matplotlib.pyplot as plt
@@ -181,3 +191,6 @@ def main():
 if __name__ == '__main__':
     main()
 
+
+
+