Merge pull request #20 from wpreimes/develop

Fix grib and era5-land grid

CHANGES.rst

@@ -6,6 +6,13 @@ Unreleased
 ==========
 -
 
+Version 0.6.1
+===========
+
+- Fix bug when creating 0.1 deg grid cells (floating point precision)
+- Missing variables in grib files are now replaced by empty images.
+- Read variable names from grib files from cfVarNameECMF instead of short_name field
+
 Version 0.6
 ===========
 

ecmwf_models/era5/reshuffle.py

@@ -191,4 +191,4 @@ def main(args):
 
 
 def run():
-    main(sys.argv[1:])
+    main(sys.argv[1:])

ecmwf_models/grid.py

@@ -52,7 +52,6 @@ def get_grid_resolution(lats, lons):
         lon_res = lons_res[0]
     return float(lat_res), float(lon_res)
 
-
 def ERA5_RegularImgLandGrid(res_lat=0.25, res_lon=0.25):
     """
     Uses the 0.25 DEG ERA5 land mask to create a land grid of the same size,
@@ -72,12 +71,10 @@ def ERA5_RegularImgLandGrid(res_lat=0.25, res_lon=0.25):
     land_mask = ds.variables['land'][:].flatten().filled(0.) == 1.
     land_points = np.ma.masked_array(global_grid.get_grid_points()[0], ~land_mask)
 
-    land_grid = global_grid.subgrid_from_gpis(land_points[~land_points.mask].filled())
+    land_grid = global_grid.subgrid_from_gpis(land_points[~land_points.mask].filled().astype('int'))
 
     return land_grid.to_cell_grid(5., 5.)
 
-
-
 def ERA_RegularImgGrid(res_lat=0.25, res_lon=0.25):
     """
     Create ECMWF regular cell grid
@@ -94,22 +91,24 @@ def ERA_RegularImgGrid(res_lat=0.25, res_lon=0.25):
     CellGrid : pygeogrids.CellGrid
         Regular, global CellGrid with 5DEG*5DEG cells
     """
-
-    lon = np.arange(0., 360. - res_lon / 2., res_lon)
-    lat = np.arange(90., -1. * 90. - res_lat / 2., -1. * res_lat)
-    lons_gt_180 = np.where(lon > 180.0)
-    lon[lons_gt_180] = lon[lons_gt_180] - 360.
+    # np.arange is not precise...
+    f_lon = (1. / res_lon)
+    f_lat = (1. / res_lat)
+    res_lon = res_lon * f_lon
+    res_lat = res_lat * f_lat
+    lon = np.arange(0., 360. * f_lon, res_lon)
+    lat = np.arange(90. * f_lat, -90 * f_lat - res_lat, -1 * res_lat)
+    lons_gt_180 = np.where(lon > (180. * f_lon))
+    lon[lons_gt_180] = lon[lons_gt_180] - (360. *f_lon)
 
     lon, lat = np.meshgrid(lon, lat)
 
-    glob_basic_grid = BasicGrid(lon.flatten(), lat.flatten())
+    glob_basic_grid = BasicGrid(lon.flatten() / f_lon, lat.flatten()  / f_lat)
     glob_cell_grid = glob_basic_grid.to_cell_grid(cellsize=5.)
 
     return glob_cell_grid
 
-
 def ERA_IrregularImgGrid(lons, lats):
     lons_gt_180 = np.where(lons > 180.0)
     lons[lons_gt_180] = lons[lons_gt_180] - 360
-    return BasicGrid(lons.flatten(), lats.flatten()).to_cell_grid(cellsize=5.)
-
+    return BasicGrid(lons.flatten(), lats.flatten()).to_cell_grid(cellsize=5.)

ecmwf_models/interface.py

@@ -192,8 +192,8 @@ class ERANcDs(MultiTemporalImageBase):
     array_1D: bool, optional (default: False)
         Read data as list, instead of 2D array, used for reshuffling.
     """
