The previous blog showed two different ways to add coordinate information to WRF netCDF output files, one using Rich Signell's CF modifications, and the other using the CDM library's internal WRFConvention Java code. Let's look at the details of both with an eye towards advising the WRF group on how to make their files CF compliant.
1. Horizontal Coordinates
Rich's solution uses the existing 2D lat/lon fields for the horizontal coordinate reference system. There are 3 such sets: XLAT/XLONG on the center points of the grid, XLAT_U/XLON_U for the U wind component, XLAT_U/XLONG_U for the U wind component, XLAT_V/XLONG_V for the V wind component. These already have the correct CF units degree_north and degree_east, so the only thing one has to do is to reference these coordinates in the corresponding data variables. The NcML needed to do this looks like:
<variable name="U" >
<attribute name="coordinates" value="XLONG_U XLAT_U ZNU XTIME"/>
</variable>
<variable name="V" >
<attribute name="coordinates" value="XLONG_V XLAT_V ZNU XTIME"/>
</variable>
<variable name="W" >
<attribute name="coordinates" value="XLONG XLAT ZNW XTIME"/>
</variable>
...
Here we are taking 3 existing variables and adding the coordinates attribute to each. When added to the existing attributes, these variables now look like:
float U(Time=1, bottom_top=27, south_north=60, west_east_stag=74);
:FieldType = 104; // int
:MemoryOrder = "XYZ";
:description = "x-wind component";
:units = "m s-1";
:stagger = "X";
:coordinates = "XLONG_U XLAT_U ZNU XTIME";
float V(Time=1, bottom_top=27, south_north_stag=61, west_east=73);
:FieldType = 104; // int
:MemoryOrder = "XYZ";
:description = "y-wind component";
:units = "m s-1";
:stagger = "Y";
:coordinates = "XLONG_V XLAT_V ZNU XTIME";
float W(Time=1, bottom_top_stag=28, south_north=60, west_east=73);
:FieldType = 104; // int
:MemoryOrder = "XYZ";
:description = "z-wind component";
:units = "m s-1";
:stagger = "Z";
:coordinates = "XLONG XLAT ZNW XTIME";
The CF coordinates attribute simply lists the coordinates for the variable. To be fully compliant, one must add a coordinates attribute to each data variable. So the above adds the correct horizontal coordinates to these 3 variables (we will cover the Z and time coordinate below), and the rest need to be done also.
The other possibility is to add the projection and projection coordinates to the file. The projection information is in the global attributes, and the CDM WRF Convention Builder (ucar.nc2.dataset.conv.WRFConvention.java on our GitHub source site) takes this approach. This code actually uses a variation of CF called the _Coordinate Attribute Conventions:
char Lambert;
:grid_mapping_name = "lambert_conformal_conic";
:latitude_of_projection_origin = 34.83000564575195; // double
:longitude_of_central_meridian = -98.0; // double
:standard_parallel = 30.0, 60.0; // double
:earth_radius = 6371229.0; // double
:_CoordinateTransformType = "Projection";
1) :_CoordinateAxisTypes = "GeoX GeoY";
double x(west_east=73);
:units = "km";
:long_name = "synthesized GeoX coordinate from DX attribute";
2) :_CoordinateAxisType = "GeoX";
3) :_CoordinateAliasForDimension = "west_east";
double x_stag(west_east_stag=74);
:units = "km";
:long_name = "synthesized GeoX coordinate from DX attribute";
2) :_CoordinateAxisType = "GeoX";
3) :_CoordinateAliasForDimension = "west_east_stag";
double y(south_north=60);
:units = "km";
:long_name = "synthesized GeoY coordinate from DY attribute";
2) :_CoordinateAxisType = "GeoY";
3) :_CoordinateAliasForDimension = "south_north";
double y_stag(south_north_stag=61);
:units = "km";
:long_name = "synthesized GeoY coordinate from DY attribute";
2) :_CoordinateAxisType = "GeoY";
3) :_CoordinateAliasForDimension = "south_north_stag";
The Lambert variable is a dummy "container variable" for the projection. The next 4 variables are projection coordinates synthesized by the code from the WRF global attributes. The above attributes have the following meanings:
1. The _CoordinateAxisTypes means use this projection on any variable that has a GeoX and GeoY coordinate.
2. The _CoordinateAxisType unambiguously identifies the type of the projection coordinate, in this case either GeoX or GeoY which just mean the projection x or y coordinate.
3. The _CoordinateAliasForDimension means associate this coordinate with any variable that uses the named dimension. This is a convenience so we don't have to annotate all the data variables by hand.
Since we actually want to use CF Conventions instead of the _Coordinate Conventions, here's the equivalent projection coordinates using CF Conventions:
char Lambert;
:grid_mapping_name = "lambert_conformal_conic";
:latitude_of_projection_origin = 34.83000564575195; // double
:longitude_of_central_meridian = -98.0; // double
:standard_parallel = 30.0, 60.0; // double
:earth_radius = 6371229.0; // double
double x(west_east=73);
:units = "km";
:long_name = "synthesized GeoX coordinate from DX attribute";
:axis = "X";
double x_stag(west_east_stag=74);
:units = "km";
:long_name = "synthesized GeoX coordinate from DX attribute";
:axis = "X";
double y(south_north=60);
:units = "km";
:long_name = "synthesized GeoY coordinate from DY attribute";
:axis = "Y";
double y_stag(south_north_stag=61);
:units = "km";
:long_name = "synthesized GeoY coordinate from DY attribute";
:axis = "Y";
Each data variables then needs a coordinates attribute and a grid_mapping attribute:
float U(Time=1, bottom_top=27, south_north=60, west_east_stag=74);
:FieldType = 104; // int
:MemoryOrder = "XYZ";
:description = "x-wind component";
:units = "m s-1";
:stagger = "X";
:grid_mapping = "Lambert";
:coordinates = "y x_stag z Time";
float V(Time=1, bottom_top=27, south_north_stag=61, west_east=73);
:FieldType = 104; // int
:MemoryOrder = "XYZ";
:description = "y-wind component";
:units = "m s-1";
:stagger = "Y";
:grid_mapping = "Lambert";
:coordinates = "y_stag x z Time";
float W(Time=1, bottom_top_stag=28, south_north=60, west_east=73);
:FieldType = 104; // int
:MemoryOrder = "XYZ";
:description = "z-wind component";
:units = "m s-1";
:stagger = "Z";
:grid_mapping = "Lambert";
:coordinates = "y x z_stag Time";
...
You have to go through each data variable and add the grid_mapping attribute and the correct coordinates attribute (the z and time coordinates are covered below). The main difference between _Conventions and CF is that in CF, you must annotate each data variable separately.
2. Vertical Coordinates
Rich's NcML adds these variables to the file:
float ZNU(bottom_top=27);
:FieldType = 104; // int
:MemoryOrder = "Z ";
:description = "eta values on half (mass) levels";
:units = "layer";
:stagger = "";
:positive = "down";
:standard_name = "atmosphere_sigma_coordinate";
:formula_terms = "ptop: P_TOP sigma: ZNU ps: PSFC";
float ZNW(bottom_top_stag=28);
:FieldType = 104; // int
:MemoryOrder = "Z ";
:description = "eta values on full (w) levels";
:units = "level";
:stagger = "Z";
:positive = "down";
:standard_name = "atmosphere_sigma_coordinate";
:formula_terms = "ptop: P_TOP sigma: ZNW ps: PSFC";
These are CF Dimensionless Vertical Coordinates. The CF spec explains in some detail how to use them, and application code can use these formulas to calculate pressure or height coordinates.
The CDM code doesnt know enough to turn the z coordinates into CF Dimensionless Vertical Coordinates, instead it just creates generic Z coordinates:
double z(bottom_top=27);
:units = "";
:long_name = "eta values from variable ZNU";
:_CoordinateAxisType = "GeoZ";
:_CoordinateAliasForDimension = "bottom_top";
double z_stag(bottom_top_stag=28);
:units = "";
:long_name = "eta values from variable ZNW";
:_CoordinateAxisType = "GeoZ";
:_CoordinateAliasForDimension = "bottom_top_stag";
double soilDepth(soil_layers_stag=4);
:units = "units";
:long_name = "soil depth";
:_CoordinateAxisType = "GeoZ";
:_CoordinateAliasForDimension = "soil_layers_stag";
However it includes the soilDepth z coordinate, which the NcML missed. So it looks like the right way to add CF vertical coordinates to this WRF files is:
float ZNU(bottom_top=27);
:FieldType = 104; // int
:MemoryOrder = "Z ";
:description = "eta values on half (mass) levels";
:units = "layer";
:stagger = "";
:positive = "down";
:standard_name = "atmosphere_sigma_coordinate";
:formula_terms = "ptop: P_TOP sigma: ZNU ps: PSFC";
float ZNW(bottom_top_stag=28);
:FieldType = 104; // int
:MemoryOrder = "Z ";
:description = "eta values on full (w) levels";
:units = "level";
:stagger = "Z";
:positive = "down";
:standard_name = "atmosphere_sigma_coordinate";
:formula_terms = "ptop: P_TOP sigma: ZNW ps: PSFC";
double soilDepth(soil_layers_stag=4);
:units = "units";
:long_name = "soil depth";
and also put the appropriate z coordinate name into the data variables' coordinates attribute.
3. Time coordinate
The existing time coordinate in the WRF file looks like:
float XTIME(Time=1);
:FieldType = 104; // int
:MemoryOrder = "0 ";
:description = "minutes since simulation start";
:units = "";
:stagger = "";
data:
{720.0}
For CF, we want to turn this into a udunits compatible date variable, by adding a units attribute:
<variable name="XTIME">
<attribute name="units" value="minutes since 2000-01-24 12:00:00"/>
</variable>
Coordinate variables are preferred, so i would recommend renaming the variable to Time:
float Time(Time=1);
:units = "minutes since 2000-01-24 12:00:00";
:FieldType = 104; // int
:MemoryOrder = "0 ";
:description = "minutes since simulation start";
:stagger = "";
4. Udunit compatible units
CF requires that data variables' units be udunit compatible. The CDM checks the units of the data variables and converts, if possible. This particular file had only a few such corrections:
1. "-" indicating dimensionless units must be an empty string in udunits
2. "W m{-2}" on the NOAHRES variable gets converted to "W m-2"
In conclusion,we have outlined the changes needed to make WRF into CF-compatible netCDF. Next time we will look at some of the features of WRF that aren't well handled in CF and/or by the CDM library.
