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 <jl...@ualg.pt> Sent: Thursday, October 14, 2021 9:03 PM To: Joaquim Manuel Freire Luís <jl...@ualg.pt>; Even Rouault <even.roua...@spatialys.com>; gdal-dev@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-boun...@lists.osgeo.org<mailto:gdal-dev-boun...@lists.osgeo.org>> On Behalf Of Joaquim Manuel Freire Luís Sent: Thursday, October 14, 2021 8:42 PM To: Even Rouault <even.roua...@spatialys.com<mailto:even.roua...@spatialys.com>>; gdal-dev@lists.osgeo.org<mailto:gdal-dev@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.roua...@spatialys.com<mailto:even.roua...@spatialys.com>> Sent: Thursday, October 14, 2021 8:14 PM To: Joaquim Manuel Freire Luís <jl...@ualg.pt<mailto:jl...@ualg.pt>>; gdal-dev@lists.osgeo.org<mailto:gdal-dev@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.roua...@spatialys.com><mailto:even.roua...@spatialys.com> Sent: Sunday, October 10, 2021 10:48 PM To: Joaquim Manuel Freire Luís <jl...@ualg.pt><mailto:jl...@ualg.pt>; gdal-dev@lists.osgeo.org<mailto:gdal-dev@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@lists.osgeo.org<mailto:gdal-dev@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.
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev