TransWikia.com

Create CF-compliant NetCDF from ESRI Grid?

Geographic Information Systems Asked on August 30, 2021

Been trying to make a CF-compliant netCDF file from an ESRI Grid. I want to serve it w/the THREDDS server. Been using gdal_translate. I know the netcdf output option is definitely not CF compliant, so I’ve been using the GMT output option. This doesn’t seem to be functioning properly w/in THREDDS, so I’m guessing that’s not CF-compliant either.

As an aside, I’m trying to validate the CF-compliance of my GMT-type netcdf. Been waiting for a really long time for it to load up to http://puma.nerc.ac.uk/cgi-bin/cf-checker.pl. In the meantime, trying to set up the cf-checker.py script but dependencies on deprecated CDAT libraries is making that a pain.

Anybody got a cleaner way to do this?

—EDIT—

ncdump of my GMT-conforming NetCDF file. Seems odd that I have "x_range" and "y_range" but no "x" or "y". Also seems like I’m providing almost no metadata via this process–making CF compliance seem unlikely to me…

netcdf E:/software/apache-tomcat-6.0.39/content/thredds/public/wk_gmt/wk_awc/al.nc {
  dimensions:
    side = 2;
    xysize = 1902411900;
  variables:
    double x_range(side=2);
      :units = "meters";

    double y_range(side=2);
      :units = "meters";

    double z_range(side=2);
      :units = "meters";

    double spacing(side=2);

    int dimension(side=2);

    float z(xysize=1902411900);
      :scale_factor = 1.0; // double
      :add_offset = 0.0; // double
      :node_offset = 1; // int

  // global attributes:
  :title = "";
  :source = "";
}

—EDIT # 2—

Just a little more content to show the difference in the NetCDF files produced by gdal vs. ESRI based on comment from @mdsumner. I totally agree that gdal should be doing it right. Quite possible I’m not using the tool correctly (gdal_translate -ot Float32 -of NetCDF ${in} ${out} in my shell script, fwiw).

As @signell notes below, the ESRI output ain’t perfect, but it seems closer to CF compliance…the following is the ncdump of the file created from a *.nc from the same ESRI Grid. This file was created using ESRI’s Raster to NetCDF.

netcdf E:/software/apache-tomcat-6.0.39/content/thredds/public/wk_esri/wk_awc/al.nc {
  dimensions:
    x = 11334;
    y = 18650;
  variables:
    double x(x=11334);
      :long_name = "x coordinate of projection";
      :standard_name = "projection_x_coordinate";
      :units = "Meter";

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

    float awc(y=18650, x=11334);
      :long_name = "awc";
      :esri_pe_string = "PROJCS["NAD_1983_Albers",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

    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

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

One Answer

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: enter image description here

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

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: enter image description here

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";
}

Correct answer by Rich Signell on August 30, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP