Merge pull request #146 from geocryology/grid-elev-in-scaled

Include grid elevation in scaled results

globsim/_version.py

@@ -1,3 +1,3 @@
-__version__ = "3.9.1"
+__version__ = "3.9.2"
 
 

globsim/scale/ERA5scale.py

@@ -122,6 +122,9 @@ def process(self):
                                     self.nc_pl_sur, self.times_out_nc,
                                     t_unit=self.scaled_t_units,
                                     station_names=stations)
+        # add surface height
+        self.add_grid_elevation(self.rg, self.getValues(self.nc_to, 'z')[0, :] / const.G)
+
         self.indProcess()
 
         logger.info(f"Created scaled output file {self.output_file}")

globsim/scale/GenericScale.py

@@ -201,6 +201,25 @@ def run_kernels(self):
     def _rh(self) -> Callable:
         rh_function = getattr(met, self._rh_function_name)
         return rh_function
+
+    def add_grid_elevation(self, ncf: "nc.Dataset", data: np.ndarray) -> None:
+        """Add station elevation to the netCDF file"""
+        elev_var = ncf.createVariable('grid_elevation', 'f4', ('station',))
+        elev_var.units = 'm'
+        elev_var.long_name = 'grid elevation'
+        elev_var.comment = 'Elevation of the grid cell at the station location'
+        elev_var[:] = data
+
+        self.warn_station_elevation(data)
+
+    def warn_station_elevation(self, data: "np.ndarray") -> None:
+        """Warn if there are stations with elevation below grid level"""
+        stn_elev = self.stations['elevation_m']
+        stn_name = self.stations['station_name']
+
+        for i, (stn, grid) in enumerate(zip(stn_elev, data)):
+            if stn < grid:
+                logger.warning(f" {stn_name[i]} site elevation ({stn} m) is below reanalysis grid elevation ({grid} m). Results may be unreliable.")
 
 
 def _check_timestep_length(nctime: "nc.Variable", source:str) -> None:

globsim/scale/JRAscale.py

@@ -86,6 +86,10 @@ def process(self):
         self.rg = new_scaled_netcdf(self.output_file, self.nc_pl_sur,
                                     self.times_out_nc, self.nc_pl_sur['time'].units)
 
+        # add station names to netcdf
+        self.add_grid_elevation(self.rg, self.get_values("to", "Geopotential", 
+                                                         (0, slice(None,None,1))) / const.G)  # [m]
+
         # add station names to netcdf
         # first convert to character array
         names_out = nc.stringtochar(np.array(self.stations['station_name'], 'S32'))

globsim/scale/MERRAscale.py

@@ -182,7 +182,10 @@ def process(self):
                                     self.nc_pl_sur,
                                     self.times_out_nc,
                                     t_unit=self.scaled_t_units)
-
+
+        # add grid elevation to netCDF
+        self.add_grid_elevation(self.rg, self.nc_sc['PHIS'][0, :] / const.G)
+
         # add station names to netcdf
         # first convert to character array
         names_out = nc.stringtochar(np.array(self.stations['station_name'], 'S32'))