Method `extract_tilefrac` and `average_tile` in `calc_land.py` lack of a key argument `lev`

[rhaegar325]

Compare with old version Mopper, our current extract_tilefrac and average_tile methods in Moppy lack of a specific argument lev, which will led to miss a coordinate type:

Coordinates:
  * time     (time) datetime64[ns] 16kB 1850-01-16T12:00:00 ... 2014-12-16T12...
  * lat      (lat) float64 1kB -90.0 -88.75 -87.5 -86.25 ... 87.5 88.75 90.0
  * lon      (lon) float64 2kB 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1
    type     |S11 11B ...

according to CMIP6_Lmon.json, there is an extra dimention for some land variables, for example:

"baresoilFrac": {
    ...
    "dimensions": "longitude latitude time typebare",
    ...
},

"c3PftFrac": {
    ...
    "dimensions": "longitude latitude time typec3pft",
    ...
}

Those different type are stored in MOPPeR/src/mopdata/landtype.yaml file in the old Mopper:

mod_mapping:
  typebare: bare_ground
  typeburnt: burnt_vegetation
  typec3pft: c3_plant_functional_types
  typec3crop: crops_of_c3_plant_functional_types
  typec3natg: natural_grasses_of_c3_plant_functional_types
  typec3pastures: pastures_of_c3_plant_functional_types
  typec4pft: c4_plant_functional_types
  typec4crop: crops_of_c4_plant_functional_types
  typec4natg: natural_grasses_of_c4_plant_functional_types
  typec4pastures: pastures_of_c4_plant_functional_types
  typecloud: cloud
  typecrop: crops
  typefis: floating_ice_shelf
  typegis: grounded_ice_sheet
  typeland: land
  typeli: land_ice
  typemp: sea_ice_melt_pond
  typenatgr: natural_grasses
  typenwd: herbaceous_vegetation
  typepasture: pastures
  typepdec: primary_deciduous_trees
  typepever: primary_evergreen_trees
  typeresidual: residual
  typesdec: secondary_decidous_trees
  typesea: sea
  typesever: secondary_evergreen_trees
  typeshrub: shrubs
  typesi: sea_ice
  typesirdg: sea_ice_ridges
  typetree: trees
  typeveg: vegetation
  typewetla: wetland

Maybe we should consider to add the machenism back to calc_land.py to make it compatible with CMIP6 convention.

[rhaegar325]

Hi, @rbeucher, while generating the mappings for CM3, I noticed that there doesn’t seem to be a method for generating the type coordinates for land fraction variables. I may have missed something, but I couldn’t find an existing implementation covering this case.

It might be worth considering adding a dedicated mechanism for generating or registering these type coordinates, as they appear to be required for the mappings to be fully consistent.

[rbeucher]

Yep I also came across that one today

[rbeucher]

Those are scalar dimensions and should be converted to a "character" coordinate.

I think we should have a more generic way of handling coordinates by using the CMIP tables.
Coordinates definitions are stored in CMIP6_coordinate.json.

For example:

        "latitude": {
            "standard_name": "latitude",
            "units": "degrees_north",
            "axis": "Y",
            "long_name": "Latitude",
            "climatology": "",
            "formula": "",
            "must_have_bounds": "yes",
            "out_name": "lat",
            "positive": "",
            "requested": "",
            "requested_bounds": "",
            "stored_direction": "increasing",
            "tolerance": "",
            "type": "double",
            "valid_max": "90.0",
            "valid_min": "-90.0",
            "value": "",
            "z_bounds_factors": "",
            "z_factors": "",
            "bounds_values": "",
            "generic_level_name": ""
        },

       "typebare": {
            "standard_name": "area_type",
            "units": "",
            "axis": "",
            "long_name": "surface type",
            "climatology": "",
            "formula": "",
            "must_have_bounds": "no",
            "out_name": "type",
            "positive": "",
            "requested": "",
            "requested_bounds": "",
            "stored_direction": "",
            "tolerance": "",
            "type": "character",
            "valid_max": "",
            "valid_min": "",
            "value": "bare_ground",
            "z_bounds_factors": "",
            "z_factors": "",
            "bounds_values": "",
            "generic_level_name": ""
        },

CMOR would go through the metadata and add non-null values to the type coordinates.
There are other similar examples, we also need to handle hybrid axis.
My point is that the extract_tile ... etc functions might not be the best place for that

[rbeucher]

Or we need a generic function for adding dimension coordinates which will use the CMIP vocabulary.

something like: add_CMIP_coordinate(dataArray, "typebare") etc...

I am still figuring out the best/cleanest approach. Open to suggestions @rhaegar325

[rhaegar325]

I agree with you @rbeucher, we don’t need to use those old functions from the legacy MOPPeR code, such as extract_tile and average_tile, to generate the extra dimensions.

However, we could add a mechanism to automatically add coordinates when the raw data lack some of the coordinates defined in the CMIP tables. Alternatively, as you suggested, we could wrap this logic into a generic function like add_CMIP_coordinate(dataArray, "typebare"), place it in calc_land, and call it during the mapping calculations.

[rhaegar325]

After multiple test, I found we actually have coorespounding process to deal with the type coords in atmosphere.py:

elif "value" in meta:
                val = meta["value"]
                # Handle character type (e.g., string coordinate)
                if meta["type"] == "character":
                    arr = xr.DataArray(
                        np.array(
                            val, dtype="S"
                        ),  # ensure type is character (byte string)
                        dims=(),
                        attrs={
                            k: v
                            for k, v in {
                                "standard_name": meta.get("standard_name"),
                                "long_name": meta.get("long_name"),
                                "units": meta.get("units"),
                                "axis": meta.get("axis"),
                                "positive": meta.get("positive"),
                                "valid_min": meta.get("valid_min"),
                                "valid_max": meta.get("valid_max"),
                            }.items()
                            if v is not None
                        },
                    )
                else:
                    arr = xr.DataArray(
                        dtype(val),
                        dims=(),
                        attrs={
                            k: v
                            for k, v in {
                                "standard_name": meta.get("standard_name"),
                                "long_name": meta.get("long_name"),
                                "units": meta.get("units"),
                                "axis": meta.get("axis"),
                                "positive": meta.get("positive"),
                                "valid_min": dtype(meta["valid_min"])
                                if "valid_min" in meta
                                else None,
                                "valid_max": dtype(meta["valid_max"])
                                if "valid_max" in meta
                                else None,
                            }.items()
                            if v is not None
                        },
                    )
                self.ds = self.ds.assign_coords({name: arr})

the problem is when we use cmoriser.write() to write the dataset, the type coords and other coords with dtype char or str will be treat as a varibale instead of a coords. I will open a new issue to describe the issue.