Hi,
How to specify the desired product (band) in the GFS NOAA GRIB file? Something like
wa.grb2=gd:HGT:1000 or something else
I wanted to do it without preliminary conversion (grdconvert). if there is such an opportunity.
File to test:
wa.grb2.zip (437,9 КБ)
use GMT 6.1.0 (windows 7) grdmath
Thank you
The +b flag seems to work
C:\v>grdinfo wa.grb2+b0
wa.grb2: Title: Grid imported via GDAL
wa.grb2: Command:
wa.grb2: Remark:
wa.grb2: Gridline node registration used [Geographic grid]
wa.grb2: Grid file format: gd = Import/export through GDAL
wa.grb2: x_min: -180 x_max: 179.5 x_inc: 0.5 (30 min) name: x n_columns: 720
wa.grb2: y_min: 0 y_max: 90 y_inc: 0.5 (30 min) name: y n_rows: 181
wa.grb2: v_min: 4931.02441406 v_max: 5920.3046875 name: z
wa.grb2: scale_factor: 1 add_offset: 0
wa.grb2: Default CPT:
+proj=longlat +R=6371229 +no_defs
C:\v>grdinfo wa.grb2+b1
wa.grb2: Title: Grid imported via GDAL
wa.grb2: Command:
wa.grb2: Remark:
wa.grb2: Gridline node registration used [Geographic grid]
wa.grb2: Grid file format: gd = Import/export through GDAL
wa.grb2: x_min: -180 x_max: 179.5 x_inc: 0.5 (30 min) name: x n_columns: 720
wa.grb2: y_min: 0 y_max: 90 y_inc: 0.5 (30 min) name: y n_rows: 181
wa.grb2: v_min: 4931.02441406 v_max: 5920.3046875 name: z
wa.grb2: scale_factor: 1 add_offset: 0
wa.grb2: Default CPT:
+proj=longlat +R=6371229 +no_defs
Compare, please, with gdalinfo min/max output for 2nd dataset
STATISTICS_MAXIMUM=450.89669799805
STATISTICS_MINIMUM=-201.15130615234
Alex
Ops, looks that I read the output too quickly. Specially because I’d tried with GMT.jl before and it worked
julia> gmtread("wa.grb2", band=2)
A GMTgrid object with 1 layers of type Float32
title: Grid imported via GDAL
Gridline node registration used
x_min: -180.0 x_max :179.5 x_inc :0.5 n_columns :720
y_min: 0.0 y_max :90.0 y_inc :0.5 n_rows :181
z_min: -201.15130615234375 z_max :450.8966979980469
Mem layout: BCB
PROJ: +proj=longlat +R=6371229 +no_defs
julia> gmtread("wa.grb2", band=1)
A GMTgrid object with 1 layers of type Float32
title: Grid imported via GDAL
Gridline node registration used
x_min: -180.0 x_max :179.5 x_inc :0.5 n_columns :720
y_min: 0.0 y_max :90.0 y_inc :0.5 n_rows :181
z_min: 4927.90478515625 z_max :5923.28466796875
Mem layout: BCB
PROJ: +proj=longlat +R=6371229 +no_defs
So we may have a bug in the +b, or it’s because the grib asks for a different procedure.
Hi,
Is there a difference between grdinfo file+b0 and grdinfo file=gd+b0 ?
It shouldn’t be. The =gd form was introduced to force reading through GDAL when an alternative way exists in GMT itself. Otherwise depending on the format the decision is taken by a guess algorithm. And in the GRIB case the only possibility is to read via GDAL.
What I suspect is happening here is that the +b, which has been designed to work with images only, is being ignored in the GRIB (a grid) case.
This explains the lack of mention of the option in the documentation. 
There is documentation saying more or less exactly what Joaquim wrote:
Modifiers to read and write grids and images via GDAL
It’s just not in the manual pages for individual modules.
Yes, its working
grdmath w2.grb2=gd+b0 w2.grb2=gd+b1 SUB = sub.nc=nf
So I was confused by the incorrect grdinfo output and i found very simple and compact syntax
grdinfo wa.grb2=gd:1
Thank You
Yes, I got surprised too that grdinfo doesn’t work but the other modules do.
Joaquim,
i started with grdcontour and does not find the answer. All this forms falls with the same error
grdcontour [ERROR]: Computed -srcwin falls outside raster size of 720x181.
grdcontour [ERROR]: ERROR reading file with gdalread.
grdcontour (gmtapi_import_grid): Could not open file [wa.grb2=gd:1]
Error returned from GMT API: GMT_GRID_READ_ERROR (18)
grdcontour “wa.grb2=gd:1” -JS75/90/61c -R-30/180/30/85 -S6 -Gd20 -W1p -A5+s26 -P >foo.ps
or gd=1 or =gd+b1 or +b1 and with the other projections and for gmt 6.6.0 version too
Alex
It is probably more surprising some things are working rather than the opposite.
In a situation like this I would rather stick to something that is known working. E.g.
gmt math w2.grb2=gd+b0 = w2_b0.nc
gmt grdcontour w2_b0.nc ...
PS I tried extracting band 2 using gdal_translate. gmt grdinfo provides wrong min/max values on the resulting wa_band2.grb2:
gdal_translate -b 2 wa.grb2 wa_band2.grb2
Input file size is 720, 181
0...10...20...30...40...50...60...70...80...90...100 - done.
gmt grdinfo wa_band2.grb2
wa_band2.grb2: Title: Grid imported via GDAL
wa_band2.grb2: Command:
wa_band2.grb2: Remark:
wa_band2.grb2: Gridline node registration used [Geographic grid]
wa_band2.grb2: Grid file format: gd = Import/export through GDAL
wa_band2.grb2: x_min: -180 x_max: 179.5 x_inc: 0.5 (30 min) name: x n_columns: 720
wa_band2.grb2: y_min: 0 y_max: 90 y_inc: 0.5 (30 min) name: y n_rows: 181
wa_band2.grb2: v_min: -169.839309692 v_max: 320.112701416 name: z
wa_band2.grb2: scale_factor: 1 add_offset: 0
wa_band2.grb2: Default CPT:
+proj=longlat +R=6371229 +no_defs
then I run gdalinfo -stats that generates wa_band2.grb2.aux.xml and suddenly gmt grdinfo wa_band2.grb2 is providing correct min/max values - only until the .aux.xml file is removed:
gdalinfo wa_band2.grb2 -nomd -stats
...
Minimum=-201.151, Maximum=450.897, Mean=104.480, StdDev=75.212
gmt grdinfo wa_band2.grb2
...
wa_band2.grb2: v_min: -201.151306152 v_max: 450.896697998 name: z
wa_band2.grb2: scale_factor: 1 add_offset: 0
...
rm wa_band2.grb2.aux.xml
gmt grdinfo wa_band2.grb2
...
wa_band2.grb2: v_min: -169.839309692 v_max: 320.112701416 name: z
...
So -R is failing here (grdcut fails with same error message). Regarding the failure of grdinfo that is the GDAL side I think. Sometimes has to be told to recompute the file stats when saving files. grdinfo also has an option to force computing the stats instead of relying on what is in the grid’s header.
I’ll try to see what is going on with the -R thing when I have some opportunity.