Tag Archives: GMT

The world in (geological) colors

The color palettes I described in my recent post can also be used with GPlates to color loaded features. Here are three examples which use the static polygon dataset which comes with GPlates. The first is using the GTS2012 chronostratigraphic time scale to color the static polygons by epoch and fill the polygons:

The world colored by geological epoch.

The world colored by geological epoch, filling the static polygons.

The second is just using the same data set and color the line features by geological era:

Coloring the polygons by era in GTS2012 colors.

Coloring the polygons by era in GTS2012 colors.

And the last is a bit Bauhaus-y by coloring everything in black and white according to the Gee & Kent geomagnetic polarity time scale (not that this would make a lot of sense, but as we can do it…):

The world in normal polarity and reverse polarity according to the Gee & Kent geomagnetic polarity timescale.

The world in normal polarity and reverse polarity according to the Gee & Kent geomagnetic polarity timescale (using the filled static polygons).

That planet looks way to boring -  hemispherical view of the Earth colored according to a normal (white) and reverse (black) geomagnetic polarity. Again using the Gee & Kent 2007 timescale.

That planet looks way to boring – hemispherical view of the Earth with the static polygons colored according to a normal (white) and reverse (black) geomagnetic polarity. Again using the Gee & Kent 2007 timescale.

In order to color loaded features by age (and timescale) just add the colorpalettes to the “Draw style” settings (Features -> Manage colouring) like this:

The geological age color palettes can be added to the Draw style (Manage colouring).

The geological age color palettes can be added to the Draw style (Manage colouring).

Once the new colour palettes are available, they can be assigned to the individual layers either through the layer window or through “Features -> “Manage colouring” .

Advertisements

Plate boundaries for GMT5’s psxy

GMT5‘s psxy offers a nice feature which lets the user designate fronts on a line segment basis – for example if you have a number of normal faults and subduction zone line segments in a multi-segment file, you can give each line a different symbol, depending on to which side the fault dips (left or right in line direction).

I have taken Peter Bird’s plate boundary file  and added these plotting instructions to the indivdual plate boundary line segment headers, for example:

> JF\NA by Peter Bird 1999
> -Sf0.45/3p+t+r

for the Juan de Fuca plate subducting under NorthAmerica. With the new header, psxy will now plot a triangle (+t) with a 3pt size (3p) to the right side (+r) of the front, fronts being separated by 0.45cm. And this is how it looks on a Robinson map centered on the dateline, generated with:


pscoast -Rg -JN180/22 -Dl -A5000 -V -G200 -K -Y4 > Bird_PlateBoundaries.ps
psxy -J -R PB2002_boundaries.gmt -W0.5p,red -O -Sf0.25/3p -Gred -Ba30:."Bird's 2003 Plate Boundaries" >> Bird_PlateBoundaries.ps

Plate boundaries plotted using GMT5’s psxy and fronts specified on line header segments


The file can be downloaded here (please make sure you cite Peter Bird’s paper as original source of the data!). The file is released under the same license as the original data.

Here’s the relevant snippet from the psxy man page for the use of -Sf in segment headers:

-Sf front. -Sfgap/size[+l|+r][+b+c+f+s+t][+ooffset]. Supply distance gap between symbols and symbol size.
If gap is negative, it is interpreted to mean the number of symbols along the front instead. Append +l
or BD+r) to plot symbols on the left or right side of the front [Default is centered]. Append +type to
specify which symbol to plot: box, circle, fault, slip, or triangle. [Default is fault]. Slip means
left-lateral or right-lateral strike-slip arrows (centered is not an option). Append +ooffset to offset
the first symbol from the beginning of the front by that amount [0]. Note: By placing -Sf options in the
segment header you can change the front types on a segment-by-segment basis.

GIS and p(l)ain text

So you are working with geospatial data. You are  collaborating with several people on the same dataset. People in your team are on different OS (Mac, Win or Linux) and want to use different geospatial tools, like QGIS, GPlates, GMT, OpenJUMP or ArcGIS or Matlab as they all have different requirements and  used to different workflows in their geoscientific research. You would like to keep track of the changes made to that specific dataset and snapshot it at different stages — ideally through SVN, git or any other revisioning tool. You don’t have any money and probably even less time.

So (unless I am mistaken) there are some options right at hand:

  1. Who cares about money: stuff all that open source software (who uses that anyway…), convert everyone to M$ Windoze and force them to use the one and only mighty ESRI ArcGIS. Yeah…NOT really.
  2. Put a lot of effort into setting up PostGIS plus versioning and then spend the rest of your life on figuring out how to connect ArcGIS in a way that you can read and write to that DB.
  3. Put the good old shp file into a a revisioning system. Hm, not too great for binary files and if you’d like to check differences on a single feature between two versions…
  4. Give up on the idea of revisioning and just make every user to save snapsots manually and store them in a central location?
Strangely, in 2012 there seems to be not a single “open”, non-binary file format which can be edited and read across the whole FOSS and proprietary GIS world – at least to my knowledge. A potential candidate, with little overhead is the GMT OGR format (see the documentation in the GMT5 cookbook – PDF file link here, which is produced by ogr2ogr, when converting shapefiles to GMT’s plain text format. So I guess this is what I will be doing:
  • Set up the common dataset as shapefile, add all attributes and geometries so far
  • Modify it in your GIS application of choice
  • Once you have done your changes, save the shapefile.
  • Convert the shapefile to GMT OGR using ogr2ogr
  • Put the GMT OGR file into the revisioning system.
  • Revert the process when needing to modify the data.
One could potentially write a Python script to do this in Arc but running ogr2ogr on the command line once shouldn’t be too hard…

ogr2ogr in the 3rd Dimension

I had to convert shapefile data from 2D lines/points to files consisting of 3 fields (x,y,z) for further processing using the GMT5. I was about to write some Python code to do that myself when I realised a new (GDAL >1.8) option in the ogr2ogr tool:

-zvalue Uses the specified field to fill the Z coordinate of geometries.

Turning an attribute from a shapefile into a z-Value for GMT files ogr2ogr and this flag is easy and helps to avoid ArcGIS’ 3D Analyst’s “Feature Attribute to Z-Value” function.

Simple usage:
ogr2ogr -F"GMT" OutFileWithZValues.gmt InFileWithZAttribute.shp -zfield"Thickness"

Done. Even works the other way around – from an OGR GMT file with Z-values in the 3rd column to a 2.5D shapefile:

ogr2ogr -F"ESRI Shapefile" OutFileWithZAttribute.shp InFileWithZValues.gmt -nlt LINESTRING25D