Fix problematic ancillary pixels

src/common_utilities.py

@@ -0,0 +1,110 @@
+# Copyright 2024 ACCESS-NRI (https://www.access-nri.org.au/)
+# See the top-level COPYRIGHT.txt file for details.
+#
+# SPDX-License-Identifier: Apache-2.0
+#
+# Created by: Chermelle Engel <Chermelle.Engel@anu.edu.au>
+
+import mule
+import iris
+import xarray as xr
+import numpy as np
+
+class ReplaceOperator(mule.DataOperator):
+    """ Mule operator for replacing the data"""
+    def __init__(self):
+        pass
+    def new_field(self, sources):
+        print('new_field')
+        return sources[0]
+    def transform(self, sources, result):
+        print('transform')
+        return sources[1]
+
+def replace_in_ff_problematic(f, mf_out, replace, stashcode, canopy_pixels, landsea_pixels):
+
+    current_data = f.get_data()
+    data=current_data.copy()
+
+    if stashcode == 218:
+        for j in range(len(canopy_pixels)):
+            iy=canopy_pixels[j][0]
+            ix=canopy_pixels[j][1]
+            if np.isnan(data[iy,ix]):
+                data[iy,ix]=1.
+    elif stashcode == 33:
+        for j in range(len(landsea_pixels)):
+            data[landsea_pixels[j][0],landsea_pixels[j][1]]=0.
+
+    mf_out.fields.append(replace([f, data]))
+
+def problematic_pixels(infile):
+    """
+    Function to get locations of pixels that have incomplete canopy height information or 
+    problematic land/sea definitions.
+
+    Parameters
+    ----------
+    infile : string
+        The name of the start dump file to read
+
+    Returns
+    -------
+    2d numpy array
+        A 2-D numpy array containg the field data for the date/time and spatial extent
+    """
+
+    # defining the variable names of surface altitude, soil porosity and canopy height in the start dump files
+    var_altitude="surface_altitude"
+    var_soil="soil_porosity"
+    var_canopy="canopy_height"
+
+    # Read in the relevant data
+    # One single read of the relevant variables is faster than individual reads
+    d = iris.load(infile,[var_altitude,var_soil,var_canopy])
+
+    # Creating an orography mask and storing the latitude and longitude data for reporting purposes
+    orog_data=xr.DataArray.from_iris(d[0])
+    lats=orog_data['latitude'].data
+    lons=orog_data['longitude'].data
+    orog_mask=np.where(orog_data.data[:,:]>0.,1,0) 
+
+    # Creating a soil parameter based mask
+    soil_data=xr.DataArray.from_iris(d[1])
+    soil_mask=np.where(np.isnan(soil_data.data[:,:]),1,0)
+
+    # Creating a canopy mask - screening across all five different plant types
+    canopy_data=xr.DataArray.from_iris(d[2])
+    canopy_mask=np.where(np.isnan(canopy_data.data[:,:,:]),1,0)
+    canopy_mask=np.sum(canopy_mask,axis=0)
+    canopy_mask=np.where(canopy_mask>0,1,0)
+
+    # Creating a compound mask that will indicate 1 for just canopy nan or 2 for canopy and soil nan
+    compound_data=canopy_mask+soil_mask
+
+    # Removing any sea points (points with altitude 0 or less
+    compound_data=compound_data*orog_mask
+
+    # Finding locations of problematic pixels in two sets
+    # One for missing canopy alone and one for misclassified land
+    try:
+        canopy_pixels=np.argwhere(compound_data==1).compute()
+        landsea_pixels=np.argwhere(compound_data==2).compute()
+    except:
+        canopy_pixels=np.argwhere(compound_data==1)
+        landsea_pixels=np.argwhere(compound_data==1)
+
+    # Printing information to the standard output for reporting purposes (so scientists can be aware)
+    npoints,nxy = (canopy_pixels.shape)
+    if npoints>0:
+        for i in range(npoints):
+            print("%.1f, %.1f, Nan canopy"%(lons[canopy_pixels[i,1]],lats[canopy_pixels[i,0]]))
+
+    npoints,nxy = (landsea_pixels.shape)
+    if npoints>0:
+        for i in range(npoints):
+            print("%.1f, %.1f, Misplaced Orography"%(lons[landsea_pixels[i,1]],lats[landsea_pixels[i,0]]))
+
+    # returning pixel locations so they can be processed appropriately
+    return canopy_pixels,landsea_pixels
+

