THREDDS and HDF

Progress in the matter of serving LSU/ESL and USF/IMaRS HDF satellite products via a THREDDS WCS server.



LSU Example

  • NetCDF Subset Service

  • WCS Service

    Python Software for Reading and Writing HDF Files
    Converting HDF4 to HDF5
    Annotated Version of albers.xml Example



    The LSU Examples

    LSU file examples can be obtained from:

    http://www.esl.lsu.edu/webdata/poes/GCOOS/2010-06/

    The LSU THREDDS server is at:

    http://www.esl.lsu.edu:8080/thredds/catalog.html

    A directory containing files being served via the THREDDS OpenDAP server is:

    http://www.esl.lsu.edu:8080/thredds/catalog/GCOOS/2010-03/catalog.html

    An example "OPeNAP Dataset Access Form" is:

    http://www.esl.lsu.edu:8080/thredds/dodsC/GCOOS/2010-03/n17.100302.1536.cloudmask.hdf.html

    The Source HDF Files

    LSU is serving files with and without a cloudmask.

    The first test file from LSU is n15.1006161124.cloudmask.hdf.

    The dimensions, variables and global attributes sections of this file as obtained with the h4dump command are:

    netcdf n15.100616.1124.cloudmask {
    dimensions:
    	line = 1600 ;
    	sample = 2400 ;
    
    variables:
    	float mcsst(line, sample) ;
    		mcsst:long_name = "mcsst" ;
    		mcsst:units = "temp_deg_c" ;
    		mcsst:valid_range = -3.4028235e+38f, 3.4028235e+38f ;
    		mcsst:_FillValue = -3.4028235e+38f ;
    		mcsst:scale_factor = 1. ;
    		mcsst:scale_factor_err = 0. ;
    		mcsst:add_offset = 0. ;
    		mcsst:add_offset_err = 0. ;
    		mcsst:calibrated_nt = 5 ;
    	float line(line) ;
    		line:long_name = "line" ;
    	float sample(sample) ;
    		sample:long_name = "sample" ;
    
    // global attributes:
    		:line\coord = "y" ;
    		:line\scale = 1. ;
    		:line\offset = 0. ;
    		:sample\coord = "x" ;
    		:sample\scale = 1. ;
    		:sample\offset = 0. ;
    		:projection_name = "rectangular" ;
    		:satellite = "noaa-15" ;
    		:sensor_name = "avhrr" ;
    		:pass_date = "2010/06/16" ;
    		:start_time = "11:27:54.398" ;
    		:center_lat = 25. ;
    		:center_lon = -87. ;
    		:upper_left_lat = 32.18652227295618 ;
    		:upper_right_lat = 32.18652227295618 ;
    		:upper_left_lon = -98.89417498516796 ;
    		:upper_right_lon = -75.10582501483204 ;
    		:lower_left_lat = 17.81347772704383 ;
    		:lower_right_lat = 17.81347772704383 ;
    		:lower_left_lon = -98.89417498516796 ;
    		:lower_right_lon = -75.10582501483204 ;
    		:spatial_resolution = 1000 ;
    		:spatial_units = "m" ;
    		:capture_site = "LSU Earth Scan Laboratory" ;
    		:capture_system = "Seaspace Terascan" ;
    
    data:
    
     mcsst =
    ...
    

    The second test file is n15.1006161124.no_cloudmask.hdf.

    The header for this file is:

    netcdf n15.100616.1124.no_cloudmask {
    dimensions:
    	line = 1600 ;
    	sample = 2400 ;
    
    variables:
    	float mcsst_no_cloudmask(line, sample) ;
    		mcsst_no_cloudmask:long_name = "mcsst_no_cloudmask" ;
    		mcsst_no_cloudmask:units = "temp_deg_c" ;
    		mcsst_no_cloudmask:valid_range = -3.4028235e+38f, 3.4028235e+38f ;
    		mcsst_no_cloudmask:_FillValue = -3.4028235e+38f ;
    		mcsst_no_cloudmask:scale_factor = 1. ;
    		mcsst_no_cloudmask:scale_factor_err = 0. ;
    		mcsst_no_cloudmask:add_offset = 0. ;
    		mcsst_no_cloudmask:add_offset_err = 0. ;
    		mcsst_no_cloudmask:calibrated_nt = 5 ;
    	float line(line) ;
    		line:long_name = "line" ;
    	float sample(sample) ;
    		sample:long_name = "sample" ;
    
    // global attributes:
    		:line\coord = "y" ;
    		:line\scale = 1. ;
    		:line\offset = 0. ;
    		:sample\coord = "x" ;
    		:sample\scale = 1. ;
    		:sample\offset = 0. ;
    		:projection_name = "rectangular" ;
    		:satellite = "noaa-15" ;
    		:sensor_name = "avhrr" ;
    		:pass_date = "2010/06/16" ;
    		:start_time = "11:27:54.398" ;
    		:center_lat = 25. ;
    		:center_lon = -87. ;
    		:upper_left_lat = 32.18652227295618 ;
    		:upper_right_lat = 32.18652227295618 ;
    		:upper_left_lon = -98.89417498516796 ;
    		:upper_right_lon = -75.10582501483204 ;
    		:lower_left_lat = 17.81347772704383 ;
    		:lower_right_lat = 17.81347772704383 ;
    		:lower_left_lon = -98.89417498516796 ;
    		:lower_right_lon = -75.10582501483204 ;
    		:spatial_resolution = 1000 ;
    		:spatial_units = "m" ;
    		:capture_site = "LSU Earth Scan Laboratory" ;
    		:capture_system = "Seaspace Terascan" ;
    
    data:
    
     mcsst_no_cloudmask =
    ...
    

    The ToolsUI NcML File Dump of the HDF4 File

    If we ingest the preceding HDF file into ToolsUI, we obtain the following NCML file:

    <?xml version="1.0" encoding="UTF-8"?>
    <netcdf xmlns="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2" location="file:/d1/THREDDS/MODIS.2009037.191026.gcoos.seadas_sst.hdf">
      <dimension name="fakeDim0" length="1430" />
      <dimension name="fakeDim1" length="2090" />
      <dimension name="fakeDim2" length="1430" />
      <dimension name="fakeDim3" length="2090" />
      <dimension name="fakeDim4" length="1430" />
      <dimension name="fakeDim5" length="2090" />
      <attribute name="Map Projection Category" value="Proj 4.4.7" />
      <attribute name="Data Source Type" value="Direct Broadcast" />
      <attribute name="L3 Software Name" value="/usr/bin/map_modis.pl" />
      <attribute name="L3 Software Version" value="4.2" />
      <attribute name="Map Limit" value="18, -98, 31, -79" />
      <attribute name="Processing Status" value="final" />
      <attribute name="Collection Location" value="USF Institute for Marine Remote Sensing" />
      <attribute name="Satellite" value="Aqua" />
      <attribute name="Input File(s)" value="/tmp/s4p_run_map_modis_2.3_ ..." />
      <attribute name="Map Projection Name" value="eqc" />
      <attribute name="Proj Projection Parameters" value="#Equidistant Cylindrical (Plate Caree)
    #	Cyl, Sph
    #	lat_ts=
    # +proj=eqc +ellps=WGS84" />
      <attribute name="Temporal Resolution" value="pass" />
      <attribute name="Spatial Resolution" value="1000" />
      <attribute name="History" value="Direct read of HDF4 file through CDM library" />
      <attribute name="HDF4_Version" value="4.2.1 (NCSA HDF Version 4.2 Release 1, February 17, 2005)" />
      <variable name="seadas_sst" shape="fakeDim0 fakeDim1" type="short">
        <attribute name="long_name" value="Sea Surface Temperature" />
        <attribute name="slope" value="[  0.005]" />
        <attribute name="units" value="degrees-C" />
        <attribute name="intercept" value="[0]" />
      </variable>
      <variable name="l2_flags" shape="fakeDim2 fakeDim3" type="int">
        <attribute name="f02_name" value="LAND" />
        <attribute name="long_name" value="Level-2 Processing Flags" />
        <attribute name="f32_name" value="OCEAN" />
        <attribute name="f03_name" value="BADANC" />
        <attribute name="f21_name" value="MODGLINT" />
        <attribute name="f10_name" value="CLDICE" />
        <attribute name="f05_name" value="HILT" />
        <attribute name="f06_name" value="HISATZEN" />
        <attribute name="f07_name" value="COASTZ" />
        <attribute name="f08_name" value="NEGLW" />
        <attribute name="f25_name" value="SEAICE" />
        <attribute name="f20_name" value="MAXAERITER" />
        <attribute name="f13_name" value="HISOLZEN" />
        <attribute name="f31_name" value="SPARE" />
        <attribute name="f19_name" value="TRICHO" />
        <attribute name="f18_name" value="ABSAER" />
        <attribute name="f30_name" value="HIPOL" />
        <attribute name="f04_name" value="HIGLINT" />
        <attribute name="f17_name" value="NAVWARN" />
        <attribute name="f01_name" value="ATMFAIL" />
        <attribute name="f26_name" value="NAVFAIL" />
        <attribute name="f24_name" value="DARKPIXEL" />
        <attribute name="f16_name" value="CHLFAIL" />
        <attribute name="f28_name" value="SSTWARN" />
        <attribute name="f15_name" value="LOWLW" />
        <attribute name="f11_name" value="COCCOLITH" />
        <attribute name="f12_name" value="TURBIDW" />
        <attribute name="f23_name" value="ATMWARN" />
        <attribute name="f09_name" value="STRAYLIGHT" />
        <attribute name="f22_name" value="CHLWARN" />
        <attribute name="f29_name" value="SSTFAIL" />
        <attribute name="f27_name" value="FILTER" />
        <attribute name="f14_name" value="HITAU" />
        <attribute name="units" value="dimensionless" />
      </variable>
      <variable name="cloud_mask" shape="fakeDim4 fakeDim5" type="byte">
        <attribute name="_Unsigned" value="true" />
      </variable>
    </netcdf>
    

    Here is successful attempt by Rich Signell to accomplish the task:

    The NcML File for Serving the SeaDAS SST Files

    <?xml version="1.0" encoding="UTF-8"?>
    
    <catalog name="GCOOS Tests"
      xmlns="http://www.unidata.ucar.edu/namespaces/thredds/InvCatalog/v1.0"
      xmlns:xlink="http://www.w3.org/1999/xlink">
      <service name="allServices" serviceType="Compound" base="">
        <service name="ncdods" serviceType="OPENDAP" base="/thredds/dodsC/"/>
        <service name="wcs" serviceType="WCS" base="/thredds/wcs/"/>
        <service name="ncss" serviceType="NetcdfSubset" base="/thredds/ncss/grid/"/>
        <service name="HTTPServer" serviceType="HTTPServer" base="/thredds/fileServer/"/>
      </service>
      
      <dataset name="GCOOS SeaDAS SST" ID="gcoos/seadas_sst" serviceName="allServices" urlPath="gcoos/seadas_sst">
        
        <netcdf xmlns="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2">
    
          <aggregation dimName="time" type="joinNew" recheckEvery="5 min">
            <variableAgg name="seadas_sst"/>
            <variableAgg name="cloud_mask"/>
            <variableAgg name="l2_flags"/>
            <scan location="/data1/TEST/SeaDAS_SST" regExp=".*MODIS.*\.gcoos.seadas_sst.hdf$" dateFormatMark="MODIS.#yyyyDDD'.'hhmmss" subdirs="false"/>
          </aggregation>
    
          <dimension name="lon" orgName="fakeDim1"/>
          <dimension name="lat" orgName="fakeDim0"/>
    
          <variable name="lon" shape="lon" type="double">
              <attribute name="units" value="degrees_east"/>
              <values start="-98.0" increment="0.00909"/>
          </variable>
    
          <variable name="lat" shape="lat" type="double">
              <attribute name="units" value="degrees_north"/>
              <values start="18.0" increment="0.00909"/>
          </variable>
    
          <variable name="seadas_sst" shape="time lat lon" type="short">
              <attribute name="standard_name" value="sea_surface_temperature"/>
              <attribute name="units" value="Celcius"/>
              <attribute name="scale_factor" type="float" value="0.005"/>
              <attribute name="add_offset" type="float" value="0.0"/>
          </variable>
    
          <variable name="l2_flags" shape="time lat lon" type="int">
              <attribute name="units" value="1"/>
          </variable>
    
          <variable name="cloud_mask" shape="time lat lon" type="byte">
              <attribute name="units" value="1"/>
          </variable>
    
        </netcdf>
    
      </dataset>
    </catalog>
    


    The NetCDF Subset Interface

    We now access a file via the NetCDF subset interface described at:

    http://www.unidata.ucar.edu/projects/THREDDS/tech/interfaceSpec/NetcdfSubsetService.html

    The TDS interface to the NetCDF subset interface to the SeaDAS SST example is at:

    http://megara.tamu.edu:8080/thredds/ncss/grid/gcoos/seadas_sst/dataset.html

    This appears to be working, at least for today.


    The WCS Interface

    We will now access a file via the WCS interface as described in the TDS Web Coverage Service documentation at:

    http://www.unidata.ucar.edu/projects/THREDDS/tech/reference/WCS.html

    The GetCapabilities Request

    We can obtain basic information about the available capabilities via the GetCapabilities request:

    http://megara.tamu.edu:8080/thredds/wcs/gcoos/seadas_sst?request=GetCapabilities&version=1.0.0&service=WCS

    will obtain the ContentMetadata information for seadas_sst:

    <ContentMetadata>
      <CoverageOfferingBrief>
        <description>
          seadas_sst     Celsius        false     Sea Surface Temperature
        </description>
        <name>seadas_sst</name>
        <label>Sea Surface Temperature</label>
        <lonLatEnvelope srsName="urn:ogc:def:crs:OGC:1.3:CRS84">
          <gml:pos>-98.0 18.0</gml:pos>
          <gml:pos>-79.01098999999999 30.98961</gml:pos>
          <gml:timePosition>2009-10-23T16:25:27Z</gml:timePosition>
          <gml:timePosition>2009-10-23T16:25:27Z</gml:timePosition>
        </lonLatEnvelope>
      </CoverageOfferingBrief>
    </ContentMetadata>
    

    The DescribeCoverage Request

    We can then obtain more detailed information about the available variables via the DescribeCoverage request.

    http://megara.tamu.edu:8080/thredds/wcs/gcoos/seadas_sst?request=DescribeCoverage&version=1.0.0&service=WCS

    This provides the ConverageOffering information for seadas_sst:

    <CoverageOffering>
      <description>
        seadas_sst     Celsius        false     Sea Surface Temperature
      </description>
      <name>seadas_sst</name>
      <label>Sea Surface Temperature</label>
      <lonLatEnvelope srsName="urn:ogc:def:crs:OGC:1.3:CRS84">
        <gml:pos>-98.0 18.0</gml:pos>
        <gml:pos>-79.01098999999999 30.98961</gml:pos>
        <gml:timePosition>2009-10-23T16:25:27Z</gml:timePosition>
        <gml:timePosition>2009-10-23T16:25:27Z</gml:timePosition>
      </lonLatEnvelope>
      <domainSet>
        <spatialDomain>
          <EnvelopeWithTimePeriod srsName="urn:ogc:def:crs:OGC:1.3:CRS84">
            <gml:pos dimension="2">-98.0 18.0</gml:pos>
            <gml:pos dimension="2">-79.01098999999999 30.98961</gml:pos>
            <gml:timePosition>2009-10-23T16:25:27Z</gml:timePosition>
            <gml:timePosition>2009-10-23T16:25:27Z</gml:timePosition>
          </EnvelopeWithTimePeriod>
          <gml:RectifiedGrid srsName="OGC:CRS84" dimension="2">
            <gml:limits>
              <gml:GridEnvelope>
                <gml:low>0 0</gml:low>
                <gml:high>2089 1429</gml:high>
              </gml:GridEnvelope>
            </gml:limits>
            <gml:axisName>x</gml:axisName>
            <gml:axisName>y</gml:axisName>
            <gml:origin>
              <gml:pos>-98.0 18.0</gml:pos>
            </gml:origin>
            <gml:offsetVector>0.009090000000000004 0.0</gml:offsetVector>
            <gml:offsetVector>0.0 0.009089999999999999</gml:offsetVector>
          </gml:RectifiedGrid>
        </spatialDomain>
        <temporalDomain>
          <gml:timePosition>2009-10-23T16:25:27Z</gml:timePosition>
        </temporalDomain>
      </domainSet>
      <rangeSet>
        <RangeSet>
          <description>
            seadas_sst     Celsius        false     Sea Surface Temperature
          </description>
          <name>seadas_sst</name>
          <label>Sea Surface Temperature</label>
        </RangeSet>
      </rangeSet>
      <supportedCRSs>
        <requestCRSs>OGC:CRS84</requestCRSs>
        <responseCRSs>OGC:CRS84</responseCRSs>
      </supportedCRSs>
      <supportedFormats>
        <formats>GeoTIFF</formats>
        <formats>GeoTIFF_Float</formats>
        <formats>NetCDF3</formats>
      </supportedFormats>
      <supportedInterpolations>
        <interpolationMethod>none</interpolationMethod>
      </supportedInterpolations>
    </CoverageOffering>
    

    The GetCoverage Request

    A GetCoverage request following the information thus obtained is:

    curl -o subset.nc 'http://megara.tamu.edu:8080/thredds/wcs/gcoos/seadas_sst?service=WCS&version=1.0.0&request=GetCoverage&format=NetCDF3&coverage=seadas_sst&bbox=-95.0,23.0,-85.0,28.0'

    This request is for the single HDF file in the /data1/TEST/SeaDAS_SST subdirectory, i.e.

    MODIS.2009296.162527.gcoos.seadas_sst.hdf

    and returns the file:

    subset.nc

    The corresponding request for a GeoTIFF file:

    curl -o subset.tiff 'http://megara.tamu.edu:8080/thredds/wcs/gcoos/seadas_sst?service=WCS&version=1.0.0&request=GetCoverage&format=GeoTIFF&coverage=seadas_sst&bbox=-95.0,23.0,-85.0,28.0'

    obtains the following file:

    subset.tiff

    This was developed analogous to the following WCS request produced by Rich.

    curl -o subset.nc 'http://blackburn.whoi.edu:8081/thredds/wcs/coastwatch/albers3?service=WCS&version=1.0&request=GetCoverage&format=NetCDF3&coverage=chlor_a&bbox=-92.7,28,-88,30.'


    Python Software for Reading and Writing HDF4 Files

    Since the SeaDAS HDF4 files cannot be served as is via the THREDDS WCS server, a program must be built that can transform them from their present into a form that can be served via THREDDS. A good candidate for this task is pyhdf, a package that can read and write HDF4 files from within Python scripts.


    Installing pyhdf

    First download pyhdf from:

    http://sourceforge.net/projects/pysclint/files/

    The latest version is 0.8.3. Unpack this via:

    
    tar xzvf pyhdf-0.8.3.tar.gz
    cd pyhdf-0.8.3
    
    
    The installation script must know the location of the HDF4 libraries. This is done by setting a couple of environmental variables:
    
    export INCLUDE_DIRS=/opt/HDF/hdf4.2r4/include
    export LIBRARY_DIRS=/opt/HDF/hdf4.2r4/lib
    
    The installation can proceed via:

    su
    python setup.py install


    Using pyhdf

    The main documentation page at:

    http://pysclint.sourceforge.net/pyhdf/documentation.html

    contains a link to the documentation for the module that implements the SD or scientific dataset API of the HDF4 library at:

    http://pysclint.sourceforge.net/pyhdf/pyhdf.SD.html

    It is with this module that we will be able to read the original HDF4 file and transform it into a new file that is compatible with the THREDDS WCS server.

    Reading an HDF4 File

    The following code script is a model demonstrating how to read an HDF4 using pyhdf.

    
    # Import SD and numpy.
    from pyhdf.SD import *
    from numpy import *
     
    fileName = 'template.hdf'
    # Open file in read-only mode (default)
    hdfFile = SD(fileName)
    # Display attributes.
    print "file:", fileName
    print "author:", hdfFile.author
    print "priority:", hdfFile.priority
    # Open dataset 'd1'
    d1 = hdfFile.select('d1')
    # Display dataset attributes.
    print "dataset:", 'd1'
    print "description:",d1.description
    print "units:", d1.units
    # Display dimensions info.
    dim1 = d1.dim(0)
    dim2 = d1.dim(1)
    print "dimensions:"
    print "dim1: name=", dim1.info()[0], 
    print "length=", dim1.length(), 
    print "units=", dim1.units
    print "dim2: name=", dim2.info()[0], 
    print "length=", dim2.length(), 
    print "units=", dim2.units
    # Show dataset values
    print d1[:]
    # Close dataset
    d1.endaccess()
    # Close file
    hdfFile.end()
    

    Protytpe Code for Converting from SeaDAS HDF to THREDDS WCS HDF

    
    # Create output file.
    outname = 'test.hdf'
    outfile = SD(outname ,SDC.WRITE|SDC.CREATE)
    
    # Open input file.
    infile = SD('MODIS.2009037.191026.gcoos.seadas_sst.hdf')
    
    # Open the seadas_sst dataset.
    sst = infile.select('seadas_sst')
    
    # Obtain the number of lon and lat points.
    dim1 = sst.dim(0)
    dim2 = sst.dim(1)
    nlat = dim1.info()[1]
    nlon = dim2.info()[1]
    
    # Extract the data values from sst.
    seadas_sst = sst[:]
    
    # Extract all global attributes.
    #attr = infile.attributes(full=1)
    
    # Extract the attribute names.
    #attrnames = attr.keys()
    
    # Print an attribute name.
    #attnames[0]
    #'Map Projection Name'
    
    # Get the map limits from the attributes.
    geolimits = attr['Map Limit'][0]
    #geolimits
    #'18, -98, 31, -79'
    geolimits = string.split(geolimits)
    #geolimits
    #['18,', '-98,', '31,', '-79']
    geolimits = geolimits.replace(',','')
    #geolimits
    #'18 -98 31 -79'
    
    # Obtain the lon and lat min and max.
    latmin = float(geolimits[0])
    latmax = float(geolimits[2])
    lonmin = float(geolimits[1])
    lonmax = float(geolimits[3])
    
    # Create 1-D arrays for lon and lat.
    lon1d = linspace(lonmin,lonmax,nlon)
    lat1d = linspace(latmax,latmin,nlat)
    
    # Open HDF datasets for 2D lon and lat arrays.
    lon2d = outfile.create('lon2d', SDC.FLOAT32, (1430,2090))
    lat2d = outfile.create('lat2d', SDC.FLOAT32, (1430,2090))
    #lon2d = outfile.create('lon2d', SDC.FLOAT32, (nlat,nlon))
    #lat2d = outfile.create('lat2d', SDC.FLOAT32, (nlat,nlon))
    
    # Name the dimensions and assign attributes.
    d1_lon = lon2d.dim(0)
    d2_lon = lon2d.dim(1)
    d1_lon.setname('rows')
    d2_lon.setname('cols')
    d1_lat = lat2d.dim(0)
    d2_lat = lat2d.dim(1)
    d1_lat.setname('rows')
    d2_lat.setname('cols')
    
    # Create 2-D arrays for lon and lat by assigning values to lon2d and lat2d.
    lon2d = vstack((nlat*[lon1d[:]]))
    lat2dpre = vstack((nlon*[lat1d[:]]))
    lat2d = swapaxes(lat2dpre,0,1)
    
    # Create attributes for lon2d and lat2d.
    lon2d.long_name = 'longitude'
    lon2d.units = 'degrees'
    lon2d.coordsys = 'Equidistant Cylindrical'
    lat2d.long_name = 'latitude'
    lat2d.units = 'degrees'
    lat2d.coordsys = 'Equidistant Cylindrical'
    
    # Close datasets lon2d and lat2d.
    lon2d.endaccess()
    lat2d.endaccess()
    
    # ============================================================
    # OUTPUT FILE STUFF
    # ============================================================
    
    # Open HDF dataset for 2D seadas_sst field.
    sst_out = outfile.create('sst_out', SDC.FLOAT32, (1430,2090))
    
    # Name the dimensions and assign attributes.
    d1_sst = sst_out.dim(0)
    d2_sst = sst_out.dim(1)
    d1_sst.setname('rows')
    d2_sst.setname('cols')
    
    # Write SST field to output file by assigning values to sst_out:
    sst_out[:] = seadas_sst
    
    # Create attributes for sst_out.
    sst_out.long_name = 'Sea Surface Temperature'
    sst_out.units = 'degrees-C'
    sst_out.slope = '[0.005]'
    sst_out.intercept = '[0]'
    
    # Close dataset sst_out.
    sst_out.endaccess()
    
    # Create attributes for the file.
    outfile.satellite = 'Aqua'
    outfile.projection = 'eqc'
    outfile.projection_parameters = 'Equidistant Cylindrical'
    outfile.proj4_version = '4.4.7'
    outfile.data_source = 'Direct Broadcast'
    outfile.L3_software_name = '/usr/bin/map_modis.pl'
    outfile.L3_software_version = '4.2'
    outfile.processing_status = 'final'
    outfile.collection_location = 'USF Institute for Marine Remote Sensing'
    outfile.temporal_resolution = 'pass'
    outfile.spatial_resolution = '1000 m'
    
    # Close input and output files.
    infile.end()
    outfile.end()
    
    

    Writing an HDF4 File

    The following code script is a model demonstrating how to write an HDF4 file using pyhdf.

    # Import SD and numpy.
    from pyhdf.SD import *
    from numpy import *
    
    fileName = 'template.hdf'
    # Create HDF file.
    hdfFile = SD(fileName ,SDC.WRITE|SDC.CREATE)
    # Assign a few attributes at the file level
    hdfFile.author = 'It is me...'
    hdfFile.priority = 2
    # Create a dataset named 'd1' to hold a 3x3 float array.
    d1 = hdfFile.create('d1', SDC.FLOAT32, (3,3))
    # Set some attributs on 'd1'
    d1.description = 'Sample 3x3 float array'
    d1.units = 'celsius'
    # Name 'd1' dimensions and assign them attributes.
    dim1 = d1.dim(0)
    dim2 = d1.dim(1)
    dim1.setname('width')
    dim2.setname('height')
    dim1.units = 'm'
    dim2.units = 'cm'
    # Assign values to 'd1'
    d1[0]  = (14.5, 12.8, 13.0)  # row 1
    d1[1:] = ((-1.3, 0.5, 4.8),  # row 2 and
              (3.1, 0.0, 13.8))  # row 3
    # Close dataset
    d1.endaccess()
    # Close file
    hdfFile.end()
    



    Converting HDF4 to HDF5

    Note: While originally considered as an option when it was thought that the latest version of THREDDS could only process HDF5 files, this section is now relegated to the "might be somehow useful in the future" category since HDF4 files can indeed be processed by the latest THREDDS version.

    A basic problem is that while THREDDS Data Server 4 can handle HDF5 datasets via the included NetCDF Java Library, it cannot handle HDF4 datasets. A possible shortcut to nirvana is the conversion library for HDF4/HDF5, which can supposedly convert any HDF4 file to HDF5. The conversion utility h4toh5 is built using this library and automatically converts HDF4 files to HDF5 files. It can also convert HDF-EOS2 files to HDF5 files that can be accessed by NetCDF4.

    Installing the Conversion Library

    The conversion library can be obtained at:

    http://www.hdfgroup.org/h4toh5/

    with the present source code version being h4h5tools-2.1.tar.gz. The installation instructions are slightly buried in the subdirectory release_docs in the file INSTALL_Unix.txt. As one might suspect, both the HDF4 and HDF5 libraries are need to compile this, but they are included in a non-standard way as explained in the file. An additional compilation parameter is also needed if you're using HDF 1.8 or higher. The three flags for, respectively, the the HDF 1.8 or higher version, the HDF4 library, and the HDF5 library are:

    CPPFLAGS=-DH5_USE_16_API

    CC=<hdf4-directory>/bin/h4cc

    --with-hdf5=%lt;hdf5-directory>

    In our case, the HDF4 library is in /opt/HDF/hdf4.2r4 and the HDF5 library is in /opt/HDF/hdf5-1.8.1, so the last two parameters will be:

    CC=/opt/HDF/hdf4.2r4/bin/h4cc

    --with-hdf5=/opt/HDF/hdf5-1.8.1

    The entire configuration string - assuming additionally that /usr is the installation directory - is:

    ./configure --prefix=/usr CPPFLAGS=-DH5_USE_16_API CC=/opt/HDF/hdf4.2r4/bin/h4cc --with-hdf5=/opt/HDF/hdf5-1.8.1

    After compilation, install this via:

    su; make install

    to obtain the two binaries h5toh4 and h4toh5 and the library libh4toh5.a.

    Conversion Example

    We now convert the HDF4 file (shown below) MODIS.2009037.191026.gcoos.seadas_sst.hdf into HDF5 using h4toh5 via:

    h4toh5 MODIS.2009037.191026.gcoos.seadas_sst.hdf test.hdf5

    There are significant differences in the way the header metainformation is handled between HDF4 and HDF5. An example of the differences can be see with the seadas_sst variable.

    HDF4 Header In the HDF4 file, the header information is of the form:

    dimensions:
            fakeDim0 = 1430 ;
            fakeDim1 = 2090 ;
    ...
    variables:
            short seadas_sst(fakeDim0, fakeDim1) ;
                    seadas_sst:long_name = "Sea Surface Temperature" ;
                    seadas_sst:slope = "[  0.005]" ;
                    seadas_sst:units = "degrees-C" ;
                    seadas_sst:intercept = "[0]" ;
    ...
     seadas_sst =
      -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999,
        -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999,
    ...
    

    HDF5 Header

    In the HDF5 file, the header information for the same variable is of the form:

    ...
       DATASET "seadas_sst" {
          DATATYPE  H5T_STD_I16BE
          DATASPACE  SIMPLE { ( 1430, 2090 ) / ( 1430, 2090 ) }
          DATA {
          (0,0): -999, -999, -999, -999, -999, -999, -999, -999, -999, -999,
          (0,10): -999, -999, -999, -999, -999, -999, -999, -999, -999, -999,
    ...
          (1429,2085): 3538, 3623, 3623, 4716, 4716
          }
          ATTRIBUTE "HDF4_OBJECT_NAME" {
             DATATYPE  H5T_STRING {
                   STRSIZE 10;
                   STRPAD H5T_STR_SPACEPAD;
                   CSET H5T_CSET_ASCII;
                   CTYPE H5T_C_S1;
                }
             DATASPACE  SCALAR
             DATA {
             (0): "seadas_sst"
             }
          }
    ...
          ATTRIBUTE "long_name" {
             DATATYPE  H5T_STRING {
                   STRSIZE 24;
                   STRPAD H5T_STR_SPACEPAD;
                   CSET H5T_CSET_ASCII;
                   CTYPE H5T_C_S1;
                }
             DATASPACE  SCALAR
             DATA {
             (0): "Sea Surface Temperature\000"
             }
          }
          ATTRIBUTE "slope" {
             DATATYPE  H5T_STRING {
                   STRSIZE 9;
                   STRPAD H5T_STR_SPACEPAD;
                   CSET H5T_CSET_ASCII;
                   CTYPE H5T_C_S1;
                }
             DATASPACE  SCALAR
             DATA {
             (0): "[  0.005]"
             }
          }
          ATTRIBUTE "units" {
             DATATYPE  H5T_STRING {
                   STRSIZE 10;
                   STRPAD H5T_STR_SPACEPAD;
                   CSET H5T_CSET_ASCII;
                   CTYPE H5T_C_S1;
                }
             DATASPACE  SCALAR
             DATA {
             (0): "degrees-C\000"
             }
          }
          ATTRIBUTE "intercept" {
             DATATYPE  H5T_STRING {
                   STRSIZE 3;
                   STRPAD H5T_STR_SPACEPAD;
                   CSET H5T_CSET_ASCII;
                   CTYPE H5T_C_S1;
                }
             DATASPACE  SCALAR
             DATA {
             (0): "[0]"
             }
          }
    ...
    



    Applicable Correspondence:


    Xiaoming,

    I need to get two Gulf of Mexico satellite groups (LSU/ESL and USF/IMaRS) serving selected satellite products with a WCS protocol. The TDS approach seems most promising at the present moment. (I know the OPeNDAP group is working on a OPeNDAP-WCS gateway).

    Both of my satellite groups use HDF (HDF4 I think) as their native storage format. I'm sure that the complication factor for me will be that their georeferencing metadata is either stored in the wrong place or is not CF, or some such.

    How did you get around this? I know that you (Rich) have been using NcML to get numerical model data and/or georeferencing metadata in the right way/place. Model data - NetCDF, NcML is a NetCDF tool.

    How do I coerce HDF contents into the proper form without wholesale recreation/duplication of their satellite archives? NcML can't be made to work on HDF files - can it?

    Is there a write up of what Coastwatch did to serve their products via a WCS protocol?

    Matt


    Notes:

    Sample files from USF/IMaRS

    The sample file locations can be found in an emial from Feb. 12, 2009 (forwarded Jun. 29, 2009) wherein we find:

    Matt & Gary,

    Actually, the HDF files contain far more than just the "chlor_a" data. This may be very difficult. I have this uneasy feeling that "if the glove does not fit, you cannot convict". I do not know of other conversations on this matter.

    For example, here is a representative HDF:

    http://www.imars.usf.edu/dods-bin/nph-dods/modis/level3/husf/gcoos/2009/037/1km/pass/final/MODIS.2009037.191026.gcoos.seadas.hdf.html

    In each GCOOS SeaDAS HDF file there are the globals:
    ...

    Brock Murch

    That URL is for the DODS interface. The file directory can be found at the corresponding URL:

    http://www.imars.usf.edu/dods-bin/nph-dods/modis/level3/husf/gcoos/2009/037/1km/pass/final/

    The SeaDAS files themselves can be found at corresponding URLs, e.g. for the above example:

    http://cyclops.marine.usf.edu/modis/level3/husf/gcoos/2009/037/1km/pass/final/MODIS.2009037.191026.gcoos.seadas.hdf

    and the corresponding SeaDAS SST file at:

    http://cyclops.marine.usf.edu/modis/level3/husf/gcoos/2009/037/1km/pass/final/MODIS.2009037.191026.gcoos.seadas_sst.hdf

    The SeaDAS file variables are:

                    senz:long_name = "Sensor Zenith Angle" ;
                    nLw_443:long_name = "Normalized water-leaving radiance at 443 nm" ;
                    aph_443_qaa:long_name = "Absorption due to phytoplankton at 443 nm, QAA algorithm" ;
                    flh:long_name = "Fluorescence Line Height" ;
                    seadas_K_490:long_name = "Diffuse attenuation coefficient at 490 nm (OBPG)" ;
                    chl_carder:long_name = "Chlorophyll Concentration, Carder model" ;
                    nLw_531:long_name = "Normalized water-leaving radiance at 531 nm" ;
                    nLw_412:long_name = "Normalized water-leaving radiance at 412 nm" ;
                    taua_551:long_name = "Aerosol optical thickness at 551 nm" ;
                    nLw_488:long_name = "Normalized water-leaving radiance at 488 nm" ;
                    taua_869:long_name = "Aerosol optical thickness at 869 nm" ;
                    Kd_488_lee:long_name = "Diffuse Attenuation Coefficient at 488 nm, Lee Algorithm" ;
                    chlor_a:long_name = "Chlorophyll Concentration, OC3 Algorithm" ;
                    nLw_869:long_name = "Normalized water-leaving radiance at 869 nm" ;
                    adg_443_carder:long_name = "Absorption due to gelbstof and detrital material at 443 nm, Carder algorithm" ;
                    nLw_748:long_name = "Normalized water-leaving radiance at 748 nm" ;
                    adg_443_qaa:long_name = "Absorption due to gelbstof and detrital material at 443 nm, QAA algorithm" ;
                    nLw_551:long_name = "Normalized water-leaving radiance at 551 nm" ;
                    cfe:long_name = "Fluorescence Efficiency" ;
                    sena:long_name = "Sensor Azimuth Angle" ;
                    nLw_678:long_name = "Normalized water-leaving radiance at 678 nm" ;
                    angstrom_551:long_name = "Angstrom coefficient, 551 to 869 nm" ;
                    nLw_667:long_name = "Normalized water-leaving radiance at 667 nm" ;
                    tsm_clark:long_name = "Total Suspended Matter, D. Clark, 2003" ;
                    bbp_551_qaa:long_name = "Particulate backscatter at 551 nm, QAA algorithm" ;
                    eps_78:long_name = "Epsilon of aerosol correction at 765 and 865 nm" ;
    

    The global attributes are:

    // global attributes:
                    :Map Projection Category = "Proj 4.4.7" ;
                    :Data Source Type = "Direct Broadcast" ;
                    :L3 Software Name = "/usr/bin/map_modis.pl" ;
                    :L3 Software Version = "4.2" ;
                    :Map Limit = "18, -98, 31, -79" ;
                    :Processing Status = "final" ;
                    :Collection Location = "USF Institute for Marine Remote Sensing" ;
                    :Satellite = "Aqua" ;
    		:Input File(s) = "..." ;
                    :Map Projection Name = "eqc" ;
                    :Proj Projection Parameters = "#Equidistant Cylindrical (Plate Caree)\n",
        "#\tCyl, Sph\n",
        "#\tlat_ts=\n",
        "# +proj=eqc +ellps=WGS84" ;
                    :Temporal Resolution = "pass" ;
                    :Spatial Resolution = "1000" ;
    

    The SeaDAS SST HDF file looks like this:

    netcdf MODIS.2009037.191026.gcoos.seadas_sst {
    dimensions:
            fakeDim0 = 1430 ;
            fakeDim1 = 2090 ;
            fakeDim2 = 1430 ;
            fakeDim3 = 2090 ;
            fakeDim4 = 1430 ;
            fakeDim5 = 2090 ;
    
    variables:
            short seadas_sst(fakeDim0, fakeDim1) ;
                    seadas_sst:long_name = "Sea Surface Temperature" ;
                    seadas_sst:slope = "[  0.005]" ;
                    seadas_sst:units = "degrees-C" ;
                    seadas_sst:intercept = "[0]" ;
            long l2_flags(fakeDim2, fakeDim3) ;
                    l2_flags:f02_name = "LAND" ;
                    l2_flags:long_name = "Level-2 Processing Flags" ;
                    l2_flags:f32_name = "OCEAN" ;
                    l2_flags:f03_name = "BADANC" ;
                    l2_flags:f21_name = "MODGLINT" ;
                    l2_flags:f10_name = "CLDICE" ;
                    l2_flags:f05_name = "HILT" ;
                    l2_flags:f06_name = "HISATZEN" ;
                    l2_flags:f07_name = "COASTZ" ;
                    l2_flags:f08_name = "NEGLW" ;
                    l2_flags:f25_name = "SEAICE" ;
                    l2_flags:f20_name = "MAXAERITER" ;
                    l2_flags:f13_name = "HISOLZEN" ;
                    l2_flags:f31_name = "SPARE" ;
                    l2_flags:f19_name = "TRICHO" ;
                    l2_flags:f18_name = "ABSAER" ;
                    l2_flags:f30_name = "HIPOL" ;
                    l2_flags:f04_name = "HIGLINT" ;
                    l2_flags:f17_name = "NAVWARN" ;
                    l2_flags:f01_name = "ATMFAIL" ;
                    l2_flags:f26_name = "NAVFAIL" ;
                    l2_flags:f24_name = "DARKPIXEL" ;
                    l2_flags:f16_name = "CHLFAIL" ;
                    l2_flags:f28_name = "SSTWARN" ;
                    l2_flags:f15_name = "LOWLW" ;
                    l2_flags:f11_name = "COCCOLITH" ;
                    l2_flags:f12_name = "TURBIDW" ;
                    l2_flags:f23_name = "ATMWARN" ;
                    l2_flags:f09_name = "STRAYLIGHT" ;
                    l2_flags:f22_name = "CHLWARN" ;
                    l2_flags:f29_name = "SSTFAIL" ;
                    l2_flags:f27_name = "FILTER" ;
                    l2_flags:f14_name = "HITAU" ;
                    l2_flags:units = "dimensionless" ;
            byte cloud_mask(fakeDim4, fakeDim5) ;
    
    // global attributes:
                    :Map Projection Category = "Proj 4.4.7" ;
                    :Data Source Type = "Direct Broadcast" ;
                    :L3 Software Name = "/usr/bin/map_modis.pl" ;
                    :L3 Software Version = "4.2" ;
                    :Map Limit = "18, -98, 31, -79" ;
                    :Processing Status = "final" ;
                    :Collection Location = "USF Institute for Marine Remote Sensing" ;
                    :Satellite = "Aqua" ;
                    :Input File(s) = "..."
                    :Map Projection Name = "eqc" ;
                    :Proj Projection Parameters = "#Equidistant Cylindrical (Plate Caree)\n",
        "#\tCyl, Sph\n",
        "#\tlat_ts=\n",
        "# +proj=eqc +ellps=WGS84" ;
                    :Temporal Resolution = "pass" ;
                    :Spatial Resolution = "1000" ;
    
    data:
    
     seadas_sst =
      -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999,
        -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999,
    ...
        4270, 3538, 3538, 3538, 3623, 3623, 4716, 4716 ;
    
     l2_flags =
      4294934546, 4294934546, 4294934546, 4294934546, 4294934546, 4294934546,
        4294934546, 4294934546, 4294934546, 4294934546, 4294934546, 4294934546,
    ...
        4294935056, 4294935056 ;
    
     cloud_mask =
      0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    ...
        0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;
    }
    

    Sample files from LSU/ESL

    Tricky.


    Matt,

    It's the similar situation we had at CoastWatch, and Rich gave us great help to solve the problem.

    I have a CGI-based OPeNDAP server that serves the CoastWatch data in HDF4 format, e.g.,

    http://coastwatch.chesapeakebay.noaa.gov/dap/edac/modis_chlora/data/

    And then refer to the OPeNDAP URLs when construct the THREDDS catalog, e.g.

    .....
            <aggregation dimName="time" type="joinNew">
              <variableAgg name="chlor_a"/>
    <netcdf location="http://coastwatch.chesapeakebay.noaa.gov/dap/edac/modis_chlora/data/2009_164_1700_aqua_ex1_chlora.hdf"/>
    <netcdf location="http://coastwatch.chesapeakebay.noaa.gov/dap/edac/modis_chlora/data/2009_165_1600_aqua_ex1_chlora.hdf"/>
    <netcdf location="http://coastwatch.chesapeakebay.noaa.gov/dap/edac/modis_chlora/data/2009_166_1640_aqua_ex1_chlora.hdf"/>
    .....
    

    In this way you don't need to convert the HDF4 data to NetCDF format. If you need further explanation about what I did in more detail, please let me know. I can also provide some sample code if you need. Hope this helps.

    Xiaoming


    Matt & Xiaoming,

    With Thredds Data Server v4, you can use NcML on the HDF4 files directly (you don't need another opendap server). Here's an example catalog that Phil Keegstra and I worked on. The HDF4 file it points to is

    http://coast-enviro.er.usgs.gov/models/share/MODSCW_P2009168_C4_1805_1810_1940_1945_GM03_closest_chlora.hdf

    And if you execute this WCS command you will get back a subsetted NetCDF file:

    curl -o subset.nc 'http://blackburn.whoi.edu:8081/thredds/wcs/coastwatch/albers3?service=WCS&version=1.0&request=GetCoverage&format=NetCDF3&coverage=chlor_a&bbox=-92.7,28,-88,30.'

    Rich


    Notes:

    Local version of the HDF and extracted NetCDF subset files are here:

    MODSCW_P2009168_C4_1805_1810_1940_1945_GM03_closest_chlora.hdf
    subset.nc

    And if we choose the snag a GeoTIFF file instead with:

    curl -o subset.tiff 'http://blackburn.whoi.edu:8081/thredds/wcs/coastwatch/albers3?service=WCS&version=1.0&request=GetCoverage&format=GeoTIFF&coverage=chlor_a&bbox=-92.7,28,-88,30.'

    we get the following TIFF file:

    subset.tiff

    If we put the HDF file on our local server, we should be able to fetch it with a similar command, i.e.:

    curl -o subset2.nc 'http://megara.tamu.edu:8080/thredds/wcs/SeaDAS_SST/MODSCW.hdf?service=WCS&version=1.0&request=GetCoverage&format=NetCDF3&coverage=chlor_a&bbox=-92.7,28,-88,30.'

    The dimensions and variables from the HDF file are:

    File Detail:

    netcdf MODSCW_P2009168_C4_1805_1810_1940_1945_GM03_closest_chlora {
    dimensions:
            rows = 1534 ;
            cols = 2122 ;
    
    variables:
            float chlor_a(rows, cols) ;
                    chlor_a:long_name = "Chlorophyll Concentration, OC3 Algorithm" ;
                    chlor_a:units = "mg m^-3" ;
                    chlor_a:coordsys = "Albers Conical Equal Area" ;
                    chlor_a:_FillValue = -1.f ;
                    chlor_a:missing_value = -1.f ;
                    chlor_a:scale_factor = 1. ;
                    chlor_a:scale_factor_err = 0. ;
                    chlor_a:add_offset = -0. ;
                    chlor_a:add_offset_err = 0. ;
                    chlor_a:calibrated_nt = 0 ;
                    chlor_a:fraction_digits = 4 ;
            float latitude(rows, cols) ;
                    latitude:long_name = "latitude" ;
                    latitude:units = "degrees" ;
                    latitude:coordsys = "Albers Conical Equal Area" ;
                    latitude:_FillValue = -999.f ;
                    latitude:missing_value = -999.f ;
                    latitude:scale_factor = 1. ;
                    latitude:scale_factor_err = 0. ;
                    latitude:add_offset = -0. ;
                    latitude:add_offset_err = 0. ;
                    latitude:calibrated_nt = 0 ;
                    latitude:fraction_digits = 4 ;
            float longitude(rows, cols) ;
                    longitude:long_name = "longitude" ;
                    longitude:units = "degrees" ;
                    longitude:coordsys = "Albers Conical Equal Area" ;
                    longitude:_FillValue = -999.f ;
                    longitude:missing_value = -999.f ;
                    longitude:scale_factor = 1. ;
                    longitude:scale_factor_err = 0. ;
                    longitude:add_offset = -0. ;
                    longitude:add_offset_err = 0. ;
                    longitude:calibrated_nt = 0 ;
                    longitude:fraction_digits = 4 ;
            byte graphics(rows, cols) ;
                    graphics:long_name = "graphics overlay planes" ;
                    graphics:coordsys = "Albers Conical Equal Area" ;
                    graphics:_FillValue = '\0' ;
                    graphics:missing_value = '\0' ;
                    graphics:scale_factor = 1. ;
                    graphics:scale_factor_err = 0. ;
                    graphics:add_offset = 0. ;
                    graphics:add_offset_err = 0. ;
                    graphics:calibrated_nt = 0 ;
                    graphics:fraction_digits = 0 ;
    

    with the global attributes:

    // global attributes:
                    :satellite = "Aqua" ;
                    :sensor = "MODIS" ;
                    :origin = "USDOC/NOAA/NESDIS CoastWatch" ;
                    :history = "PGE01:5.0.29;PGE ..." ;
                    :cwhdf_version = "3.4" ;
                    :pass_date = 14412, 14412, 14412, 14412 ;
                    :start_time = 70808., 65408., 65108., 71108. ;
                    :temporal_extent = 299., 298., 299., 299. ;
                    :projection_type = "mapped" ;
                    :projection = "Albers Conical Equal Area" ;
                    :gctp_sys = 3 ;
                    :gctp_zone = 0 ;
                    :gctp_parm = 0., 0., 20030000., 27030000., -89000000., 24000000., 0., 0., 0., 0., 0., 0., 0., 0., 0. ;
                    :gctp_datum = 12 ;
                    :et_affine = 0., -1008.748000000001, 1008.748000000001, 0., -1069777.254000001, 809520.2699999964 ;
                    :rows = 1534 ;
                    :cols = 2122 ;
                    :polygon_latitude = 30.9284265727064, ... ;
                    :polygon_longitude = -100.140439404871,... ;
                    :pass_type = "day" ;
                    :orbit_type = "ascending" ;
                    :raster_type = "RasterPixelIsArea" ;
                    :swath_sync_lines = 10, 10, 10, 10 ;
                    :composite = "true" ;
    

    and the dimensions and variables from the subsetted NetCDF file are:

    subset {
    dimensions:
            rows = 226 ;
            cols = 459 ;
    variables:
            double chlor_a(rows, cols) ;
                    chlor_a:long_name = "Chlorophyll Concentration, OC3 Algorithm" ;
                    chlor_a:units = "mg m^-3" ;
                    chlor_a:coordsys = "Albers Conical Equal Area" ;
                    chlor_a:scale_factor_err = 0. ;
                    chlor_a:add_offset_err = 0. ;
                    chlor_a:calibrated_nt = 0 ;
                    chlor_a:fraction_digits = 4 ;
                    chlor_a:grid_mapping = "myMapping" ;
                    chlor_a:coordinates = "yc xc " ;
            double yc(rows) ;
                    yc:units = "km" ;
                    yc:long_name = "y-coordinate of albers projection" ;
                    yc:standard_name = "projection_y_coordinate" ;
            double xc(cols) ;
                    xc:units = "km" ;
                    xc:long_name = "x-coordinate of albers projection" ;
                    xc:standard_name = "projection_x_coordinate" ;
            int myMapping ;
                    myMapping:grid_mapping_name = "albers_conical_equal_area" ;
                    myMapping:standard_parallel = 20.5, 27.5 ;
                    myMapping:longitude_of_central_meridian = -89. ;
                    myMapping:latitude_of_projection_origin = 24. ;
                    myMapping:units = "m" ;
                    myMapping:false_easting = 0. ;
                    myMapping:false_northing = 0. ;
                    myMapping:_CoordinateTransformType = "Projection" ;
                    myMapping:_CoordinateAxisTypes = "GeoX GeoY" ;
            double lat(rows, cols) ;
                    lat:units = "degrees_north" ;
                    lat:long_name = "latitude coordinate" ;
                    lat:standard_name = "latitude" ;
                    lat:_CoordinateAxisType = "Lat" ;
            double lon(rows, cols) ;
                    lon:units = "degrees_east" ;
                    lon:long_name = "longitude coordinate" ;
                    lon:standard_name = "longitude" ;
                    lon:_CoordinateAxisType = "Lon" ;
    


    The catalog.xml file used to serve the sample file perform the transformations follows, and is heavily annotated so the reason for and documentation of each part is readily available for perusal.

    albers.xml

    <?xml version="1.0" encoding="us-ascii"?>
    <catalog name="Coastwatch Tests" xmlns="http://www.unidata.ucar.edu/namespaces/thredds/InvCatalog/v1.0" xmlns:xlink="http://www.w3.org/1999/xlink">
        <service name="allServices" serviceType="Compound" base="">
            <service name="ncdods" serviceType="OPENDAP" base="/thredds/dodsC/"/>
            <service name="wcs" serviceType="WCS" base="/thredds/wcs/"/>
            <service name="ncss" serviceType="NetcdfSubset" base="/thredds/ncss/grid/"/>
            <service name="HTTPServer" serviceType="HTTPServer" base="/thredds/fileServer/"/>
        </service>
    
    According to Rich:
        <dataset name="Albers Test 1: with lon/lat and coordinates attribute" ID="coastwatch/albers1" serviceName="allServices" urlPath="coastwatch/albers1">
            <netcdf xmlns="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2" location="http://coast-enviro.er.usgs.gov/models/share/MODSCW_P2009168_C4_1805_1810_1940_1945_GM03_closest_chlora.hdf">
                <variable name="chlor_a">
                    <attribute name="coordinates" value="longitude latitude"/>
                </variable>
            </netcdf>
        </dataset>
    
         <dataset name="Albers Test 3: without lon/lat and coordinates attribute" ID="coastwatch/albers3" serviceName="allServices" urlPath="coastwatch/albers3">
    
    The destination NetCDF dataset is created from the source or referenced dataset using the capabilities of NcML. To give a view of the forest before we start on the trees, the XML below:
    • Opens a container into which the finished file will be placed.
    • Reads the variable chlor_a and all of its attributes from the source file, and assigns additional attributes to it.
    • Deletes the longitude and latitude coordinate variables.
    • Creates a new variable myMapping containing the projection information needed to construct new coordinate variables.
    • Creates new coordinate variables xc and yc for calculating destination file x and y distance equivalents of the longitudes and latitudes in the source file.
    Now we move on to lengthy and tedious explanations of the individual trees.

    The netcdf element represents a generic NetCDF dataset, i.e. a container for data that conforms to the NetCDF data model. The subelements used here, as defined in the NcML Schema, are:

            <netcdf xmlns="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2" location="http://coast-enviro.er.usgs.gov/models/share/MODSCW_P2009168_C4_1805_1810_1940_1945_GM03_closest_chlora.hdf" enhance="true">
    
    Here we assign the variable chlor_a the attributes grid_mapping and coordinates, with the accompanying variables defined further down the way. All of the chlor_a attributes in the source file are also read and included in the destination dataset.
                <variable name="chlor_a">
                    <attribute name="grid_mapping" value="myMapping"/>
                    <attribute name="coordinates" value="xc yc"/>
                </variable>
    
    Here we remove the longitude and latitude variables that were automatically read when the netcdf element was used to open the dataset.
                <remove type="variable" name="longitude"/>
                <remove type="variable" name="latitude"/>
    
    In this section we create a MyMapping variable containing the projection information as attributes. Detailed information about the required attributes can be found in the appropriate section in the CF 1.4 Conventions Handbook, i.e. at:

    http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#grid-mappings-and-projections

    The Grid Mappings Appendix of that document supplies all the needed information for each of the supported projections.

                <variable name="myMapping" shape="" type="int">
                    <attribute name="grid_mapping_name" value="albers_conical_equal_area"/>
                    <attribute name="standard_parallel" type="double" value="20.5  27.5"/>
                    <attribute name="longitude_of_central_meridian" type="double" value="-89.0"/>
                    <attribute name="latitude_of_projection_origin" type="double" value="24.0"/>
                    <attribute name="units" value="m"/>
                    <attribute name="false_easting" type="double" value="0.0"/>
                    <attribute name="false_northing" type="double" value="0.0"/>
                </variable>
    
    In this section we create xc and yc variables containing the x and y coordinate information that can be calculated using the information about the projection available in the source HDF file.

    The attributes immediately following the variable name are a part of the NcML Schema that were not contained within CDL, the predecessor to NcML. They are:

    • shape - A list of the names of the dimensions on which the variable depends. In this case, xc depends on the cols dimension, and yc depends on the rows dimension.
    • type - The datatype of the variable, selected from the list of available DataTypes.
    The attributes below the variable name are taken from the Attributes Appendix of the CF definition document. They are:
    • units - A required attribute for all variables that represent dimensional quantities. The Units section of the CF standard supplies the recommendations for and restrictions on these.
    • long_name - A descriptive name indicating a variable's content. This is optional but recommended for supplying a lengthier and more fully descriptive alternative to the required standard_name.
    • standard_name - A required standardized name for the variable. A Standard Name Table supplies a list of permissable standard names along with their associated required canonical units and a description when additional information is needed for clarification. The standard names projection_x_coordinate and projection_y_coordinate indicate that x and y are the projection coordinates.
    • axis - The Coordinate Types section of the CF standard supplies this attribute for identifying generic spatial coordinates, specifying that X and Y be used for horizontal coordinate variables.
    The optional values element is provided under the variable element in the NcML Schema. It can be used to specify the actual data values of a variable. This can be done by either providing the entire list or by specifying a starting value and an increment, with the latter route taken in this example for the variables xc and yc. In this case, the inimitable Rich Signell carefully scrutinized the source HDF file and figured out that the et_affine global attribute was in fact a specification of the x and y limits of the data domain. The attribute string is:

    et_affine = 0., -1008.748000000001, 1008.748000000001, 0., -1069777.254000001, 809520.2699999964

    wherein the second and third values are the increments and the fifth and sixth the starting values. For a reason I can't currently fathom, the increments were added to each of the starting values to obtain the starting values actually seen below.

                <variable name="xc" shape="cols" type="double">
                    <attribute name="units" value="m"/>
                    <attribute name="long_name" value="x-coordinate of albers projection"/>
                    <attribute name="standard_name" value="projection_x_coordinate"/>
                    <attribute name="axis" value="X"/>
                    <!-- X values from the "et_affine" attribute in the HDF4 file:
                        start=et_affine(5)   increment=et_affine(2)  -->
                    <values start="-1070786" increment="1008.748"/>
                </variable>
    
                <variable name="yc" shape="rows" type="double">
                    <attribute name="units" value="m"/>
                    <attribute name="long_name" value="y-coordinate of albers projection"/>
                    <attribute name="standard_name" value="projection_y_coordinate"/>
                    <attribute name="axis" value="Y"/>
                    <!-- Y values from the "et_affine" attribute in the HDF4 file:
                        start=et_affine(6)   increment=et_affine(3)  -->
                     <values start="810529.0" increment="-1008.748"/>
                </variable>
                <attribute name="Conventions" value="CF-1.0"/>
            </netcdf>
        </dataset>
    
    </catalog>
    


    Rich,

    Thanks for the good point. Phil mentioned the THREDDS v4 release to me the other day, and I haven't got into detail yet. For Matt, I think it's now even easier to build a THREDDS/WCS service with existing HDF4 data. In the sample code, I noticed that the HDF4 data is still referred as a http link. The HDF4 data still need to be served through a web server, so that it can be referred in NcML? (if the data are on the same computer as the THREDDS server, can you can use a local link to a HDF4 file in NcML?)

    Xiaoming


    Xiaoming,

    Sure, you can use local HDF4 files. That is the preferred method.

    I only used a remote HDF4 file

    1. to show that it was possible -- pretty cool, no?
    2. so you could easily grab the HDF4 file from that web address if you wanted to

    Rich


    Matt,

    As Rich mentioned, you don't need another Opendap server if you use THREDDS v4. You can build a THREDDS catalog (add CF compliant meta data) directly using the local address of HDF4 files if they are on the same computer as the THREDDS server. Attached please find the shell scripts on linux that I used to generate a THREDDS catalog. But with THREDDS v4, you can modify the code to refer to the HDF4 data as a local address instead, e.g.:

    ...
    <netcdf location="/home/xiaoming/modis_chlora/data/2009_165_1600_aqua_ex1_chlora.hdf"/>
    ...
    

    Xiaoming

    gen_catalog.sh
    footer_avhrr.xml
    footer_goes.xml
    footer_modis_chla.xml
    footer_modis_k490.xml
    footer_modis_rrs667.xml
    header_avhrr.xml
    header_goes.xml
    header_modis_chla.xml
    header_modis_k490.xml
    header_modis_rrs667.xml
    header.xml
    variableAgg_avhrr.xml
    variableAgg_goes.xml
    variableAgg_modis_chla.xml
    variableAgg_modis_k490.xml
    variableAgg_modis_rrs667.xml
    variableAgg.xml


    All,

    FWIW, I have a test case on our public THREDDS server demonstrating this. (The other part is that the catalog is generated by a perl script, so it is suitable for actual use.) The announcement is one of the things on my "to-do" list. I'm curious to know if this is "good enough" to meet people's needs so that we can indeed dispense with the cwhdf2netcdf converter.

    http://coastwatch.noaa.gov/thredds/catalog_hdf4_test.html

    n.b. we still haven't solved the config issues with content-length in HEAD, so THREDDS-to-THREDDS chaining won't work for now.

    PBK


    All,

    OK, I think the perl code I described is "the rest of the story" you're asking about.

    It, along with a C language metadata querier binary which pulls the metadata out of the CWHDF, parses out the metadata in our CWHDF format and writes out NCML representing the extended NetCDF/CF format THREDDS is expecting.

    PBK


    Phil,

    Could you please pass along the perl script? I'm curious what you needed to do and it might serve as a starting point for similar scripts we might have to write down the line.

    It would also be interesting to see what more would need to be added to NcML to obviate the need for such scripts.

    I just found out recently that NcML in NJ4 has the capability to promote global attributes to variables, which might solve some of the issues that you are solving via the perl script, I imagine. Check out the last example on this page:

    http://www.unidata.ucar.edu/software/netcdf/ncml/v2.2/Cookbook.html

    Congrats on the getting the WCS w/ HDF4 via TDS4 into production!

    Rich


    Rich,

    What the perl does a lot more of is mapping attributes from one semantics to another.

    Perl script is attached.

    It's not production yet, just a demo to get feedback. If noone misses the few things we can do with real NetCDF, we can put it into production.

    PBK