[gdal-dev] Writing multi-bands vs multi-datasets

Joaquim Manuel Freire Luís jluis at ualg.pt
Fri Oct 15 06:59:23 PDT 2021


Even, sorry to keep bothering you with this but found no docs on this.

The CreateCopy() solution created me a file with the coordinates system info stored in the "transverse_mercator" attribute (why?)

Metadata:
  Band1#grid_mapping=transverse_mercator
  Band1#long_name=GDAL Band Number 1
  Band1#_FillValue=65535
  NC_GLOBAL#Conventions=CF-1.5
  NC_GLOBAL#GDAL=GDAL 3.4.0dev, released 2021/99/99
  NC_GLOBAL#history=Thu Oct 14 20:59:24 2021: GDAL CreateCopy( gd_cube.nc, ... )
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={7,6}
  NETCDF_DIM_time_VALUES={1,2,3,4,5,6,7}
  time#axis=t
  transverse_mercator#false_easting=500000
  transverse_mercator#false_northing=0
  transverse_mercator#GeoTransform=481845 30 0 4287765 0 -30
  transverse_mercator#grid_mapping_name=transverse_mercator
  transverse_mercator#inverse_flattening=298.257223563
  transverse_mercator#latitude_of_projection_origin=0
  transverse_mercator#longitude_of_central_meridian=-9
  transverse_mercator#longitude_of_prime_meridian=0
  transverse_mercator#long_name=CRS definition
  transverse_mercator#scale_factor_at_central_meridian=0.9996
  transverse_mercator#semi_major_axis=6378137
  transverse_mercator#spatial_ref=PROJCS["unknown",GEOGCS["unknown",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]


The problem is that either GMT and Mirone look for that info in global attribute called "grid_mapping" as in

