diff --git a/src/common_utilities.py b/src/common_utilities.py new file mode 100644 index 0000000..57b0b89 --- /dev/null +++ b/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 + +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 + diff --git a/src/hres_eccb.py b/src/hres_eccb.py index ddc8646..cd5b6b5 100755 --- a/src/hres_eccb.py +++ b/src/hres_eccb.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python3 # Copyright 2024 ACCESS-NRI (https://www.access-nri.org.au/) # See the top-level COPYRIGHT.txt file for details. # @@ -5,7 +6,6 @@ # # Created by: Chermelle Engel -#!/usr/bin/env python3 """ Replace the land/surface fields in the ec_cb000 file with higher-resolution diff --git a/src/hres_ic.py b/src/hres_ic.py index 09fdae0..60c2cfb 100755 --- a/src/hres_ic.py +++ b/src/hres_ic.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python3 # Copyright 2024 ACCESS-NRI (https://www.access-nri.org.au/) # See the top-level COPYRIGHT.txt file for details. # @@ -5,7 +6,6 @@ # # Created by: Chermelle Engel -#!/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: diff --git a/src/replace_landsurface_with_BARRA2R_IC.py b/src/replace_landsurface_with_BARRA2R_IC.py index e041382..d93be25 100755 --- a/src/replace_landsurface_with_BARRA2R_IC.py +++ b/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) diff --git a/src/replace_landsurface_with_ERA5land_IC.py b/src/replace_landsurface_with_ERA5land_IC.py index 46fdd71..c3f28c8 100755 --- a/src/replace_landsurface_with_ERA5land_IC.py +++ b/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) diff --git a/src/replace_landsurface_with_FF_IC.py b/src/replace_landsurface_with_FF_IC.py index a5b3bc5..10508a8 100755 --- a/src/replace_landsurface_with_FF_IC.py +++ b/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)