src/hres_eccb.py

@@ -1,11 +1,11 @@
+#!/usr/bin/env python3
 # Copyright 2024 ACCESS-NRI (https://www.access-nri.org.au/)
 # See the top-level COPYRIGHT.txt file for details.
 #
 # SPDX-License-Identifier: Apache-2.0
 #
 # Created by: Chermelle Engel <Chermelle.Engel@anu.edu.au>
 
-#!/usr/bin/env python3
 
 """
 Replace the land/surface fields in the ec_cb000 file with higher-resolution

src/hres_ic.py

@@ -1,11 +1,11 @@
+#!/usr/bin/env python3
 # Copyright 2024 ACCESS-NRI (https://www.access-nri.org.au/)
 # See the top-level COPYRIGHT.txt file for details.
 #
 # SPDX-License-Identifier: Apache-2.0
 #
 # Created by: Chermelle Engel <Chermelle.Engel@anu.edu.au>
 
-#!/usr/bin/env python3
 
 """
 Replace the land/surface fields in the astart file with higher-resolution
@@ -59,13 +59,13 @@ def main():
 
     # If necessary replace ERA5 land/surface fields with higher-resolution options
     if "era5land" in args.type:
-        replace_landsurface_with_ERA5land_IC.swap_land_era5land(args.mask, args.file, t)
+        replace_landsurface_with_ERA5land_IC.swap_land_era5land(args.mask, args.file, t, fix_problematic_pixels="yes")
         shutil.move(args.file.as_posix(), args.file.as_posix().replace('.tmp', ''))
     elif "barra" in args.type:
-        replace_landsurface_with_BARRA2R_IC.swap_land_barra(args.mask, args.file, t)
+        replace_landsurface_with_BARRA2R_IC.swap_land_barra(args.mask, args.file, t, fix_problematic_pixels="yes")
         shutil.move(args.file.as_posix(), args.file.as_posix().replace('.tmp', ''))
     elif "astart" in args.type:
-        replace_landsurface_with_FF_IC.swap_land_ff(args.mask, args.file, args.hres_ic,t)
+        replace_landsurface_with_FF_IC.swap_land_ff(args.mask, args.file, args.hres_ic,t, fix_problematic_pixels="yes")
         shutil.move(args.file.as_posix(), args.file.as_posix().replace('.tmp', ''))
 
     else:

src/replace_landsurface_with_BARRA2R_IC.py

@@ -20,18 +20,6 @@
 BARRA_DIR = os.path.join(ROSE_DATA, 'etc', 'barra_r2')
 
 
-class ReplaceOperator(mule.DataOperator):
-    """ Mule operator for replacing the data"""
-    def __init__(self):
-        pass
-    def new_field(self, sources):
-        print('new_field')
-        return sources[0]
-    def transform(self, sources, result):
-        print('transform')
-        return sources[1]
-
-
 class bounding_box(): 
     """ Container class to hold spatial extent information."""
     def __init__(self, ncfname, maskfname, var):
@@ -156,7 +144,7 @@ def get_BARRA_nc_data(ncfname, FIELDN, wanted_dt, NLAYERS, bounds):
     return data.data
 
 
