|
| 1 | +#!/usr/bin/env python3 |
| 2 | + |
| 3 | +from netCDF4 import Dataset |
| 4 | +from scipy.io import FortranFile |
| 5 | +import numpy as np |
| 6 | +import sys |
| 7 | +import os |
| 8 | +import glob |
| 9 | + |
| 10 | +def get_i_j_pf(til_file): |
| 11 | + ii = [] |
| 12 | + jj = [] |
| 13 | + pfaf = [] |
| 14 | + with open(til_file, 'r') as f: |
| 15 | + for line in f: |
| 16 | + data = line.split() |
| 17 | + if data[0] == '100' : |
| 18 | + ii.append(int(data[4])) |
| 19 | + jj.append(int(data[5])) |
| 20 | + pfaf.append(int(data[7])) |
| 21 | + return ii, jj, pfaf |
| 22 | + |
| 23 | +def create_reduced_nc4(input_file, output_file, indices): |
| 24 | + """ |
| 25 | + Create a new nc4 file with with indices of for the input file. |
| 26 | + |
| 27 | + Parameters: |
| 28 | + input_file (str): Path to the input nc4 file |
| 29 | + output_file (str): Path to the output nc4 file |
| 30 | + indices(int array): |
| 31 | + """ |
| 32 | + |
| 33 | + # Open the input file in read mode |
| 34 | + try: |
| 35 | + src = Dataset(input_file, 'r') |
| 36 | + print(f"Successfully opened input file: {input_file}") |
| 37 | + except Exception as e: |
| 38 | + print(f"Error opening input file: {e}") |
| 39 | + return |
| 40 | + |
| 41 | + # Create the output file in write mode |
| 42 | + try: |
| 43 | + dst = Dataset(output_file, 'w', format='NETCDF4') |
| 44 | + print(f"Creating output file: {output_file}") |
| 45 | + except Exception as e: |
| 46 | + print(f"Error creating output file: {e}") |
| 47 | + src.close() |
| 48 | + return |
| 49 | + |
| 50 | + # Copy global attributes |
| 51 | + dst.setncatts(src.__dict__) |
| 52 | + |
| 53 | + # Copy dimensions with reduced size for unlimited/time dimensions |
| 54 | + for name, dimension in src.dimensions.items(): |
| 55 | + if name == 'tile': |
| 56 | + dst.createDimension(name, len(indices)) |
| 57 | + else: |
| 58 | + dst.createDimension(name, len(dimension)) |
| 59 | + # Copy variables with reduced data |
| 60 | + for name, variable in src.variables.items(): |
| 61 | + # Create variable with same attributes but possibly reduced dimensions |
| 62 | + dims = variable.dimensions |
| 63 | + new_dims = [] |
| 64 | + for dim in dims: |
| 65 | + new_dims.append(dim) |
| 66 | + # Create the variable in the output file |
| 67 | + dst_var = dst.createVariable( |
| 68 | + name, |
| 69 | + variable.datatype, |
| 70 | + new_dims, |
| 71 | + ) |
| 72 | + |
| 73 | + # Copy variable attributes |
| 74 | + dst[name].setncatts(variable.__dict__) |
| 75 | + |
| 76 | + # Get the data and reduce it |
| 77 | + if len(dims) == 2 : |
| 78 | + dst_var[:,:] = variable[:,indices] |
| 79 | + if len(dims) == 1 : |
| 80 | + dst_var[:] = variable[indices] |
| 81 | + |
| 82 | + |
| 83 | + # Close both files |
| 84 | + src.close() |
| 85 | + dst.close() |
| 86 | + print(f"Successfully created reduced file: {output_file}") |
| 87 | + |
| 88 | +def create_reduced_bin(input_file, output_file, indices): |
| 89 | + """ |
| 90 | + Create a new ibinary file with with indices of for the input file. |
| 91 | + |
| 92 | + Parameters: |
| 93 | + input_file (str): Path to the input binary file |
| 94 | + output_file (str): Path to the output binaryfile |
| 95 | + indices(int array): |
| 96 | + """ |
| 97 | + nland = len(indices) |
| 98 | + fout = FortranFile(output_file, 'w') |
| 99 | + with FortranFile(input_file, 'r') as f: |
| 100 | + while True : |
| 101 | + try: |
| 102 | + a = f.read_reals(dtype=np.float32) |
| 103 | + b = f.read_reals(dtype=np.float32) |
| 104 | + a[12] = np.float32(nland) |
| 105 | + fout.write_record(a) |
| 106 | + b = b[indices] |
| 107 | + fout.write_record(b) |
| 108 | + except : |
| 109 | + break |
| 110 | + fout.close() |
| 111 | + print(f"Successfully created reduced file: {output_file}") |
| 112 | + |
| 113 | +# Example usage |
| 114 | +if __name__ == "__main__": |
| 115 | + # Replace these with your actual file paths |
| 116 | + # input_nc4 = "/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles/v13/land/CF1440x6C/clsm/catchcn_params.nc4" |
| 117 | + # output_nc4 = "output_file_reduced.nc" |
| 118 | + bcs_dir = sys.argv[1] |
| 119 | + bcs_ver = sys.argv[2] |
| 120 | + air_res = sys.argv[3] |
| 121 | + ocn_res = sys.argv[4] |
| 122 | + agname = air_res + '_' + air_res |
| 123 | + orig_tile = bcs_dir + '/' + bcs_ver + '/geometry/' + agname + '/' + agname +'-Pfafstetter.til' |
| 124 | + if not os.path.exists(orig_tile) : |
| 125 | + print( "The original tile file must exist " + orig_tile) |
| 126 | + exit() |
| 127 | + |
| 128 | + mom_tile = 'til/' + air_res + '_' + ocn_res + '-Pfafstetter.til' |
| 129 | + if not os.path.exists(mom_tile) : |
| 130 | + print( "The MOM tile file must exist " + mom_tile) |
| 131 | + exit() |
| 132 | + |
| 133 | + ii1, jj1, pf1 = get_i_j_pf(orig_tile) |
| 134 | + ii2, jj2, pf2 = get_i_j_pf(mom_tile) |
| 135 | + i1 = 0 |
| 136 | + indices =[] |
| 137 | + for i2 in range(len(ii2)): |
| 138 | + match = ii1[i1] == ii2[i2] and jj1[i1] == jj2[i2] and pf1[i1] == pf2[i2] |
| 139 | + while not match : |
| 140 | + i1 = i1 +1 |
| 141 | + match = ii1[i1] == ii2[i2] and jj1[i1] == jj2[i2] and pf1[i1] == pf2[i2] |
| 142 | + indices.append(i1) |
| 143 | + i1 = i1 + 1 |
| 144 | + |
| 145 | + catch_params_file = bcs_dir + '/' + bcs_ver + '/land/' + air_res + '/clsm/catch_params.nc4' |
| 146 | + catchcn_params_file = bcs_dir + '/' + bcs_ver + '/land/' + air_res + '/clsm/catchcn_params.nc4' |
| 147 | + mom_catch_params_file = 'clsm/catch_params.nc4' |
| 148 | + mom_catchcn_params_file = 'clsm/catchcn_params.nc4' |
| 149 | + |
| 150 | + # create_reduced_nc4(catch_params_file, mom_catch_params_file, indices) |
| 151 | + # create_reduced_nc4(catchcn_params_file, mom_catchcn_params_file, indices) |
| 152 | + |
| 153 | + files = glob.glob(bcs_dir + '/' + bcs_ver + '/land/' + air_res + '/*.da*') |
| 154 | + for file in files: |
| 155 | + fname = os.path.basename(file) |
| 156 | + if 'vegdyn' in fname: |
| 157 | + create_reduced_nc4( file,'clsm/vegdyn.dat', indices) |
| 158 | + else: |
| 159 | + shortname = 'clsm/' + fname.split('_')[0] + '.dat' |
| 160 | + create_reduced_bin(file, shortname, indices) |
0 commit comments