[gdal-dev] Writing multi-bands vs multi-datasets
Even Rouault
even.rouault at spatialys.com
Fri Oct 15 07:16:20 PDT 2021
Le 15/10/2021 à 15:59, Joaquim Manuel Freire Luís a écrit :
>
> Even, sorry to keep bothering you with this but found no docs on this.
>
Well, the ultimate doc is
https://github.com/OSGeo/gdal/blob/master/gdal/frmts/netcdf/netcdfdataset.cpp
>
> The CreateCopy() solution created me a file with the coordinates
> system info stored in the “transverse_mercator” attribute (why?)
>
why not ? :-) This is like this before I was born.
If you look at "Example 5.6. Rotated pole grid" of
https://cfconventions.org/Data/cf-conventions/cf-conventions-1.9/cf-conventions.html,
you can see that the name of the variable that stores the CRS definition
can be anything. The link to it is done through the "grid_mapping"
attribute of the variable of interest.
>
> How can I control the name of the attribute where the CRS is tored so
> I can find it from other programs?
>
You can't. It is set at
https://github.com/OSGeo/gdal/blob/40d18b06d78d5b97f9c6387ef6a08a9c0142c884/gdal/frmts/netcdf/netcdfdataset.cpp#L5607
with the same value as the #grid_mapping_name attribute of the CRS
definition vriable.
>
> Also, is there a way of setting the “actual_range” and “_Fillvalue”
> data attributes?
>
actual_range: perhaps by setting a "variable#actual_range={min,max}"
metadata item ?
_Fillvalue: through the SetNoDataValue() method of the raster band object
> 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/
> <https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/>
>
> and https://docs.julialang.org/en/v1/base/c/#ccall
> <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
> <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 <https://lists.osgeo.org/mailman/listinfo/gdal-dev>
>
> --
>
> http://www.spatialys.com <http://www.spatialys.com>
>
> My software is free, but my time generally not.
>
> --
> http://www.spatialys.com <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/46cae483/attachment-0001.html>
More information about the gdal-dev
mailing list