-def swap_land_barra(mask_fullpath, ec_cb_file_fullpath, ic_date):
+def swap_land_barra(mask_fullpath, ec_cb_file_fullpath, ic_date, fix_problematic_pixels="no"):
     """
     Function to get the BARRA2-R data for all land/surface variables.
 
@@ -187,6 +175,9 @@ def swap_land_barra(mask_fullpath, ec_cb_file_fullpath, ic_date):
     # Path to input file
     ff_in = ec_cb_file_fullpath.as_posix().replace('.tmp', '')
 
+    if fix_problematic_pixels == "yes":
+        canopy_pixels,landsea_pixels=problematic_pixels(ff_in)
+
     # Path to output file
     ff_out = ec_cb_file_fullpath.as_posix()
     print(ff_in, ff_out)
@@ -195,7 +186,7 @@ def swap_land_barra(mask_fullpath, ec_cb_file_fullpath, ic_date):
     mf_in = mule.load_umfile(ff_in)
 
     # Create Mule Replacement Operator
-    replace = ReplaceOperator() 
+    replace = common_utilities.ReplaceOperator() 
 
     # Read in the surface temperature data from the archive
     BARRA_FIELDN = 'ts'
@@ -232,7 +223,6 @@ def swap_land_barra(mask_fullpath, ec_cb_file_fullpath, ic_date):
     # For each field in the input write to the output file (but modify as required)
     for f in mf_in.fields:
 
-      print(f.lbuser4, f.lblev, f.lblrec, f.lbhr, f.lbcode)
       if f.lbuser4 == 9:
         # replace coarse soil moisture with high-res information
         current_data = f.get_data()
@@ -251,6 +241,9 @@ def swap_land_barra(mask_fullpath, ec_cb_file_fullpath, ic_date):
         data = surface_temp
         data = np.where(np.isnan(data), current_data, data)
         mf_out.fields.append(replace([f, data]))
+      elif ((f.lbuser4 == 33) or (f.lbuser4 == 218)) and fix_problematic_pixels == "yes":
+        # replace surface altitude and canopy_height
+        common_utilities.replace_in_ff_problematic(f, mf_out, replace,f.lbuser4,canopy_pixels,landsea_pixels)
       else:
         mf_out.fields.append(f)
 

src/replace_landsurface_with_ERA5land_IC.py

@@ -14,6 +14,8 @@
 from pathlib import Path
 import xarray as xr, sys, argparse
 from datetime import datetime,timedelta
+import common_utilities
+
 
 ROSE_DATA = os.environ.get('ROSE_DATA')
 # Base directory of the ERA5-land archive on NCI
@@ -23,17 +25,6 @@
 ##########multipliers=[7.*10., 21.*10., 72.*10., 189.*10.]
 multipliers = [10.*10., 25.*10., 65.*10., 200.*10.]
 
-class ReplaceOperator(mule.DataOperator):
-    """ Mule operator for replacing the data"""
-    def __init__(self):
-        pass
-    def new_field(self, sources):
-        print('new_field')
-        return sources[0]
-    def transform(self, sources, result):
-        print('transform')
-        return sources[1]
-
 class bounding_box():
     """ Container class to hold spatial extent information."""
     def __init__(self, ncfname, maskfname, var):
@@ -196,7 +187,7 @@ def replace_in_ff(f, generic_era5_fname, ERA_FIELDN, multiplier, ic_z_date, mf_o
     data = np.where(np.isnan(data), current_data, data)
     mf_out.fields.append(replace([f, data]))
 
-def swap_land_era5land(mask_fullpath, ic_file_fullpath, ic_date):
+def swap_land_era5land(mask_fullpath, ic_file_fullpath, ic_date, fix_problematic_pixels="no"):
     """
     Function to get the ERA5-land data for all land/surface variables.
 
@@ -236,6 +227,9 @@ def swap_land_era5land(mask_fullpath, ic_file_fullpath, ic_date):
     # Path to input file 
     ff_in = ic_file_fullpath.as_posix().replace('.tmp', '')
 
+    if fix_problematic_pixels == "yes":
+        canopy_pixels,landsea_pixels=common_utilities.problematic_pixels(ff_in)
+
     # Path to output file 
     ff_out = ic_file_fullpath.as_posix()
     print(ff_in, ff_out)
