Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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";
}
share|improve this question
1  
GMT is a really weird variant that drops the array dims and requires the client to reinstate it, I doubt that it's what you want –  mdsumner Feb 7 at 21:50
    
+1 @mdsumner for confirming that. Thanks. It was a guess based on a 6-year old post to another wikispace. Desperate measure at best! –  Roland Feb 7 at 22:11
    
I think you should demonstrate your evidence for this and what version of GDAL you are using and how you use it, " I know the netcdf output option is definitely not CF compliant" - it's not something I've looked at in detail but a lot of work has been done to improve NetCDF in GDAL in the last few years. –  mdsumner Feb 7 at 22:48
    
although I'm no expert, I felt the ncdump that I edited into the question earlier shows at least one major issue, which is that there are no "x" and "y" variables. More heuristically, THREDDS wasn't able to produce a usable WMS with that. I'll also include the ncdump of file created using Raster to NetCDF for comparison. –  Roland Feb 7 at 23:08
1  
GDAL version 1.9.2 works perfectly when you use -of netCDF, generating not only the coordinate variables x and y, but also the CF grid_mapping variable complete with ellipsoid information semi_major_axis and inverse_flattening. –  Rich Signell Feb 8 at 14:23

1 Answer 1

up vote 3 down vote accepted

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";
}
share|improve this answer
    
and extra +1 for the workaround that can be used until ESRI patches this tool! –  Roland Feb 7 at 22:11
    
Thanks Rich, great answer. Hopefully something that can be fixed in the next release. –  scw Feb 13 at 2:35

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.