...
Metadata:
  grid_mapping#spatial_ref=PROJCS["unknown",
    GEOGCS["unknown",
        DATUM["Unknown based on WGS84 ellipsoid",
            SPHEROID["WGS 84",6378137,298.257223563,


How can I control the name of the attribute where the CRS is tored so I can find it from other programs?

Also, is there a way of setting the "actual_range" and "_Fillvalue" data attributes?

Joaquim

From: Joaquim Manuel Freire Luís <jluis at ualg.pt>
Sent: Thursday, October 14, 2021 9:03 PM
To: Joaquim Manuel Freire Luís <jluis at ualg.pt>; Even Rouault <even.rouault at spatialys.com>; gdal-dev at lists.osgeo.org
Subject: RE: [gdal-dev] Writing multi-bands vs multi-datasets

On more. The warning about "No 1D variable is indexed by dimension time" turned out to be because I was calling

Gdal.setmetadataitem(ds, "NETCDF_DIM_EXTRA", "{time}")

(like you call '{time,Z}' in your Python example). Removing the curly braces here made that warning go away.

From: gdal-dev <gdal-dev-bounces at lists.osgeo.org<mailto:gdal-dev-bounces at lists.osgeo.org>> On Behalf Of Joaquim Manuel Freire Luís
Sent: Thursday, October 14, 2021 8:42 PM
To: Even Rouault <even.rouault at spatialys.com<mailto:even.rouault at spatialys.com>>; gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>
Subject: Re: [gdal-dev] Writing multi-bands vs multi-datasets

Much better (I was confused by those curly braces thinking it was part of python syntax).
Though still get this warning

Warning 1: No 1D variable is indexed by dimension time

But otherwise gdalinfo does not complain on anything wrong

The "-co" is confusing. I had to add in other situations otherwise some options were not set.

Julia does not use bindings, it can access the exported functions in shared libraries directly
(if you are curious see https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/
and https://docs.julialang.org/en/v1/base/c/#ccall)
so passing that vector of strings should be translated automatically to a "char *string[]" and GDAL expects the "-co" (I think)

From: Even Rouault <even.rouault at spatialys.com<mailto:even.rouault at spatialys.com>>
Sent: Thursday, October 14, 2021 8:14 PM
To: Joaquim Manuel Freire Luís <jluis at ualg.pt<mailto:jluis at ualg.pt>>; gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>
Subject: Re: [gdal-dev] Writing multi-bands vs multi-datasets



Le 14/10/2021 à 21:03, Joaquim Manuel Freire Luís a écrit :
Even,

There are some issues with this way.

My draft Julia function looks like bellow  (tested with a MxNx7 array) but get this errors

julia> GMT.write_multiband(cube, "gd_cube.nc")
Warning 1: No 1D variable is indexed by dimension time
Warning 6: creation option '-co' is not formatted with the key=value format
Warning 6: creation option '-co' is not formatted with the key=value format
Warning 6: creation option '-co' is not formatted with the key=value format
I don't know the Julia bindings, but the error message seems to suggest you can remove the "-co" element in your options[] array since it presumably must only contain the "key=value" creation options, like in the Python bindings
ERROR 1: netcdf error #-45 : NetCDF: Not a valid data type or _FillValue type mismatch .
at (C:\programs\compa_libs\gdal_GIT\gdal\frmts\netcdf\netcdfdataset.cpp,netCDFDataset::CreateCopy,9012)

ERROR 1: netcdf error #-49 : NetCDF: Variable not found .
at (C:\programs\compa_libs\gdal_GIT\gdal\frmts\netcdf\netcdfdataset.cpp,NCDFPutAttr,10454)

Warning 1: No 1D variable is indexed by dimension time

And indeed the "time" variable shows up filled with 7 zeros when I open the nc file that is nevertheless created.
Also checked with getmetadataitem(ds, "NETCDF_DIM_time_VALUES")  that the "1,2,3,4,5,6,7" info was correctly set.
You need to surround the values of the _DEF and _VALUES items with { and } braces to indicate this is a list (this is the particular syntax expected by the netCDF driver)

I's also strange the warnings about the "-co" options also because they were in fact taken into account. Omitting them creates a classic nc file.


function write_multiband(cube, fname::AbstractString, v=nothing; dim_name::String="time")
               nbands = size(cube,3)
               ds = gmt2gd(cube)
               Gdal.setmetadataitem(ds, "NETCDF_DIM_EXTRA", dim_name)
               Gdal.setmetadataitem(ds, "NETCDF_DIM_" * dim_name * "_DEF", "$(nbands),6")
               Gdal.setmetadataitem(ds, "NETCDF_DIM_" * dim_name * "_VALUES", "1,2,3,4,5,6,7")
               Gdal.setmetadataitem(ds, dim_name * "#axis", string(dim_name[1]))
               crs = Gdal.getproj(cube, wkt=true)
               (crs != "" ) && Gdal.setproj!(ds, crs)
               Gdal.unsafe_copy(ds, filename=fname, driver=getdriver("netCDF"), options=["-co", "FORMAT=NC4", "-co", "COMPRESS=DEFLATE", "-co", "ZLEVEL=4"])
end

From: Even Rouault <even.rouault at spatialys.com><mailto:even.rouault at spatialys.com>
Sent: Sunday, October 10, 2021 10:48 PM
To: Joaquim Manuel Freire Luís <jluis at ualg.pt><mailto:jluis at ualg.pt>; gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>
Subject: Re: [gdal-dev] Writing multi-bands vs multi-datasets


Joaquim,

https://github.com/OSGeo/gdal/pull/4634 should hopefully help

Even
Le 10/10/2021 à 22:26, Joaquim Manuel Freire Luís a écrit :
Hi,

I have a Julia version of gdal_translate (that ends up calling GDALDatasetRasterIOEx) that works in saving a MxNxP array, but not the way I wanted.
It creates a multi-subdatasets file, whilst what I wished is to write a multi-band file.

For example, the file written with this Julia function looks like below

How can I write a multi-band instead of a multi-datasets?
I've been looking at the docs of GDALDatasetRasterIOEx but can't see any way of picking one or the other.

Joaquim

gdalinfo cube_rad.nc
Driver: netCDF/Network Common Data Format
Files: cube_rad.nc
Size is 512, 512
Metadata:
  NC_GLOBAL#Conventions=CF-1.5
  NC_GLOBAL#GDAL=GDAL 3.4.0dev, released 2021/99/99
  NC_GLOBAL#history=Sun Oct 10 13:01:06 2021: GDAL CreateCopy( cube_rad.nc, ... )
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"cube_rad.nc":Band1
  SUBDATASET_1_DESC=[1567x1519] Band1 (32-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:"cube_rad.nc":Band2
  SUBDATASET_2_DESC=[1567x1519] Band2 (32-bit floating-point)
...

Saving the same data (or similar) from a Matlab code that I wrote long ago to stack grids, shows

gdalinfo evi_cube.nc
Driver: netCDF/Network Common Data Format
Files: evi_cube.nc
Size is 37, 59
Origin = (580140.000000000000000,4418550.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
...
  latitude#actual_range={4416795,4418535}
  latitude#long_name=latitude
  latitude#units=degrees_north
  longitude#actual_range={580155,581235}
  longitude#long_name=longitude
  longitude#units=degrees_east
  NC_GLOBAL#Conventions=COARDS
  NC_GLOBAL#description=File written from Matlab
  NC_GLOBAL#node_offset=0
  NC_GLOBAL#Source_Software=Mirone
  NC_GLOBAL#title=Grid created by Mirone
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={67,6}
  NETCDF_DIM_time_VALUES={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67}
  z#actual_range={-856.4703369140625,176.8457641601562}
  z#long_name=z
  z#units=unknown
  z#_FillValue=-nan(ind)
Corner Coordinates:
Upper Left  (  580140.000, 4418550.000)
Lower Left  (  580140.000, 4416780.000)
Upper Right (  581250.000, 4418550.000)
Lower Right (  581250.000, 4416780.000)
Center      (  580695.000, 4417665.000)
Band 1 Block=37x59 Type=Float32, ColorInterp=Undefined
  NoData Value=nan
  Unit Type: unknown
  Metadata:
    actual_range={-856.4703369140625,176.8457641601562}
    long_name=z
    NETCDF_DIM_time=1
    NETCDF_VARNAME=z
    units=unknown
    _FillValue=-nan(ind)
Band 2 Block=37x59 Type=Float32, ColorInterp=Undefined
  NoData Value=nan
  Unit Type: unknown
  Metadata:
    actual_range={-856.4703369140625,176.8457641601562}
    long_name=z
    NETCDF_DIM_time=2
    NETCDF_VARNAME=z
    units=unknown
    _FillValue=-nan(ind)
...


_______________________________________________

gdal-dev mailing list

gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/gdal-dev

--

http://www.spatialys.com

My software is free, but my time generally not.

--

http://www.spatialys.com

My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20211015/aa3d8d90/attachment-0001.html>


More information about the gdal-dev mailing list