GRASS and GMT: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
(→‎Vector: v.out.ogr)
(→‎Importing CPT color tables: update cpt-city URL)
 
(31 intermediate revisions by 4 users not shown)
Line 8: Line 8:


=== Raster ===
=== Raster ===
* Use the <tt>grdinfo</tt> utility from GMT to check the ''Grid file format''. See the <tt>grdreformat</tt> man page for an explanation of the codes.
** Codes starting with '''n''' are "new format" (NetCDF) which can be read by [http://www.gdal.org/frmt_various.html#GMT GDAL] tools and imported with the {{cmd|r.in.gdal}} module.
** Codes starting with '''b''' are "old format" which can be imported with the {{cmd|r.in.bin}} module.


==== Old style binary grd ====
==== Old style binary grd ====


* Import using '''r.in.bin''' with the -h flag.
* Import using {{cmd|r.in.bin}} with the -h flag.
: ''Note:'' (in case you missed it) '''r.in.bin -h''' was written to be used with the old style GMT grd format.
: ''Note:'' (in case you missed it) '''r.in.bin -h''' was written to be used with the old style GMT grd format.
<!-- insert 'r.in.bin -h' command line example here -->
<!-- insert 'r.in.bin -h' command line example here -->
* ''Grid file format'' codes (see above) ending in <tt>f</tt> contain floating point data and require that ''r.in.bin'' be passed the '''-f''' flag. Format codes ending in <tt>d</tt> contain double precision floating point data and require the '''-d''' flag.


==== New style NetCDF grd ====
==== New style NetCDF grd ====


* {{cmd|r.in.gdal}} should be used for the new NetCDF GMT grd format directly.
* {{cmd|r.in.gdal}} is used to import the new NetCDF GMT grd format.
: See the [http://www.gdal.org/frmt_various.html#GMT GDAL netCDF driver] page.
: See the [http://www.gdal.org/frmt_various.html#GMT GDAL netCDF driver] page.
<!-- insert r.in.gdal command line example here -->
<!-- insert r.in.gdal command line example here -->


* If you wish to use the '''r.in.bin -h''' method with a NetCDF based grd file you will need to reformat the GMT grid into the old format:
* 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 <tt>grdreformat</tt> tool to perform the conversion:
  grdreformat in.grd out.grd=bf
  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 <tt>grdinfo</tt> utility's output may give you some help, as will the output of <tt>gdalinfo</tt>.
For an example of importing a GMT grid into GRASS, see [[Marine_Science#Bathymetric_data|this page]].


=== Vectors ===
=== Vectors ===
Line 32: Line 46:
==== Importing CPT color tables ====
==== Importing CPT color tables ====


* Import GMT style .cpt color table files with the '''r.cpt2grass''' script on the [[GRASS_AddOns#r.cpt2grass|Add-ons]] page.
* Import GMT style .cpt [[Color tables|color table]] files with the '''r.cpt2grass''' script on the [[GRASS_AddOns#r.cpt2grass|Add-ons]] page.
** [[GRASS_AddOns#r.cpt2grass|r.cpt2grass]] - Import GMT color rules
** [[GRASS_AddOns#r.cpt2grass|r.cpt2grass]] - Import GMT color rules
: For a large collection of GMT .cpt files see http://sview01.wiredworkplace.net/pub/cpt-city/
: For a large collection of GMT .cpt files see http://soliton.vm.bytemark.co.uk/pub/cpt-city/


== Exporting GRASS maps to GMT ==
== Exporting GRASS maps to GMT ==
Line 44: Line 58:


* [[GRASS_AddOns#v.out.gmt|v.out.gmt]] - Vector export for GMT (rewrite to use new [http://www.gdal.org/ogr/drv_gmt.html OGR r/w support]?)
* [[GRASS_AddOns#v.out.gmt|v.out.gmt]] - Vector export for GMT (rewrite to use new [http://www.gdal.org/ogr/drv_gmt.html OGR r/w support]?)
* [[GRASS_AddOns#v.out.gmt|r.out.gmt]] - Raster export for GMT
* [[GRASS_AddOns#r.out.gmt|r.out.gmt]] - Raster export for GMT
* [[GRASS_AddOns#v.out.gmt|r.out.gmt2]] - Raster export for GMT
* [[GRASS_AddOns#r.out.gmt2|r.out.gmt2]] - Raster export for GMT (updated version of the previous one)




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


* 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: [http://david.p.finlayson.googlepages.com/gisscripts r.out.gmt.py], [http://bambi.otago.ac.nz/hamish/grass/r.out.gmt r.out.gmt] (Hamish and Dylan), [http://169.237.35.250/~dylan/grass_user_group/r.out.gmt.sh r.out.gmt.sh] (Dylan, based Hamish's work).


=== Vector ===
=== Vector ===
Line 61: Line 73:


* OGR now supports GMT vector  (multiline) output, so it is easy to convert data to GMT format for plotting. Use the {{cmd|v.out.ogr}} module.
* OGR now supports GMT vector  (multiline) output, so it is easy to convert data to GMT format for plotting. Use the {{cmd|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'''
<!-- doesn't gmt do this automatically?
* If you are only plotting a small region of a large vector map you may want to clip the lines and polygons to be within the region with {{cmd|v.in.region}} + {{cmd|v.select}} or {{cmd|v.overlay}}.
-->
==== Examples ====
''' Export and plot areas '''
Export area vectors from a non-geographic (projected) location.
Uses the "trn_sites" map from the Spearfish sample dataset.
<source lang="bash">
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
</source>
<tt>psxy</tt> options explained:
<source lang="bash">
-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
</source>
If exporting lat/lon data you might skip the rounding outwards hack shown above and change the <tt>psxy</tt> options as follows:
<source lang="bash">
  -Jm15
  -B30m/15m:."Map Title":WeSn
  --PLOT_DEGREE_FORMAT=ddd:mmA
</source>
See GMT's <tt>psxy</tt>, <tt>psbasemap</tt>, and <tt>gmtdefaults</tt> 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 <tt>psxy</tt> to keep the PostScript file open, then the second time you run <tt>psxy</tt> 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.
<source lang="bash">
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
</source>
''' Fancy border '''
<source lang="bash">
-B30mf10m/15mf5m::WeSn
</source>
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: [http://david.p.finlayson.googlepages.com/gisscripts r.out.gmt.py], {{AddonCmd|r.out.gmt}} (Hamish and Dylan), {{AddonCmd|r.out.gmt2}} (Huidae's upgraded version of r.out.gmt), [http://169.237.35.250/~dylan/grass_user_group/r.out.gmt.sh 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 <tt>grdimage</tt> 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 <tt>grdreformat</tt> man page for a list of two-letter codes).
<source lang="bash">
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
</source>


=== Coloring ===
=== Coloring ===
Line 66: Line 219:
* Conversion of GRASS raster color data into GMT compatible CPT files
* Conversion of GRASS raster color data into GMT compatible CPT files
David Finlayson's [http://david.p.finlayson.googlepages.com/gisscripts 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.
David Finlayson's [http://david.p.finlayson.googlepages.com/gisscripts 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 {{AddonCmd|r.out.gmt}} and {{AddonCmd|r.out.gmt2}} modules automatically create a .cpt file for their grids during export.


=== Symbology ===
=== Symbology ===
* GMT has its own symbol symbol format, [http://gmt.soest.hawaii.edu/gmt/doc/gmt/html/GMT_Docs/node139.html specified here].
: It is ''very'' similar to GRASS's {{cmd|d.graph}} command format, so should be pretty easy to write something that converts them over.
: As for GRASS's [[IconSymbols|GRASS's symbol format]], perhaps export a single symbol only a blank page with {{cmd|ps.map}} or the {{cmd|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
* 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.
Ideas expressed on various mailing list, haven't seem much since. It ''should'' be a relatively simple excercise in XML parsing to convert symbology stored in a QGIS project file into something that GMT can use.


=== Development plans ===
=== Development plans ===
Line 85: Line 245:
* [http://netgis.geo.uw.edu.pl/srtm/ Dane DEM uzyskane z radarowej misji promu kosmicznego - SRTM-3]: Tutorial in Polish
* [http://netgis.geo.uw.edu.pl/srtm/ Dane DEM uzyskane z radarowej misji promu kosmicznego - SRTM-3]: Tutorial in Polish


[[Category:Documentation]]
[[Category: Documentation]]
[[Category: Tutorial]]
[[Category: Linking to other languages]]

Latest revision as of 21:38, 16 June 2016

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://soliton.vm.bytemark.co.uk/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