Create CF-compliant NetCDF from ESRI Grid?

In ArcGIS 10.1, I just tried converting two rasters (geographic and albers) using ArcToolbox=>Multidimensional Tools=>Raster to NetCDF and the resulting files seem to be CF-Compliant -- they work just fine with the Unidata THREDDS Data Server.

Try drilling down to the OPeNDAP data links or Godiva2 viewer links here: http://geoport.whoi.edu/thredds/catalog/usgs/data2/rsignell/data/bathy/arc_netcdf/catalog.html to see these two netCDF files.

And of course you can just download them using the HTTPServer links.

Here's what the geographic file looks like in Unidata's Tools-UI:

But oh no, here's what the Albers file looks like:

Things don't look good. The topography is not plotting in the right location -- it's shifted by quite a bit.

Taking a look at the ArcGIS 10.1-produced NetCDF file, we see that there is no ellipsoid information in the grid_mapping variable, in this case, the variable albers_conical_equal_area:

netcdf dods://geoport.whoi.edu/thredds/dodsC/usgs/data2/rsignell/data/bathy/arc_netcdf/traster_albers.nc {
  dimensions:
    x = 1272;
    y = 1361;
  variables:
    double x(x=1272);
      :long_name = "x coordinate of projection";
      :standard_name = "projection_x_coordinate";
      :units = "Meter";

    double y(y=1361);
      :long_name = "y coordinate of projection";
      :standard_name = "projection_y_coordinate";
      :units = "Meter";

    int albers_conical_equal_area;
      :grid_mapping_name = "albers_conical_equal_area";
      :longitude_of_central_meridian = -96.0; // double
      :latitude_of_projection_origin = 23.0; // double
      :false_easting = 0.0; // double
      :false_northing = 0.0; // double
      :standard_parallel = 29.5, 45.5; // double

    float traster_View_ProjectRaster(y=1361, x=1272);
      :_CoordinateAxes = "x y y x ";
      :long_name = "traster_View_ProjectRaster";
      :esri_pe_string = "PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-96.0],PARAMETER["Standard_Parallel_1",29.5],PARAMETER["Standard_Parallel_2",45.5],PARAMETER["Latitude_Of_Origin",23.0],UNIT["Meter",1.0]]";
      :coordinates = "x y";
      :grid_mapping = "albers_conical_equal_area";
      :units = "Meter";
      :missing_value = -3.4028235E38f; // float

  // global attributes:
  :Conventions = "CF-1.0";
  :Source_Software = "Esri ArcGIS";
}

The CF Conventions specify that the ellipsoid representation should be encoded using semi_major_axis and (semi_minor_axis or inverse_flattening). Since semi_major_axis and inverse_flattening values are right there in the esri_pe_string attribute, we can fix this file up using NcML thusly:

<netcdf xmlns="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2"
  location="traster_albers.nc">
<variable name="albers_conical_equal_area">
    <attribute name="semi_major_axis" type="double" value="6378137.0"/>
    <attribute name="inverse_flattening" type="double" value="298.257222101"/>
</variable>
</netcdf>

If we try plotting again, it now looks correct:

This should be easy for ESRI to fix because all they have to do is pull the ellipsoid values out of the esri_pe_string. The resulting CF-Compliant file would then look like this:

netcdf dods://geoport.whoi.edu/thredds/dodsC/usgs/data2/rsignell/data/bathy/arc_netcdf/traster_albers_fixed2.ncml {
  dimensions:
    x = 1272;
    y = 1361;
  variables:
    double x(x=1272);
      :long_name = "x coordinate of projection";
      :standard_name = "projection_x_coordinate";
      :units = "Meter";

    double y(y=1361);
      :long_name = "y coordinate of projection";
      :standard_name = "projection_y_coordinate";
      :units = "Meter";

    int albers_conical_equal_area;
      :grid_mapping_name = "albers_conical_equal_area";
      :longitude_of_central_meridian = -96.0; // double
      :latitude_of_projection_origin = 23.0; // double
      :false_easting = 0.0; // double
      :false_northing = 0.0; // double
      :standard_parallel = 29.5, 45.5; // double
      :semi_major_axis = 6378137.0; // double
      :inverse_flattening = 298.257222101; // double

    float traster_View_ProjectRaster(y=1361, x=1272);
      :_CoordinateAxes = "x y y x ";
      :long_name = "traster_View_ProjectRaster";
      :esri_pe_string = "PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-96.0],PARAMETER["Standard_Parallel_1",29.5],PARAMETER["Standard_Parallel_2",45.5],PARAMETER["Latitude_Of_Origin",23.0],UNIT["Meter",1.0]]";
      :coordinates = "x y";
      :grid_mapping = "albers_conical_equal_area";
      :units = "Meter";
      :missing_value = -3.4028235E38f; // float

  // global attributes:
  :Conventions = "CF-1.0";
  :Source_Software = "Esri ArcGIS";
}