-    def __init__(self, root_path, product, parameter=['swvl1', 'swvl2'],
-                 subgrid=None, mask_seapoints=False, h_steps=[0, 6, 12, 18],
+    def __init__(self, root_path, product, parameter=('swvl1', 'swvl2'),
+                 subgrid=None, mask_seapoints=False, h_steps=(0, 6, 12, 18),
                  array_1D=False):
 
         self.h_steps = h_steps
@@ -266,7 +266,7 @@ class ERAGrbImg(ImageBase):
         Mode in which to open the file, changing this can cause data loss.
         This argument should not be changed!
     """
-    def __init__(self, filename, product, parameter=['swvl1', 'swvl2'],
+    def __init__(self, filename, product, parameter=('swvl1', 'swvl2'),
                  subgrid=None, mask_seapoints=False, array_1D=True, mode='r'):
 
         super(ERAGrbImg, self).__init__(filename, mode=mode)
@@ -276,43 +276,58 @@ def __init__(self, filename, product, parameter=['swvl1', 'swvl2'],
 
         self.parameter = lookup(product, parameter)[
             'short_name'].values  # look up short names
+        self.product = product
 
         self.mask_seapoints = mask_seapoints
         self.array_1D = array_1D
         self.subgrid = subgrid
 
     def read(self, timestamp=None):
-        '''
+        """
         Read data from the loaded image file.
 
         Parameters
         ---------
         timestamp : datetime, optional (default: None)
             Specific date (time) to read the data for.
-        '''
+        """
         grbs = pygrib.open(self.filename)
 
         grid = self.subgrid
 
         return_img = {}
         return_metadata = {}
 
-        sea_mask = None
 
-        for message in grbs:
-            param_name = message.short_name
+        var_msg_lut = {p : None for p in self.parameter}
+        sea_mask = None
+        for N in range(grbs.messages):
+            n = N+1
+            message = grbs.message(n)
+            param_name = str(message.cfVarNameECMF)
 
             if param_name == 'lsm':
                 if self.mask_seapoints and sea_mask is None:
                     sea_mask = message.values.flatten()
 
-            param_name = message.short_name
             if param_name not in self.parameter:
                 continue
+            else:
+                var_msg_lut[param_name] = n
+
+        # available variables
+        shape = None
+        for param_name, n in var_msg_lut.items():
+            if n is None:
+                continue
 
             return_metadata[param_name] = {}
-            param_data = message.values.flatten()
 
+            message = grbs.message(n)
+
+            param_data = message.values.flatten()
+            if not shape:
+                shape  = param_data.shape
             return_img[param_name] = param_data
 
             if grid is None:
@@ -347,10 +362,21 @@ def read(self, timestamp=None):
 
         grbs.close()
 
+        # missing variables
+        for param_name, n in var_msg_lut.items():
+            if n is not None:
+                continue
+            param_data= np.full(shape, np.nan)
+            warnings.warn('Cannot load variable {var} from file {thefile}. '
+                          'Filling image with NaNs.'.format(var=param_name, thefile=self.filename))
+            return_img[param_name] = param_data
+            return_metadata[param_name] = {}
+            return_metadata[param_name]['long_name'] = \
+                lookup(self.product, [param_name]).iloc[0]['long_name']
+
         if self.array_1D:
             return Image(grid.activearrlon, grid.activearrlat,
                          return_img, return_metadata, timestamp)
-
         else:
             nlat = np.unique(grid.activearrlat).size
             nlon = np.unique(grid.activearrlon).size
@@ -389,8 +415,8 @@ class ERAGrbDs(MultiTemporalImageBase):
     expand_grid: bool, optional (default: True)
         If the reduced gaussian grid should be expanded to a full gaussian grid.
     """
-    def __init__(self, root_path, product, parameter=['swvl1', 'swvl2'],
-                 subgrid=None, mask_seapoints=False, h_steps=[0, 6, 12, 18],
+    def __init__(self, root_path, product, parameter=('swvl1', 'swvl2'),
+                 subgrid=None, mask_seapoints=False, h_steps=(0, 6, 12, 18),
                  array_1D=True):
 
         self.h_steps = h_steps
@@ -481,4 +507,4 @@ def __init__(self, ts_path, grid_path=None, **kwargs):
             grid_path = os.path.join(ts_path, "grid.nc")
 
         grid = load_grid(grid_path)
-        super(ERATs, self).__init__(ts_path, grid, **kwargs)
+        super(ERATs, self).__init__(ts_path, grid, **kwargs)

requirements.txt

@@ -12,6 +12,7 @@ pygeobase
 pygeogrids
 repurpose
 pynetcf
+pygrib
 configobj
 requests
 more_itertools