@@ -244,7 +238,7 @@ def swap_land_era5land(mask_fullpath, ic_file_fullpath, ic_date):
     mf_in = mule.load_umfile(ff_in)
 
     # Create Mule Replacement Operator
-    replace = ReplaceOperator() 
+    replace = common_utilities.ReplaceOperator() 
 
     # Define spatial extent of grid required
     bounds = bounding_box(era5_fname, mask_fullpath.as_posix(), "land_binary_mask")
@@ -255,8 +249,6 @@ def swap_land_era5land(mask_fullpath, ic_file_fullpath, ic_date):
     # For each field in the input write to the output file (but modify as required)
     for f in mf_in.fields:
 
-      print(f.lbuser4, f.lblev, f.lblrec, f.lbhr, f.lbcode)
-
       if f.lbuser4 == 9:
         # replace coarse soil moisture with high-res information
         if f.lblev == 4:
@@ -280,7 +272,13 @@ def swap_land_era5land(mask_fullpath, ic_file_fullpath, ic_date):
           replace_in_ff(f, generic_era5_fname, 'stl1', -1, ic_z_date, mf_out, replace, bounds)
 
       elif f.lbuser4 == 24:
+        # surface temperature
         replace_in_ff(f, generic_era5_fname, 'skt', -1, ic_z_date, mf_out, replace, bounds)
+
+      elif ((f.lbuser4 == 33) or (f.lbuser4 == 218)) and fix_problematic_pixels == "yes":
+        # surface altitude and canopy_height
+        common_utilities.replace_in_ff_problematic(f, mf_out, replace,f.lbuser4,canopy_pixels,landsea_pixels)
+
       else:
         mf_out.fields.append(f)
 

src/replace_landsurface_with_FF_IC.py

@@ -7,23 +7,12 @@
 
 import mule
 
-class ReplaceOperator(mule.DataOperator):
-    """ Mule operator for replacing the data"""
-    def __init__(self):
-        pass
-    def new_field(self, sources):
-        print('new_field')
-        return sources[0]
-    def transform(self, sources, result):
-        print('transform')
-        return sources[1]
-
 def replace_in_ff_from_ff(f, sf, mf_out, replace):
 
     replacement_data = sf.get_data()
     mf_out.fields.append(replace([f, replacement_data]))
 
-def swap_land_ff(mask_fullpath, ic_file_fullpath, source_fullpath, ic_date):
+def swap_land_ff(mask_fullpath, ic_file_fullpath, source_fullpath, ic_date, fix_problematic_pixels="no"):
     """
     Function to get the land/surface data from another fields file into the start dump.
 
@@ -49,6 +38,9 @@ def swap_land_ff(mask_fullpath, ic_file_fullpath, source_fullpath, ic_date):
     # Path to input file 
     ff_in = ic_file_fullpath.as_posix().replace('.tmp', '')
 
+    if fix_problematic_pixels == "yes":
+        canopy_pixels,landsea_pixels=common_utilities.problematic_pixels(ff_in)
+
     # Path to input file 
     sf_in = source_fullpath.as_posix()
 
@@ -61,7 +53,7 @@ def swap_land_ff(mask_fullpath, ic_file_fullpath, source_fullpath, ic_date):
     msf_in = mule.load_umfile(sf_in)
 
     # Create Mule Replacement Operator
-    replace = ReplaceOperator() 
+    replace = common_utilities.ReplaceOperator() 
 
     # Set up the output file
     mf_out = mf_in.copy()
@@ -71,6 +63,9 @@ def swap_land_ff(mask_fullpath, ic_file_fullpath, source_fullpath, ic_date):
 
         if f.lbuser4 in [9, 20, 24]:
             replace_in_ff_from_ff(f, sf, mf_out, replace)
+        elif ((f.lbuser4 == 33) or (f.lbuser4 == 218)) and fix_problematic_pixels == "yes":
+            # surface altitude and canopy_height
+            common_utilities.replace_in_ff_problematic(f, mf_out, replace,f.lbuser4,canopy_pixels,landsea_pixels)
         else:
             mf_out.fields.append(f)