[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