[gdal-dev] Accessing to a NetCDF subdataset with Python

[Paolo Corti]

On Sat, May 25, 2013 at 7:11 PM, Even Rouault
<even.rouault at mines-paris.org> wrote:
> Le samedi 25 mai 2013 17:47:59, Paolo Corti a écrit :
>> Hi list
>>
>> I am with GDAL 1.10.
>> I thought it was possible to access to a netCDF subdataset in the
>> usual manner (using GetRasterBand), but it does not seem to work in my
>>
>> local box:
>> >>> from osgeo import gdal
>> >>> from osgeo.gdalconst import *
>> >>> gdal.AllRegister()
>> >>> file_nc = 'NETCDF:"sresa1b_ncar_ccsm3_0_run1_200001.nc"'
>> >>> ds = gdal.Open(file_nc, GA_ReadOnly)
>> >>> ds.RasterCount
>>
>> 0
>
> If the netcdf file has only one variable with 2 dimensions, then you can
> directly open the dataset with the filename. If there are more than one
> variable with 2 dimensions, you need to specify the subdataset name
> (logically)
>
>>
>> while, using the subdataset name access syntax, it is working:
>> >>> from osgeo import gdal
>> >>> from osgeo.gdalconst import *
>> >>> gdal.AllRegister()
>> >>> file_nc = 'NETCDF:"sresa1b_ncar_ccsm3_0_run1_200001.nc":area'
>> >>> ds = gdal.Open(file_nc, GA_ReadOnly)
>>
>> Warning 1: Latitude grid not spaced evenly.
>> Seting projection for grid spacing is within 0.1 degrees threshold.
>>
>> >>> ds.RasterCount
>>
>> 1
>>
>> Now, do I have a local problem, or the second way is the only way to
>> access to a subdataset?
>
> See above answer
>
>> And in such a case, is the only way to programmatically access to the
>> subdatasets by reading their name using GetSubdataset or there is a
>> better way?
>
> I can't think of a better way.
>
>>
>> >>> ds.GetSubDatasets()
>>
>> [('NETCDF:"sresa1b_ncar_ccsm3_0_run1_200001.nc":area', '[128x256] area
>> (32-bit floating-point)'),
>> ('NETCDF:"sresa1b_ncar_ccsm3_0_run1_200001.nc":lat_bnds', '[128x2]
>> lat_bnds (64-bit floating-point)'),
>> ('NETCDF:"sresa1b_ncar_ccsm3_0_run1_200001.nc":lon_bnds', '[256x2]
>> lon_bnds (64-bit floating-point)'),
>> ('NETCDF:"sresa1b_ncar_ccsm3_0_run1_200001.nc":msk_rgn', '[128x256]
>> msk_rgn (32-bit integer)'),
>> ('NETCDF:"sresa1b_ncar_ccsm3_0_run1_200001.nc":pr', '[1x128x256]
>> precipitation_flux (32-bit floating-point)'),
>> ('NETCDF:"sresa1b_ncar_ccsm3_0_run1_200001.nc":tas', '[1x128x256]
>> air_temperature (32-bit floating-point)'),
>> ('NETCDF:"sresa1b_ncar_ccsm3_0_run1_200001.nc":time_bnds', '[1x2]
>> time_bnds (64-bit floating-point)'),
>> ('NETCDF:"sresa1b_ncar_ccsm3_0_run1_200001.nc":ua', '[1x17x128x256]
>> eastward_wind (32-bit floating-point)')]
>>
>> >>> sds = gdal.Open(ds.GetSubDatasets()[0][0])
>>
>> thanks in advance!
>> p
>
> --
> Geospatial professional services
> http://even.rouault.free.fr/services.html

Thanks Even!
I have not been working with dataset having subdatasets for a while
now, and I did
not consider that, while this is not the case for the netcdf file I am
working with, the subdatasets may be multibands!
Now you put me on the correct route, thanks :)
cheers
p

-- 
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti