|
| 1 | +"""ESMValTool CMORizer for GLWD data. |
| 2 | +
|
| 3 | +Tier |
| 4 | + Tier 2: other freely-available dataset. |
| 5 | +
|
| 6 | +Source |
| 7 | + https://figshare.com/articles/dataset/Global_Lakes_and_Wetlands_Database_GLWD_version_2_0/28519994 |
| 8 | +
|
| 9 | +Last access |
| 10 | + 20250701 |
| 11 | +
|
| 12 | +Download and processing instructions |
| 13 | + Download the file GLWD_v2_0_combined_classes_tif.zip |
| 14 | +
|
| 15 | +""" |
| 16 | + |
| 17 | +import logging |
| 18 | +import shutil |
| 19 | +import zipfile |
| 20 | +from datetime import datetime |
| 21 | +from pathlib import Path |
| 22 | + |
| 23 | +import iris |
| 24 | +import numpy as np |
| 25 | +from cf_units import Unit |
| 26 | +from iris.coords import AuxCoord, CellMethod, DimCoord |
| 27 | +from osgeo import gdal |
| 28 | + |
| 29 | +from esmvaltool.cmorizers.data import utilities as utils |
| 30 | + |
| 31 | +logger = logging.getLogger(__name__) |
| 32 | + |
| 33 | + |
| 34 | +TIME_UNITS = Unit("days since 1950-01-01 00:00:00", calendar="standard") |
| 35 | + |
| 36 | + |
| 37 | +def _create_time_coord(): |
| 38 | + """Create time coordinate.""" |
| 39 | + # Time bounds of the climatology are set to 1984-2020 following the |
| 40 | + # corresponding publication: https://doi.org/10.5194/essd-17-2277-2025 |
| 41 | + time_points = TIME_UNITS.date2num([datetime(2002, 7, 2)]) |
| 42 | + time_bounds = [datetime(1984, 1, 1), datetime(2020, 12, 31)] |
| 43 | + time_bounds = TIME_UNITS.date2num(time_bounds) |
| 44 | + # Add new time coordinate to cube |
| 45 | + return DimCoord( |
| 46 | + time_points, |
| 47 | + bounds=time_bounds, |
| 48 | + standard_name="time", |
| 49 | + long_name="time", |
| 50 | + var_name="time", |
| 51 | + units=TIME_UNITS, |
| 52 | + climatological=True, |
| 53 | + ) |
| 54 | + |
| 55 | + |
| 56 | +def _get_bounds(points, resolution): |
| 57 | + """Compute bounds following points and resolution.""" |
| 58 | + lower = points - resolution / 2 |
| 59 | + upper = points + resolution / 2 |
| 60 | + return np.stack([lower, upper], axis=1) |
| 61 | + |
| 62 | + |
| 63 | +def _create_lat_lon_coords(n_lat, n_lon): |
| 64 | + """Create latitude/longitude coordinates.""" |
| 65 | + # The product is covering the area from 180° West to 180° East |
| 66 | + # and from 56° South to 84° North with a resolution of 15 arc-second. |
| 67 | + lat_start = -56 |
| 68 | + lon_start = -180 |
| 69 | + res = 1 / 240 # 15 / 3600 |
| 70 | + # Create coordinate points |
| 71 | + lat_points = lat_start + res * (0.5 + np.arange(n_lat)) |
| 72 | + lon_points = lon_start + res * (0.5 + np.arange(n_lon)) |
| 73 | + # Define bounds for lat/lon coordinates |
| 74 | + lat_bounds = _get_bounds(lat_points, res) |
| 75 | + lon_bounds = _get_bounds(lon_points, res) |
| 76 | + # Define coordinates |
| 77 | + latitude = DimCoord( |
| 78 | + lat_points, |
| 79 | + bounds=lat_bounds, |
| 80 | + standard_name="latitude", |
| 81 | + var_name="lat", |
| 82 | + long_name="Longitude", |
| 83 | + units="degrees", |
| 84 | + ) |
| 85 | + longitude = DimCoord( |
| 86 | + lon_points, |
| 87 | + bounds=lon_bounds, |
| 88 | + standard_name="longitude", |
| 89 | + var_name="lon", |
| 90 | + long_name="Longitude", |
| 91 | + units="degrees", |
| 92 | + ) |
| 93 | + return latitude, longitude |
| 94 | + |
| 95 | + |
| 96 | +def _create_typewetla_coord(): |
| 97 | + """Create wetland type coordinate.""" |
| 98 | + typewetla = AuxCoord( |
| 99 | + "wetland", |
| 100 | + var_name="typewetla", |
| 101 | + standard_name="area_type", |
| 102 | + long_name="Wetland", |
| 103 | + units=Unit("no unit"), |
| 104 | + ) |
| 105 | + return typewetla |
| 106 | + |
| 107 | + |
| 108 | +def _extract_variable(var, var_info, cmor_info, attrs, filedir, out_dir, cfg): |
| 109 | + """Extract variable.""" |
| 110 | + logger.info("Loading input files...") |
| 111 | + |
| 112 | + # Load data of wetland area |
| 113 | + ds = gdal.Open(Path(filedir) / cfg["area_file"]) |
| 114 | + array = ds.ReadAsArray() |
| 115 | + n_lat, n_lon = array.shape |
| 116 | + |
| 117 | + # Get ocean/fill_value mask from main class array |
| 118 | + # Classes in [0=dry-land, 1,..., 33], fill_value(ocean) = 255 |
| 119 | + dl = gdal.Open(Path(filedir) / cfg["main_class_file"]) |
| 120 | + main_class = dl.ReadAsArray() |
| 121 | + mask = main_class == 255 |
| 122 | + |
| 123 | + logger.info("Fixing data and creating coordinates...") |
| 124 | + |
| 125 | + # Fix data: |
| 126 | + # - mask oceans (fill_value = 255) + set value to 0 |
| 127 | + # - flip latitude axis |
| 128 | + array = np.where(array <= 100, array, 0) |
| 129 | + array = np.ma.array(array, mask=mask) |
| 130 | + array = np.flip(array, axis=0) |
| 131 | + |
| 132 | + # Time coordinate |
| 133 | + time_coord = _create_time_coord() |
| 134 | + |
| 135 | + # Latitude and longitude coordinates |
| 136 | + latitude_coord, longitude_coord = _create_lat_lon_coords(n_lat, n_lon) |
| 137 | + |
| 138 | + # Type wetland coordinate |
| 139 | + typewetla_coord = _create_typewetla_coord() |
| 140 | + |
| 141 | + # Cube data |
| 142 | + logger.info("Setting up the cube for variable %s", var) |
| 143 | + cube = iris.cube.Cube( |
| 144 | + array, |
| 145 | + standard_name=cmor_info.standard_name, |
| 146 | + units=cmor_info.units, |
| 147 | + dim_coords_and_dims=[(latitude_coord, 0), (longitude_coord, 1)], |
| 148 | + ) |
| 149 | + cube.add_aux_coord(time_coord, ()) |
| 150 | + cube.add_aux_coord(typewetla_coord, ()) |
| 151 | + |
| 152 | + # Add coordinate time axis of size 1 |
| 153 | + cube = iris.util.new_axis(cube, "time") |
| 154 | + |
| 155 | + # Fix cell methods |
| 156 | + cube.add_cell_method(CellMethod("mean within years", coords=time_coord)) |
| 157 | + cube.add_cell_method(CellMethod("mean over years", coords=time_coord)) |
| 158 | + |
| 159 | + # Fix coords |
| 160 | + cube = utils.fix_coords(cube) |
| 161 | + |
| 162 | + # Fix var metadata |
| 163 | + utils.fix_var_metadata(cube, cmor_info) |
| 164 | + |
| 165 | + # Fix global metadata |
| 166 | + utils.set_global_atts(cube, attrs) |
| 167 | + |
| 168 | + # Save variable |
| 169 | + utils.save_variable(cube, var, out_dir, attrs) |
| 170 | + |
| 171 | + |
| 172 | +def _unzip(filepath, out_dir): |
| 173 | + """Unzip `*.zip` file.""" |
| 174 | + extracted_dir = Path(out_dir) / "tmp_extracted_files" |
| 175 | + logger.info("Starting extraction of %s to %s", filepath, extracted_dir) |
| 176 | + with zipfile.ZipFile(filepath, "r") as zip_ref: |
| 177 | + zip_ref.extractall(extracted_dir) |
| 178 | + logger.info("Succefully extracted file to %s", extracted_dir) |
| 179 | + return extracted_dir |
| 180 | + |
| 181 | + |
| 182 | +def cmorization(in_dir, out_dir, cfg, cfg_user, start_date, end_date): |
| 183 | + """Cmorization func call.""" |
| 184 | + cmor_table = cfg["cmor_table"] |
| 185 | + glob_attrs = cfg["attributes"] |
| 186 | + |
| 187 | + # Run the cmorization |
| 188 | + for var, var_info in cfg["variables"].items(): |
| 189 | + logger.info("CMORizing variable '%s'", var) |
| 190 | + glob_attrs["mip"] = var_info["mip"] |
| 191 | + cmor_info = cmor_table.get_variable(var_info["mip"], var) |
| 192 | + |
| 193 | + # Extract file from ZIP archive |
| 194 | + zip_file = Path(in_dir) / cfg["archive_filename"] |
| 195 | + if not Path(zip_file).is_file(): |
| 196 | + logger.debug("Skipping '%s', file '%s' not found", var, zip_file) |
| 197 | + continue |
| 198 | + logger.info("Found input file '%s'", zip_file) |
| 199 | + filedir = _unzip(zip_file, out_dir) |
| 200 | + |
| 201 | + # Extract and save variable file |
| 202 | + _extract_variable( |
| 203 | + var, var_info, cmor_info, glob_attrs, filedir, out_dir, cfg |
| 204 | + ) |
| 205 | + |
| 206 | + # Remove extracted directory |
| 207 | + shutil.rmtree(Path(filedir)) |
| 208 | + logger.info("Removed cached input directory %s", filedir) |
0 commit comments