GRASS and GMT: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
Line 62: Line 62:




To aid custom scripting, {{cmd|g.region}} in GRASS 6.5+ supports the ''-t'' flag to output the current region in GMT's <tt>W/E/S/N</tt> region format.
To aid custom scripting, {{cmd|g.region}} in GRASS 6.5+ supports the ''-t'' flag to output the current region in GMT's <tt>W/E/S/N</tt> region format. Both GMT and GRASS support the D:M:Sh notation for latitude and longitude.




Revision as of 05:54, 24 January 2010

GMT (The Generic Mapping Tools)

GMT (Generic Mapping Tools) is a Free software package for creating publication quality cartography.

GMT homepage: http://gmt.soest.hawaii.edu

Importing from GMT

Raster

  • Use the grdinfo utility from GMT to check the Grid file format. See the grdreformat man page for an explanation of the codes.
    • Codes starting with n are "new format" (NetCDF) which can be read by GDAL tools and imported with the r.in.gdal module.
    • Codes starting with b are "old format" which can be imported with the r.in.bin module.

Old style binary grd

Note: (in case you missed it) r.in.bin -h was written to be used with the old style GMT grd format.
  • Grid file format codes (see above) ending in f contain floating point data and require that r.in.bin be passed the -f flag. Format codes ending in d contain double precision floating point data and require the -d flag.

New style NetCDF grd

  • r.in.gdal is used to import the new NetCDF GMT grd format.
See the GDAL netCDF driver page.
  • If for some reason r.in.gdal doesn't work (e.g. your copy of GDAL wasn't build with GMT/NetCDF support) you can convert the new-style grid file into the old-style grid format and then proceed to import with the r.in.bin -h method as given above. Use the GMT grdreformat tool to perform the conversion:
grdreformat in.grd out.grd=bf
Sometimes the new-style NetCDF grid contains multiple layers of data and you'll have to explicitly specify the layer in the input grid that you want. For example, if the layer is called "z":
grdreformat in.grd?z out.grd=bf
  • Examination of the grdinfo utility's output may give you some help, as will the output of gdalinfo.


For an example of importing a GMT grid into GRASS, see this page.

Vectors

Coloring

Importing CPT color tables

For a large collection of GMT .cpt files see http://sview01.wiredworkplace.net/pub/cpt-city/

Exporting GRASS maps to GMT


(Supplied by the GRASS Users Group of Davis, California)

Currently there are several *.out.GMT permutations on the Add-ons page, in several different languages (bash, python, etc.), and each of which with relative pros/cons. An effort to unify these approaches would save much of the current difficulties in moving complex raster+vector data into a GMT-friendly format. A simple road map toward this goal is outlined below.


To aid custom scripting, g.region in GRASS 6.5+ supports the -t flag to output the current region in GMT's W/E/S/N region format. Both GMT and GRASS support the D:M:Sh notation for latitude and longitude.


Vector

  • Conversion of GRASS vector data to GMT compatible ascii files

