[gdal-dev] GDAL multi bands files conversion to NetCDF

Fri Jul 20 02:54:49 PDT 2018

On mercredi 18 juillet 2018 18:57:18 CEST Rousseau Lambert2, Louis-Philippe 
(EC) wrote:
> Thanks for your answer Even,
> 
> I was looking at a one step solution because I wanted MapServer to
> output the NetCDF with different bands instead of individual
> Subdatasets... I should have mentioned that first in my question.
> 
> Basically, I have layers in MapServer which are served in WCS. Since we
> can't use TIME= in WCS 2.0.1 and that I need to be able to output time
> range and lists in a single WCS request, I'm using the band names and
> the RANGESUBSET= feature to request multiple bands in the NetCDF output
> format.
> 
> But the MapServer output NetCDF have many SUBDATASETs instead of many
> bands... I was hoping to be able to add a FORMATOPTION in the NetCDF
> OUTPUTFORMAT in MapServer to make it work, but from your answer I get
> that it's not possible.

You are seing more an artifact of how the GDAL netCDF driver than a limitation 
of MapServer or the netCDF format itself.

When creating a netCDF dataset with multiple GDAL bands, it get reads back by 
GDAL as several subdatasets, but neither the concept of band or subdataset 
really exist in the netCDF format itself: you have variables in netCDF. As 
variables may be indexed by different dimensions, GDAL chooses to read them as 
subdatasets to be on the safe side. Readers that are closer to the netCDF data 
model will see those bands/subdatasets as different variables.
The GDAL netCDF driver could potentially be improved to (as an option ?) offer 
the possibility to expose those subdatasets as bands when this is possible 
(same indexing dimensions)

Even

> 
> Thanks !
> 
> LP
> 
> On 07/18/2018 06:44 PM, Even Rouault wrote:
> > Louis Philippe
> > 
> > there's no one single step to do that, but you can do that in 2 steps:
> > 
> > - use gdal_translate -sds to export all the subdatasets in separate tif
> > files - use gdal_merge.py -separate to combine them in a single tif file
> > (provided they have same resolution, extent, cordinate system), which
> > seems to be your case here
> > 
> > Even
> > 
> >> Hi everyone!
> >> 
> >> I have s question regarding the conversion of a multiband file (like
> >> GeoTIFF) to NetCDF.
> >> 
> >> First my observations:
> >> 
> >> I used gdal_translate to convert a Geotiff to NetCDF:
> >> 
> >> gdal_translate -of input.tif result.nc
> >> 
> >> The gdalinfo result is:
> >> 
> >> Driver: netCDF/Network Common Data Format
> >> Files: result.nc
> >> Size is 512, 512
> >> Coordinate System is `'
> >> 
> >> Metadata:
> >>   NC_GLOBAL#Conventions=CF-1.5
> >>   NC_GLOBAL#GDAL=GDAL 1.10.0, released 2013/04/24
> >>   NC_GLOBAL#GDAL_AREA_OR_POINT=Area
> >>   NC_GLOBAL#history=Wed Jul 18 15:33:16 2018: GDAL CreateCopy( result.nc,
> >> 
> >> ... ) Subdatasets:
> >>   SUBDATASET_1_NAME=NETCDF:"result.nc":Band1
> >>   SUBDATASET_1_DESC=[1500x1650] Band1 (32-bit floating-point)
> >>   SUBDATASET_2_NAME=NETCDF:"result.nc":Band2
> >>   SUBDATASET_2_DESC=[1500x1650] Band2 (32-bit floating-point)
> >>   SUBDATASET_3_NAME=NETCDF:"result.nc":Band3
> >>   SUBDATASET_3_DESC=[1500x1650] Band3 (32-bit floating-point)
> >>   SUBDATASET_4_NAME=NETCDF:"result.nc":Band4
> >>   SUBDATASET_4_DESC=[1500x1650] Band4 (32-bit floating-point)
> >>   SUBDATASET_5_NAME=NETCDF:"result.nc":Band5
> >>   SUBDATASET_5_DESC=[1500x1650] Band5 (32-bit floating-point)
> >> 
> >> Corner Coordinates:
> >> Upper Left  (    0.0,    0.0)
> >> Lower Left  (    0.0,  512.0)
> >> Upper Right (  512.0,    0.0)
> >> Lower Right (  512.0,  512.0)
> >> Center      (  256.0,  256.0)
> >> 
> >> As you can see, the different bands are written in different SUBDATASET.
> >> I
> >> would like the output to be only different bands, not different datasets.
> >> I
> >> don't think I saw a creation option that could do that...
> >> 
> >> I would like something like:
> >> 
> >> Driver: netCDF/Network Common Data Format
> >> Files: output.nc
> >> Size is 1068, 510
> >> Coordinate System is `'
> >> Origin = (-140.999964377091317,83.499994903744323)
> >> Pixel Size = (0.083333295198278,-0.083333313348017)
> >> 
> >> Metadata:
> >>   lat#axis=Y
> >>   lat#long_name=latitude
> >>   lat#standard_name=latitude
> >>   lat#units=degrees_north
> >>   lon#axis=X
> >>   lon#long_name=longitude
> >>   lon#standard_name=longitude
> >>   lon#units=degrees_east
> >> 
> >> ...Some metadata...
> >> Corner Coordinates:
> >> Upper Left  (-140.9999644,  83.4999949)
> >> Lower Left  (-140.9999644,  41.0000051)
> >> Upper Right ( -52.0000051,  83.4999949)
> >> Lower Right ( -52.0000051,  41.0000051)
> >> Center      ( -96.4999847,  62.2500000)
> >> Band 1 Block=1068x1 Type=Float32, ColorInterp=Undefined
> >> 
> >>     Computed Min/Max=-29.841,11.720
> >>   
> >>   NoData Value=-99
> >>   
> >>   Metadata:
> >>     long_name=Autumn mean temperature
> >>     missing_value=-99
> >>     NETCDF_DIM_time=17093376
> >>     NETCDF_VARNAME=tm_son
> >>     _FillValue=-99
> >> 
> >> Band 2 Block=1068x1 Type=Float32, ColorInterp=Undefined
> >> 
> >>     Computed Min/Max=-30.198,11.576
> >>   
> >>   NoData Value=-99
> >>   
> >>   Metadata:
> >>     long_name=Autumn mean temperature
> >>     missing_value=-99
> >>     NETCDF_DIM_time=17102136
> >>     NETCDF_VARNAME=tm_son
> >>     _FillValue=-99
> >> 
> >> Band 3 Block=1068x1 Type=Float32, ColorInterp=Undefined
> >> 
> >>     Computed Min/Max=-30.479,11.854
> >>   
> >>   NoData Value=-99
> >>   
> >>   Metadata:
> >>     long_name=Autumn mean temperature
> >>     missing_value=-99
> >>     NETCDF_DIM_time=17110920
> >>     NETCDF_VARNAME=tm_son
> >>     _FillValue=-99
> >> 
> >> Is there anything I missed ? Can you think of another way to do that ?
> >> 
> >> Thanks a lot!
> >> 
> >> LP

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com