GRIB: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
(updated)
 
(34 intermediate revisions by 5 users not shown)
Line 1: Line 1:
* ''see also the [[Meteorology]] wiki page''
== GRIB data ==
== GRIB data ==


Meteorological data from the WMO is supplied in the GRIB format. See the [http://www.gdal.org/frmt_grib.html GDAL GRIB format page]
Meteorological data from the [http://www.wmo.int WMO] is supplied in GRIB format. See the [http://www.gdal.org/frmt_grib.html GDAL GRIB format page].
 
* Some data available here:
** http://www.globalmarinenet.com/free-grib-file-downloads/
** http://www.raymarine.com/default.aspx?site=1&section=3&page=149
** http://www.navcenter.com/newdwnld.ccml
** https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs


Newer versions of [http://www.gdal.org GDAL] can read this data (i.e. newer than the 1.5 branch; SVN trunk is best), and it can be imported with {{cmd|r.in.gdal}}.
* GRIB data viewer:
** zyGRIB (GPL3)
: http://www.zygrib.org (French; translate with [http://babelfish.yahoo.com/ Babelfish])


It often contains a small number of cells but many many bands (layers).


View band metadata with <tt>gdalinfo</tt>.
[http://www.gdal.org GDAL] can read this data, and it can be imported with {{cmd|r.in.gdal}}.
* '''''Make sure you are using GDAL version 1.10 or newer due to bugs in earlier versions of GDAL's GRIB driver.'''''


GRIB files often contains a small number of cells but many many bands (layers).


* ''see also the [[Meteorology]] wiki page''
View band metadata with <tt>gdalinfo</tt>.




Line 17: Line 28:
A mini-tutorial follows.
A mini-tutorial follows.


Data downloaded from http://gribs.ocens.net
wget "http://gribs.ocens.net/TasmanSea.wind.grb.bz2"


The data is in Lat/Lon on a Sphere. It is not quite correct to import into a Lat/Lon WGS84 location, but the resolution is so crude that you might be able to hold your nose and do so anyway.
The data is in Lat/Lon on a Sphere. It is not quite correct to import into a Lat/Lon WGS84 location, but the resolution is so crude that you might be able to hold your nose and do so anyway.


* ''There may be a cell-registration bug: (grid vs cell-center convention) cells all offset to the east and south?''
* ''There was a cell-registration bug in GDAL older than version 1.10: (grid vs cell-center convention) cells were all offset to the east and south by half a cell. If you are stuck with an old version of GDAL you can fix this by hand with the {{cmd|r.region}} module. See:''
: https://trac.osgeo.org/gdal/ticket/2550
: https://trac.osgeo.org/gdal/ticket/2637 # '''Note''': these 2 bugs appear now closed and fixed
: GRASS's [[Grid registration]] wiki page for a general discussion


==== GDAL bug workaround ====


Convert latest u,v components of wind to GeoTiff (not really needed, but easier to debug) and import into GRASS:
* Versions of GDAL '''prior''' to 1.10 had a 1/2 cell registration bug. It is fixed now in newer versions of GDAL. What follows is a work-around for that old bug:


No-data (NULL) values are flagged as "9999" in the file.
<source lang="bash">
g.region rast=mapname
eval `g.region -pg`
HALF_NSRES=`echo "$nsres" | awk '{printf("%.15g", $1 / 2.)}'`
HALF_EWRES=`echo "$ewres" | awk '{printf("%.15g", $1 / 2.)}'`
NEW_N=`echo "$n" "$HALF_NSRES" | awk '{printf("%.15g", $1 + $2)}'`
NEW_S=`echo "$s" "$HALF_NSRES" | awk '{printf("%.15g", $1 + $2)}'`
NEW_E=`echo "$e" "$HALF_EWRES" | awk '{printf("%.15g", $1 - $2)}'`
NEW_W=`echo "$w" "$HALF_EWRES" | awk '{printf("%.15g", $1 - $2)}'`
r.region mapname n=$NEW_N s=$NEW_S e=$NEW_E w=$NEW_W
</source>


  $ gdal_translate -b 1 -a_nodata 9999 TasmanSea.wind.grb TasmanSea.wind.u_wind_latest.tif
* ''Currently band metadata is not automatically copied with the import. Add it in manually with'' "<tt>r.support history=</tt>".
  $ gdal_translate -b 2 -a_nodata 9999 TasmanSea.wind.grb TasmanSea.wind.v_wind_latest.tif
 
 
Convert latest u,v components of wind to GeoTiff (not really needed, but easier to debug) and import into GRASS. Use <tt>gdalinfo</tt> to view the list of available bands and determine the needed band number(s).
 
No-data (NULL) values are flagged as "9999" in this file.
 
  $ gdal_translate -b 1 -a_nodata 9999 TasmanSea.wind.grb TasmanSea.wind.u_wind_latest.tif
  $ gdal_translate -b 2 -a_nodata 9999 TasmanSea.wind.grb TasmanSea.wind.v_wind_latest.tif
   
   
  r.in.gdal in=TasmanSea.wind.u_wind_latest.tif \
  r.in.gdal -o in=TasmanSea.wind.u_wind_latest.tif out=TasmanSea.wind.u_wind_latest
  out=TasmanSea.wind.u_wind_latest -o
  r.in.gdal -o in=TasmanSea.wind.v_wind_latest.tif out=TasmanSea.wind.v_wind_latest
  r.in.gdal in=TasmanSea.wind.v_wind_latest.tif \
  out=TasmanSea.wind.v_wind_latest -o


or import into GRASS directly:
or import into GRASS directly:
Line 45: Line 74:




Convert U,V velocity component maps into magnitide,direction maps for use with d.rast.arrow:
Convert U,V velocity component maps into magnitude, direction maps for use with d.rast.arrow:
  g.copy TasmanSea.wind.u_wind_latest,U_map
  g.copy raster=TasmanSea.wind.u_wind_latest,U_map
  g.copy TasmanSea.wind.v_wind_latest,V_map
  g.copy raster=TasmanSea.wind.v_wind_latest,V_map
   
   
  g.region rast=TasmanSea.wind.u_wind_latest -p
  g.region raster=TasmanSea.wind.u_wind_latest -p
  r.mapcalc 'magnitude = sqrt(U_map^2 + V_map^2)'
  r.mapcalc 'magnitude = sqrt(U_map^2 + V_map^2)'
r.colors magnitude col=bcyr
r.support magnitude units="m/s"
  r.mapcalc 'direction = atan(U_map, V_map)'
  r.mapcalc 'direction = atan(U_map, V_map)'


Convert m/s to knots:
r.mapcalc 'magnitude_kts = magnitude * 3600/1852.0'
r.support magnitude_kts units="knots"
r.colors magnitude_kts col=bcyr


Display:
Display:
  d.mon x0
  d.mon wx0
r.colors magnitude col=bcyr
  d.rast magnitude
  d.rast magnitude
  d.vect admin98 type=area fcol=225:225:255 color=none
  d.vect admin98 type=boundary
   
   
  d.rast.arrow map=direction type=grass magnitude_map=magnitude \
  d.rast.arrow map=direction type=grass magnitude_map=magnitude grid=none color=black
  grid=none arrow=black
   
   
  d.legend map="magnitude" at=90.3,94.7,1.8,26.5
  d.legend map="magnitude" at=90.3,94.7,1.8,26.5
Line 73: Line 106:
: [[Image:TasmanSea winds latest.jpg]]
: [[Image:TasmanSea winds latest.jpg]]


Create 2 m/s contours:
==== Create contours ====
  r.contour in=magnitude out=current_wind_2ms step=2
 
Create smooth contour lines at 5 knot intervals. First resample to a higher raster resolution with {{cmd|r.resamp.interp}} so the lines come out smooth.
 
g.region res=0:15 -p
r.resamp.interp in=magnitude_kts out=magnitude_kts.resamp.cube method=bicubic
r.contour in=magnitude_kts.resamp.cube out=latest_wind_5kt step=5
 
 
Alternatively the {{cmd|v.generalize}} module may be used to smooth the contours. The ''Snakes'' spline method is used here as it smooths by minimization of the "energy" in the line, which mimics the wind's natural tendency to strive for geostrophic balance.
 
g.region rast=magnitude_kts
  r.contour in=magnitude_kts out=latest_wind_5kt_raw step=5
v.generalize in=latest_wind_5kt_raw out=latest_wind_5kt_Snakes method=snakes alpha=0.25 beta=0.25
 
 
To remove any tiny rings and slivers you can add a length column in the contours map with {{cmd|v.db.addcolumn}}, populate the column with line length using the {{cmd|v.to.db}} module, and finally extract only lines with a length greater than some threshold with a SQL query to {{cmd|v.extract}}.
 
=== Automation ===
 
You can automate the generation of a weather map with a little GRASS scripting (e.g. as a {{wikipedia|cron_job}}). After downloading the latest forecast with <tt>wget</tt> and extracting timestamp info with <tt>gdalinfo</tt>, with a <tt>GRASS_BATCH_JOB</tt> you can import the data into a lat/lon location with {{cmd|r.in.gdal}} and then in a following batch job from the master script use {{cmd|r.proj}} to reproject the data into another location/projection (if you choose) and run the appropriate d.* rendering commands. See the [[GRASS_and_Shell#GRASS_Batch_jobs]] wiki page and the main GRASS man page for more details.
 
=== TODO ===
 
* d.barb module: ''work in progress by Hamish'' (partially done; low on the priority list currently; still plan to finish it)
 
 
[[Category:Documentation]]
[[Category:Tutorial]]
[[Category:Grib]]

Latest revision as of 10:47, 31 October 2016

GRIB data

Meteorological data from the WMO is supplied in GRIB format. See the GDAL GRIB format page.

  • GRIB data viewer:
    • zyGRIB (GPL3)
http://www.zygrib.org (French; translate with Babelfish)


GDAL can read this data, and it can be imported with r.in.gdal.

  • Make sure you are using GDAL version 1.10 or newer due to bugs in earlier versions of GDAL's GRIB driver.

GRIB files often contains a small number of cells but many many bands (layers).

View band metadata with gdalinfo.


Processing

A mini-tutorial follows.


The data is in Lat/Lon on a Sphere. It is not quite correct to import into a Lat/Lon WGS84 location, but the resolution is so crude that you might be able to hold your nose and do so anyway.

  • There was a cell-registration bug in GDAL older than version 1.10: (grid vs cell-center convention) cells were all offset to the east and south by half a cell. If you are stuck with an old version of GDAL you can fix this by hand with the r.region module. See:
https://trac.osgeo.org/gdal/ticket/2550
https://trac.osgeo.org/gdal/ticket/2637 # Note: these 2 bugs appear now closed and fixed
GRASS's Grid registration wiki page for a general discussion

GDAL bug workaround

  • Versions of GDAL prior to 1.10 had a 1/2 cell registration bug. It is fixed now in newer versions of GDAL. What follows is a work-around for that old bug:
 g.region rast=mapname
 eval `g.region -pg`
 HALF_NSRES=`echo "$nsres" | awk '{printf("%.15g", $1 / 2.)}'`
 HALF_EWRES=`echo "$ewres" | awk '{printf("%.15g", $1 / 2.)}'`
 NEW_N=`echo "$n" "$HALF_NSRES" | awk '{printf("%.15g", $1 + $2)}'`
 NEW_S=`echo "$s" "$HALF_NSRES" | awk '{printf("%.15g", $1 + $2)}'`
 NEW_E=`echo "$e" "$HALF_EWRES" | awk '{printf("%.15g", $1 - $2)}'`
 NEW_W=`echo "$w" "$HALF_EWRES" | awk '{printf("%.15g", $1 - $2)}'`
 
 r.region mapname n=$NEW_N s=$NEW_S e=$NEW_E w=$NEW_W
  • Currently band metadata is not automatically copied with the import. Add it in manually with "r.support history=".


Convert latest u,v components of wind to GeoTiff (not really needed, but easier to debug) and import into GRASS. Use gdalinfo to view the list of available bands and determine the needed band number(s).

No-data (NULL) values are flagged as "9999" in this file.

$ gdal_translate -b 1 -a_nodata 9999 TasmanSea.wind.grb TasmanSea.wind.u_wind_latest.tif
$ gdal_translate -b 2 -a_nodata 9999 TasmanSea.wind.grb TasmanSea.wind.v_wind_latest.tif

r.in.gdal -o in=TasmanSea.wind.u_wind_latest.tif out=TasmanSea.wind.u_wind_latest
r.in.gdal -o in=TasmanSea.wind.v_wind_latest.tif out=TasmanSea.wind.v_wind_latest

or import into GRASS directly:

r.in.gdal in=TasmanSea.wind.grb band=1 out=TasmanSea.wind.u_wind_latest
r.in.gdal in=TasmanSea.wind.grb band=2 out=TasmanSea.wind.v_wind_latest
r.null TasmanSea.wind.u_wind_latest setnull=9999
r.null TasmanSea.wind.v_wind_latest setnull=9999


Convert U,V velocity component maps into magnitude, direction maps for use with d.rast.arrow:

g.copy raster=TasmanSea.wind.u_wind_latest,U_map
g.copy raster=TasmanSea.wind.v_wind_latest,V_map

g.region raster=TasmanSea.wind.u_wind_latest -p
r.mapcalc 'magnitude = sqrt(U_map^2 + V_map^2)'
r.colors magnitude col=bcyr
r.support magnitude units="m/s"
r.mapcalc 'direction = atan(U_map, V_map)'

Convert m/s to knots:

r.mapcalc 'magnitude_kts = magnitude * 3600/1852.0'
r.support magnitude_kts units="knots"
r.colors magnitude_kts col=bcyr

Display:

d.mon wx0
d.rast magnitude
d.vect admin98 type=boundary

d.rast.arrow map=direction type=grass magnitude_map=magnitude grid=none color=black

d.legend map="magnitude" at=90.3,94.7,1.8,26.5
echo "wind (m/s)" | d.text color=black at=13.8,80.6 size=3 align=lc

Export to a PNG image using the Cairo driver:

d.out.file -c TasmanSea_winds_latest


Resulting image:

Create contours

Create smooth contour lines at 5 knot intervals. First resample to a higher raster resolution with r.resamp.interp so the lines come out smooth.

g.region res=0:15 -p
r.resamp.interp in=magnitude_kts out=magnitude_kts.resamp.cube method=bicubic
r.contour in=magnitude_kts.resamp.cube out=latest_wind_5kt step=5


Alternatively the v.generalize module may be used to smooth the contours. The Snakes spline method is used here as it smooths by minimization of the "energy" in the line, which mimics the wind's natural tendency to strive for geostrophic balance.

g.region rast=magnitude_kts
r.contour in=magnitude_kts out=latest_wind_5kt_raw step=5
v.generalize in=latest_wind_5kt_raw out=latest_wind_5kt_Snakes method=snakes alpha=0.25 beta=0.25


To remove any tiny rings and slivers you can add a length column in the contours map with v.db.addcolumn, populate the column with line length using the v.to.db module, and finally extract only lines with a length greater than some threshold with a SQL query to v.extract.

Automation

You can automate the generation of a weather map with a little GRASS scripting (e.g. as a cron_job). After downloading the latest forecast with wget and extracting timestamp info with gdalinfo, with a GRASS_BATCH_JOB you can import the data into a lat/lon location with r.in.gdal and then in a following batch job from the master script use r.proj to reproject the data into another location/projection (if you choose) and run the appropriate d.* rendering commands. See the GRASS_and_Shell#GRASS_Batch_jobs wiki page and the main GRASS man page for more details.

TODO

  • d.barb module: work in progress by Hamish (partially done; low on the priority list currently; still plan to finish it)