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

Luigi Ponti lponti at infinito.it
Tue Apr 14 20:00:22 EDT 2009


On 14/04/2009 12.13, Hamish wrote:
> 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?
>>>       
>> [...]
>> 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)
>   

Yes, it is a bounding box around the Mediterranean Sea.

>   
>> 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.
>   

When I try to import to a latlong location, I get

    G_set_window(): Illegal latitude for North

and no raster is created. When I use a projected location, the file gets 
imported but then I have bounding coordinates in dd:mm:ss and r.region 
complains that those are invalid.

    r.region map=SRF_1958.AVG n=61.4488N s=21.29965N e=43.07446E w=12.74702W
    <n=61.4488N> ** illegal value **

Not sure what I am doing wrong now.

> 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.. :)
>   

"good at shell scripting" is a bit out of reach for me, but will do my best.

>> 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?)
>   

Yes, wind speed among other things: it is output from a regional climate 
model.

> 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.
>   

Stuck at point 3. Also tried the -a_ullr option with gdal_translate but 
imported GeoTiff does not show up. Maybe GRASS gets the coordinates in 
meters as opposed to degrees.

 +----------------------------------------------------------------------------+
 |                                                                            
|
 +----------------------------------------------------------------------------+
 | Layer:    SRF_1958.AVG                   Date: Wed Apr 15 01:17:47 
2009    |
 | Mapset:   luigi                          Login of Creator: 
Luigi           |
 | Location: 
EurLCC                                                           |
 | DataBase: 
C:/cygwin/home/andy                                              |
 | Title:     ( SRF_1958.AVG 
)                                                |
 | Timestamp: 
none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            
|
 |   Type of Map:  raster               Number of Categories: 
255             |
 |   Data Type:    
FCELL                                                      |
 |   Rows:         
148                                                        |
 |   Columns:      
158                                                        |
 |   Total Cells:  
23384                                                      |
 |        Projection: Lambert Conformal 
Conic                                 |
 |            N:    61.4488    S:   21.29965   Res: 
0.27127804                |
 |            E:   43.07446    W:  -12.74702   Res: 
0.35330051                |
 |   Range of data:    min = 258.227509  max = 
303.087402                     |
 |                                                                            
|
 |   Data 
Description:                                                        |
 |    generated by 
r.in.gdal                                                  |
 |                                                                            
|
 |   
Comments:                                                                |
 |    r.in.gdal -o input="C:/giotto/giotto_data/SRF_1958.AVG.tiff" 
output=\   |
 |    
"SRF_1958.AVG"                                                          |


> [...]
>> 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?
>   

Not sure that there is a special name, but digging around I found that 
the netCDF file has a companion file with .DDF extension which includes 
a PDEF string:

    PDEF  158  148 lcc   41.00   15.00   79.00   74.00   30.00   60.00  
    15.00  30000.  30000.

whose syntax is documented here 
http://grads.iges.org/grads/gadoc/pdef.html as

    /PDEF isize jsize LCC latref lonref iref jref Struelat Ntruelat slon
    dx dy/

which would suggest this data is pre-projected (no idea what it means) 
according to the above link.

Thanks and sorry if this thread is starting to go off topic,

Luigi
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20090415/85738080/attachment.html


More information about the grass-user mailing list