[GRASS-user] netCDF import via r.in.gdal

[Hamish]

Luigi wrote:
> GRASS complains about bands -- the file includes georeferenced
> daily climate variables. The file opens OK in viewers such as
> Panoply <http://www.giss.nasa.gov/tools/panoply/>.

Hamish:
> > what does the output of gdalinfo look like?

> This is what gdalinfo says:
> 
> C:\Documents and
> Settings\Luigi\Desktop>gdalinfo
> C:\giotto\giotto_data\SRF_195
> 
> AVG.NC
> 
> Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
> Driver: netCDF/Network Common Data Format
> Files: C:\giotto\giotto_data\SRF_1958.AVG.NC
> Size is 512, 512
> Coordinate System is `'

   ^^^ not georeferenced (in a way gdal knows about), but....

> Metadata:
>   NC_GLOBAL#domxmin=-1.262909e+001
>   NC_GLOBAL#domxmax=4.295537e+001
>   NC_GLOBAL#domymin=2.129965e+001
>   NC_GLOBAL#domymax=6.144880e+001
>   NC_GLOBAL#domzmin=1.050000e+003
>   NC_GLOBAL#domzmax=1.050000e+003

do those make sense as lat/lon bounds? (12.63W, 42.96E; 21.29N, 61.44N)

> Subdatasets:
> SUBDATASET_1_NAME=NETCDF:"C:\giotto\giotto_data\SRF_1958.AVG.NC":lon
>   SUBDATASET_1_DESC=[148x158] lon (32-bit floating-point)
> 
> SUBDATASET_2_NAME=NETCDF:"C:\giotto\giotto_data\SRF_1958.AVG.NC":lat
>   SUBDATASET_2_DESC=[148x158] lat (32-bit floating-point)

here in the first two subdatasets you have lat and lon maps.

import them, then use the values in them to set grass map bounds with
r.region, instead of the starting 0,512 x 0x512 pixel coords.

if you are good at shell scripting you could automate the import of
all the bands and post the script to the NetCDF wiki page for others
who find similar data files.. :)


> SUBDATASET_3_NAME=NETCDF:"C:\giotto\giotto_data\SRF_1958.AVG.NC":ua
> SUBDATASET_4_NAME=NETCDF:"C:\giotto\giotto_data\SRF_1958.AVG.NC":va
> SUBDATASET_5_NAME=NETCDF:"C:\giotto\giotto_data\SRF_1958.AVG.NC":drag

.... I assume these move into the actual data (wind velocities?)

Hamish:
> you might try to use gdal_translate extract a single band from
> the file and convert it into a GeoTiff. Then import the geotiff
> into grass. But anything gdal_translate can read r.in.gdal
> should be able to read too (with the right options).

Luigi: 
> Found out that it is a multi-dataset file:
> Got info on subdatasets following directions at
> http://www.gdal.org/frmt_netcdf.html
....
> There is one subdataset per climate variable. Each subdataset has one
> band per day (365 bands). So I translated to geotiff as you suggested.
> 
> C:\giotto\giotto_data>gdal_translate -of GTiff -b 1
>     NETCDF:"SRF_1958.AVG.NC":tamax SRF_1958.AVG.tiff
>
> I then tried to import in grass the geotiff. It imports and display
> fine when creating a new location, but no projection info is carried
> along so no way to use it with other layers/locations in my grass
> dataset.

good to see the data makes it in.

so the process would look like:
0. gdalinfo to find band names
1. r.in.gdal with magic netcdf input string with filename and data name
   (or gdal_translate to geotiff + r.in.gdal)
2. examine lat and lon subdatasets; instead of importing perhaps try
   to use 'gdalinfo -stats' to read max/min of lat,lon layers.
3. apply georeferencing info with r.region or gdal_translate for
   geotiffs using the -a_srs and -a_ullr options.

> (by the way r.in.gdal in wxpython gui does not work
> grass6.4.0-rc3: I used the tcltk gui).

(look/file in the bug tracker please)


> Imports fine but still unable to get projection info, so the data is
> useless. If Panoply <http://www.giss.nasa.gov/tools/panoply/>
> is able to plot the data on a global coastline, the proj info
> should be somewhere in the netCDF file:

is there a special name for this flavour of netcdf file?
ie is it standardized way of geocoding them or some single nasa in-house
method not seen anywhere else?


Hamish