[gdal-dev] HowTo open NetCDF SUBDATASET and write to new file

[Steve Gaffigan]

Roger,

Hello.  Our examples differ in that you called the ReadAsArray() on a
single band, which returns an array with dimensions (ny,nx), and I called
it on the subdataset, which returns an array with dimensions
(nbands,ny,nx).  In my example the array returned has size (539,78,141),
where the two outer dimensions are concatenated into the single band
dimension.  The reshape() method is called on the array to get to the 4-d
form that is in the netcdf file.

Note, this works if gdal python bindings use numpy.  If bindings use
Numeric, then instead it should be either ds.ReadAsArray().resize(...) or
Numeric.reshape(ds.ReadAsArray(), ...).

I arrived at my example through the "banging away at it" method, too.

Steve


> Hi Steve,
>
> Thanks very much for the code snippet.  I figured out the following via
> experimentation.
>
> - Open specific layer by using:
> layer = gdal.Open('netCDF:"annual_mean.nc":temp')
>
> - Then get at the layer data data using the following:
> netcdf_band = layer.GetRasterBand(1)
>
> - Read in the data one row at a time:
>     for iY in range(y_size):
>       netcdf_data = netcdf_band.ReadAsArray(0, iY, x_size, 1)
>
> It works, but is a complete hack in the sense that I just plugged in stuff
> from a bunch of different scripts and banged away at it until it gave me
> what I needed.  Glad to see that it is similar to your code.  I'm curious
> about the ".reshape()" method that you pass to "ReadAsArray" above.  What
> does that do?
>
> Thanks again for the help.
>
> Roger
> --
>
> On Tue, Sep 16, 2008 at 12:05 AM, Steve Gaffigan
> <gaffigan at sfos.uaf.edu>wrote:
>
>> Hello.  Here's an example of reading an entire array for a single
>> subdataset.
>>
>> [gaffigan at dp2 AEFF]$ gdalinfo latest.nc
>> Driver: netCDF/Network Common Data Format
>> ...
>>  SUBDATASET_25_NAME=NETCDF:"latest.nc":TMP_GPML
>>  SUBDATASET_25_DESC=[49x11x78x141]
>> air_temperature_at_constant_altitude_above_m
>> ean_sea_level (8-bit integer)
>> ...
>>
>> [gaffigan at dp2 AEFF]$ python
>> >>> from osgeo import gdal
>> >>> ds = gdal.Open('NETCDF:"latest.nc":TMP_GPML')
>> >>> data = ds.ReadAsArray().reshape(49,11,78,141)
>>
>> Using a more generic call of the form below you can get at particular
>> hyperslabs:
>>
>> >>>
>> data=ds.GetRasterBand(i).ReadAsArray(xoff=x0,yoff=y0,win_xsize=nx,win_ysize=ny)
>>
>> where, for the above example, i=[1,...,539].
>>
>> Steve
>>
>>
>> > Hi List,
>> >
>> > I'm would like to use the Python bindings to open a specific
>> SUBDATASET
>> in
>> > a
>> > NetCDF file and read the contents.  I read in
>> > http://www.gdal.org/gdal_datamodel.html, under the "SUBDATASETS
>> domain"
>> > section, that I should be able to pass the "_NAME" parameter into
>> > GDALOpen()
>> > in order to do this, but I'm not exactly sure how.  Could someone pass
>> me
>> > a
>> > couple lines of Python that show how this is done?
>> >
>> > Thanks in advance.
>> >
>> > Roger
>> > --
>> > _______________________________________________
>> > gdal-dev mailing list
>> > gdal-dev at lists.osgeo.org
>> > http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>>
>>
>