There is currently an effort (with some funding!), see some of the chatter on the GRASS and GMT mailing lists: GRASS-list GMT-help

  • OGR now supports GMT vector (multiline) output, so it is easy to convert data to GMT format for plotting. Use the v.out.ogr module.
  • If outputting areas be sure to explicitly request them, otherwise data attributes may not be transfered. The -c flag limits the export to areas with a centroid (i.e. so holes don't get filled)

For example: 

v.out.ogr -c input=RBD_F1v3 dsn=RBD_F1v3.gmt format=GMT type=area


Examples

Export and plot areas

Export area vectors from a non-geographic (projected) location. Uses the "trn_sites" map from the Spearfish sample dataset. 

MAP=trn_sites

v.out.ogr -c in=$MAP dsn=$MAP type=area format=GMT

# get map bounds
eval `v.info -g $MAP`
# ugly hack to round outward to the next whole kilometer
w=`expr \( \( $(echo $west | cut -f1 -d.) - 1000 \) / 1000 \) \* 1000`
e=`expr \( \( $(echo $east | cut -f1 -d.) + 1000 \) / 1000 \) \* 1000`
s=`expr \( \( $(echo $south | cut -f1 -d.) - 1000 \) / 1000 \) \* 1000`
n=`expr \( \( $(echo $north | cut -f1 -d.) + 1000 \) / 1000 \) \* 1000`

# GMT plotting command
psxy $MAP.gmt \
   -R$w/$e/$s/$n \
   -JXh \
   -G220/220/220 -W -M \
   -B2000/2000:."Map of $MAP":WeSn \
   --D_FORMAT='%.15g' > $MAP.ps

psxy options explained: 

 -R$w/$e/$s/$n  # bounding box
 -JXh           # map projection is non-geographic, set map-box width to fill a landscape page
 -G220/220/220  # area fill color to a light grey (R/G/B, 0-255) 
 -W             # draw area boundary
 -M             # break polygons
 -B2000:."Map of $MAP":WeSn   # tick interval : title : tick numbers on sides with capital letters
 --D_FORMAT='%.15g'  # grid reference tick number formatting

If exporting lat/lon data you might skip the rounding outwards hack shown above and change the psxy options as follows: 

  -Jm15
  -B30m/15m:."Map Title":WeSn
  --PLOT_DEGREE_FORMAT=ddd:mmA

See GMT's psxy, psbasemap, and gmtdefaults man pages for further details.


Overlay lines

Export a vector containing lines: 

v.out.ogr in=roads dsn=roads format=GMT type=line

GMT commands are same as above, but add the -K flag the first time you run psxy to keep the PostScript file open, then the second time you run psxy you can simplify the options a bit and give it the -O flag to tell it that this layer will be an overlay, not a new map. Then append (>> instead of >) to the end of the existing PostScript file. 

psxy $MAP.gmt \
   -R$w/$e/$s/$n \
   -JXh \
   -G220/240/220 -W -M \
   -B2000/2000:."Map of $MAP":WeSn \
   --D_FORMAT='%.15g' \
   -K > $MAP.ps

psxy roads.gmt -R -J -M -O >> ${MAP}.ps


Fancy border 

 -B30mf10m/15mf5m::WeSn

The above sets the border as follows: 

-B 30m f10m / 15m f5m : ."Title" : WeSn

-B    draw border
30m   x: tick numbers every 30 minutes (1/2 degree)
f10m  x: border flip every 10 minutes
15m   y: tick numbers every 15 minutes
f5m   y: border flip every 5 minutes
."t"  Set map title to string inside "quotes"
WeSn  tick numbers on W and S sides, border on all sides

Raster

  • Conversion of GRASS raster data to GMT compatible binary grids

A combination of r.out.bin | xyz2grd can accomplish this. Several attempts at generalizing this procedure have been proposed: r.out.gmt.py, r.out.gmt (Hamish and Dylan), r.out.gmt2 (Huidae's upgraded version of r.out.gmt), r.out.gmt.sh (Dylan, based Hamish's work).

  • Imagery data e.g. aerial photos

For multiband imagery data, the following commands will export the different colour bands: 

r.mapcalc "image.red=r#image; image.green=g#image; image.blue=b#image"
r.out.bin -h input=image.red output=image.red.grd
r.out.bin -h input=image.green output=image.green.grd
r.out.bin -h input=image.blue output=image.blue.grd

They can then be plotted with: 

grdimage image.red.grd image.green.grd image.blue.grd -J -R -B ...etc.

Examples

Export a raster and overlay lines

This example uses the Spearfish dataset (UTM projection). For the grdimage command the =bs after the .grd filename instructs GMT to expect a classic-GMT format integer grid file (for a CELL map); change it to =bf for a FCELL or r.out.gmt-exported DCELL map. (see the grdreformat man page for a list of two-letter codes). 


MAP=elevation.dem
g.region rast=$MAP

r.out.gmt $MAP -p

# set up plot (from instructions in last command)
psbasemap -JX7i/5.153239i -R590010/609000/4914020/4928000 \
   -Bpf300a3000:"Easting (m)":/f300a3000:"Northing (m)"::."elevation.dem":WeSn \
   -Xc -Yc -U"scale 1:106802" \
   --D_FORMAT='%.15g' \
   -K -V > $MAP.ps

# overplot exported raster grid
grdimage $MAP.grd=bs -Celevation.dem.cpt \
   -R -J -K -O -V >> $MAP.ps


# export roads vector
v.out.ogr in=roads dsn=roads format=GMT type=line

# overplot roads
psxy roads.gmt -R -J -M -O >> $MAP.ps

Coloring

  • Conversion of GRASS raster color data into GMT compatible CPT files

David Finlayson's r.out.gmt.py does a nice job of this. Once we decide on an optimal language to implement the routines in this may need translation.

  • The r.out.gmt and r.out.gmt2 modules automatically create a .cpt file for their grids during export.

Symbology

It is very similar to GRASS's d.graph command format, so should be pretty easy to write something that converts them over.
As for GRASS's GRASS's symbol format, perhaps export a single symbol only a blank page with ps.map or the d.mon PS graphics driver to get a stroked version?
  • To use GMT's .def symbols in GRASS probably it is easier to render them with d.graph or hack into a simple postscript file and plot the .eps with ps.map.
  • Automatic conversion of symbology data stored in a gis.m or QGIS saved state to GMT options
Ideas expressed on various mailing list, haven't seem much since. It should be a relatively simple exercise in XML parsing to convert symbology stored in a QGIS project file into something that GMT can use.

Development plans

  • General approach:

Since GMT relies on a sequence of specialized programs to "build-up" a postscript file, some thought must be put into how the conversion should take place. As usual, form should follow function- maximum flexibility, robustness, and accuracy being primary objectives. However, a simple means of creating high quality 2D maps would be a tremendous (I think) addition to the GRASS toolset. Especially since this is something frequently cited by critics. --DylanBeaudette 02:47, 10 December 2006 (CET)

1. should we continue down the well troden path of single-use, highly efficient programs for the various conversion steps: i.e v.out.GMT, r.out.GMT, etc.?

2. should there be a unified approach to the process: something akin to ps.map - GMT.map ?

Tutorials

Retrieved from "https://grasswiki.osgeo.org/w/index.php?title=GRASS_and_GMT&oldid=10227"