[GRASSLIST:864] GMT / GRASS use / shp2gmt

GRASS GRASS user
,
1

Grudgingly I'm using GMT after working through Dylan's example and am
getting it to be somewhat functional.

I will be needing to script the building of a large number of maps
which I may store either in PostGIS or GRASS.

I've noticed in the use of shp2gmt that there is no way of adding in a
z attribute for use with the thematic colour assignment feature in GMT.
If necessary I can extract bits based on different attributes and add
them in separately, but it would be nice to be able to simply define a
colour them once and be done with it. Any suggestions on how I can turn
an attribute column in GRASS into a z-element that will be converted
with shp2gmt would be appreciated.

T
--
Trevor Wiens
twiens@interbaun.com

The significant problems that we face cannot be solved at the same
level of thinking we were at when we created them.
(Albert Einstein)

2

Hi Trevor:

(sorry for any cross-posts, hoping to get some extra traction on this topic)

On Tuesday 25 April 2006 01:13 pm, Trevor Wiens wrote:

Grudgingly I'm using GMT after working through Dylan's example and am
getting it to be somewhat functional.

:slight_smile: sorry for any slightly outdated material on that site
(http://169.237.35.250/~dylan/grass_user_group/map1.html). I have a lot of
changes that I would like to make, along with some more examples...
Unfortunately I don't have a lot of spare time these days.

Another, simpler example using something other than UTM projected data can be
found here:
http://casoilresource.lawr.ucdavis.edu/drupal/node/102

Another (unfinished) collection of samples can be found here:
http://casoilresource.lawr.ucdavis.edu/drupal/node/51

I will be needing to script the building of a large number of maps
which I may store either in PostGIS or GRASS.

you might find the scripts "g.fieldsheet.sh" and "make_fieldsheets.sh" found
on the third URL above useful templates. I used these along with some
additional bash scripts to automate the process of making about 2 hundred
maps for soil survey operations. examples here:
http://169.237.35.250/~dylan/NPS/#fieldsheets

ask if you would like some of the details.

I've noticed in the use of shp2gmt that there is no way of adding in a
z attribute for use with the thematic colour assignment feature in GMT.
If necessary I can extract bits based on different attributes and add
them in separately, but it would be nice to be able to simply define a
colour them once and be done with it. Any suggestions on how I can turn
an attribute column in GRASS into a z-element that will be converted
with shp2gmt would be appreciated.

T

Yes. This is a major drawback for multi-class maps. Unfortunately I did not
have the time, nor expertise to customize shp2gmt for inclusion of symbology
parameters. There has been some discussion on the GMT and GDAL lists about
the possible creation of a OGR GMT driver ... but I haven't heard any news on
this topic for a while. I have been toying with the idea of a new GRASS
command called 'v.out.gmt' which would export simple geometry along with
selected symbology to the GMT format.

A starting point would be to use `v.out.ascii format=standard` along with
v.out.ascii.db to add an attribute on the starting line of each feature.

the resulting file for lines or closed polygons would look like this:

-W2/255/1/1 -G255/0/0

660835 4040837
660827 4040845
660819 4040846

-W2/255/1/255

660835 4040837
660827 4040845
660819 4040846

the '-W' and '-G' control pen and fill colors.

for points of varying symbol size, the third field is the symbol size:

-G255/0/0 -W2/0/255/0

660835 4040837 1

-G0/255/255 -W2/0/255/0

660827 4040845 9

An alternate method is to use an external CPT file for colors:

psxy [...] -L -Cthe_cpt_file.cpt > output.ps

each polygon stanza should have the following:

-Z<value>

661096.53862171 4041121.73826666
661093.98591761 4041116.63285845
661091.4332135 4041111.52745023
661083.77510118 4041108.97474613
660835.31190146 4040837.53720944

where <value> is a number referenced to a color in the CPT file

...So where to start? I think that it shouldn't be too hard to implement
simple vector->GMT conversion, based on the code found in v.out.ascii , and
as Hamish suggested a while ago ps.map (for attribute access).

A skeleton GRASS module could look something like:

src/vector/v.out.gmt/
---main.c----
1. get user args and init pointers
2. detect feature type
3. loop through features:
  1. extract & write symbology
  2. write vertex records
  3. optionally close polygon
  4. add to buffer, next feature

4. flush buffer to disk
5. clean-up and exit
------------

Apologies for the long post. Ideally it would be great to have GMT vector
support in GDAL/OGR. Raster support is nearly ready through r.out.bin |
xyz2grd -or- r.out.gdal format=NetCDF ...

Cheers,

--
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341

3

Trevor,

I'm a rank beginner at GMT but in my one exercise to date with it I
learned the following. Mark Fenbers has modified the shp2gmt.c file to
include the attribute info from the dbf file into the multi-segment
separator of the psxy input file. His code is given in his message here:
http://www.mechanik.tu-darmstadt.de/GMT-Help/Archiv/11021.html.

The output from his shp2gmt version looks like this for a tiny extract
from a state map of the US:

      Polygon #0 nVertices=17 nParts=1 cat=2 AREA=0 PERIMETER=0.224 STATESP020=3 STATE="Alaska" STATE_FIPS="02"

-2.86129e+06 6.01613e+06
-2.86152e+06 6.01607e+06
-2.86163e+06 6.01615e+06
-2.86178e+06 6.01631e+06
-2.86248e+06 6.01724e+06
-2.86327e+06 6.01807e+06
etc.

What I can then do is write an awk script to substitute in the correct
-G and -W plotting parameters depending on the attributes. For instance
running this bit of awk over the above input file

{if ($0 ~ /Alaska/) printf "> -G175 -Wblack\n";
else print;}

should transform it into an input file where Alaska is plotted with a grey
fill. With more complicated lookups in the awk script, one can do pretty
sophisticated thematic mapping using this approach.

Best wishes,

Allan Hollander

---

On Tue, 25 Apr 2006, Trevor Wiens wrote:

Grudgingly I'm using GMT after working through Dylan's example and am
getting it to be somewhat functional.

I will be needing to script the building of a large number of maps
which I may store either in PostGIS or GRASS.

I've noticed in the use of shp2gmt that there is no way of adding in a
z attribute for use with the thematic colour assignment feature in GMT.
If necessary I can extract bits based on different attributes and add
them in separately, but it would be nice to be able to simply define a
colour them once and be done with it. Any suggestions on how I can turn
an attribute column in GRASS into a z-element that will be converted
with shp2gmt would be appreciated.

T

--

4

Allan,

Great! This is certainly a step in the right direction.

Just for the record, to compile this new version of shp2gmt.c I added the
following to the makefile in the shputils dir:

shp2gmt_new: shp2gmt_new.c shpopen.o dbfopen.o
  $(CC) $(CFLAGS) shp2gmt_new.c shpopen.o dbfopen.o $(LINKOPT) -o shp2gmt_new

where shp2gmt_new.c is the code you referenced.

works like a charm!

Cheers,

Dylan

On Tuesday 25 April 2006 04:14 pm, Allan Hollander wrote:

Trevor,

I'm a rank beginner at GMT but in my one exercise to date with it I
learned the following. Mark Fenbers has modified the shp2gmt.c file to
include the attribute info from the dbf file into the multi-segment
separator of the psxy input file. His code is given in his message here:
http://www.mechanik.tu-darmstadt.de/GMT-Help/Archiv/11021.html.

The output from his shp2gmt version looks like this for a tiny extract

from a state map of the US:
> Polygon #0 nVertices=17 nParts=1 cat=2 AREA=0
> PERIMETER=0.224 STATESP020=3 STATE="Alaska" STATE_FIPS="02"

-2.86129e+06 6.01613e+06
-2.86152e+06 6.01607e+06
-2.86163e+06 6.01615e+06
-2.86178e+06 6.01631e+06
-2.86248e+06 6.01724e+06
-2.86327e+06 6.01807e+06
etc.

What I can then do is write an awk script to substitute in the correct
-G and -W plotting parameters depending on the attributes. For instance
running this bit of awk over the above input file

{if ($0 ~ /Alaska/) printf "> -G175 -Wblack\n";
else print;}

should transform it into an input file where Alaska is plotted with a grey
fill. With more complicated lookups in the awk script, one can do pretty
sophisticated thematic mapping using this approach.

Best wishes,

Allan Hollander

---

On Tue, 25 Apr 2006, Trevor Wiens wrote:
> Grudgingly I'm using GMT after working through Dylan's example and am
> getting it to be somewhat functional.
>
> I will be needing to script the building of a large number of maps
> which I may store either in PostGIS or GRASS.
>
> I've noticed in the use of shp2gmt that there is no way of adding in a
> z attribute for use with the thematic colour assignment feature in GMT.
> If necessary I can extract bits based on different attributes and add
> them in separately, but it would be nice to be able to simply define a
> colour them once and be done with it. Any suggestions on how I can turn
> an attribute column in GRASS into a z-element that will be converted
> with shp2gmt would be appreciated.
>
> T

--
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341

5

On Tue, 25 Apr 2006 16:14:35 -0700 (PDT)
Allan Hollander <adh@ice.ucdavis.edu> wrote:

Trevor,

I'm a rank beginner at GMT but in my one exercise to date with it I
learned the following. Mark Fenbers has modified the shp2gmt.c file to
include the attribute info from the dbf file into the multi-segment
separator of the psxy input file. His code is given in his message here:
http://www.mechanik.tu-darmstadt.de/GMT-Help/Archiv/11021.html.

Thanks for the helpful information. I do have a question however.

The output from his shp2gmt version looks like this for a tiny extract
from a state map of the US:

> Polygon #0 nVertices=17 nParts=1 cat=2 AREA=0 PERIMETER=0.224 STATESP020=3 STATE="Alaska" STATE_FIPS="02"
-2.86129e+06 6.01613e+06
-2.86152e+06 6.01607e+06
-2.86163e+06 6.01615e+06
-2.86178e+06 6.01631e+06
-2.86248e+06 6.01724e+06
-2.86327e+06 6.01807e+06
etc.

OK. Got that code and that is working as reported.

What I can then do is write an awk script to substitute in the correct
-G and -W plotting parameters depending on the attributes. For instance
running this bit of awk over the above input file

{if ($0 ~ /Alaska/) printf "> -G175 -Wblack\n";
else print;}

You lost me here. In what context am I going to process the file with
awk, in the psxy command line or before hand?

Please provide a slightly more fleshed out example for me to follow.

Thanks in advance.

T
--
Trevor Wiens
twiens@interbaun.com

The significant problems that we face cannot be solved at the same
level of thinking we were at when we created them.
(Albert Einstein)

6

Grudgingly I'm using GMT after working through Dylan's example and am
getting it to be somewhat functional.

I will be needing to script the building of a large number of maps
which I may store either in PostGIS or GRASS.

Any time someone wants to help with out with r.out.gmt, just chime in.
  http://grass.gdf-hannover.de/twiki/bin/view/GRASS/GrassAddOns#Raster_add_ons

I don't know about Dylan but I'm not really working on it currently.
  (but will accept updates of course!)
Dylan's other scripts are probably more up to date.

Hamish

7

Hamish,

I read r.out.gmt and it mentioned the need for a script for xyz import; take
a look at the attached script, would this fit the bill? I haven't posted it
on the WIKI yet because I was still testing it out, but if it's useful for
r.out.gmt, go nuts.

~ Eric.

-----Original Message-----
From: owner-GRASSLIST@baylor.edu
To: Trevor Wiens
Cc: dylan@iici.no-ip.org; grasslist@baylor.edu
Sent: 4/26/2006 7:59 AM
Subject: [GRASSLIST:895] Re: GMT / GRASS use / shp2gmt

Grudgingly I'm using GMT after working through Dylan's example and am
getting it to be somewhat functional.

I will be needing to script the building of a large number of maps
which I may store either in PostGIS or GRASS.

Any time someone wants to help with out with r.out.gmt, just chime in.

http://grass.gdf-hannover.de/twiki/bin/view/GRASS/GrassAddOns#Raster_add
_ons

I don't know about Dylan but I'm not really working on it currently.
  (but will accept updates of course!)
Dylan's other scripts are probably more up to date.

Hamish

(attachments)

r.in.xyz (2.73 KB)

8

Trevor,

I'm sorry that was a bit terse. Basically, both the use of shp2gmt and the
awk script are preprocessing steps before you run psxy. The processing
flow is like this:

1) Run the shp2gmt script over the shapefile to produce a
preliminary input file for psxy. In a bash shell, this would be something
like:

% shp2gmt usstates > usstates0.xy

2) Run the awk script over the preliminary input file to produce the final
version of the input file for psxy which is tagged for polygon-by-polygon
symbology. Assuming the awk script below is in a file named alaska.awk,
the command line here would be:

% awk -f alaska.awk usstates0.xy > usstates.xy

3) You can then plot this input file using psxy, as follows:

% psxy usstates.xy [other GMT parameters go here...] >> outputplot.ps

By the way, there's nothing magical about using awk -- most any scripting
language can be used for this sort of processing (a lot of people would
use perl, for instance).

Hope this helps,

Allan

On Tue, 25 Apr 2006, Trevor Wiens wrote:

On Tue, 25 Apr 2006 16:14:35 -0700 (PDT)
Allan Hollander <adh@ice.ucdavis.edu> wrote:

> Trevor,
>
> I'm a rank beginner at GMT but in my one exercise to date with it I
> learned the following. Mark Fenbers has modified the shp2gmt.c file to
> include the attribute info from the dbf file into the multi-segment
> separator of the psxy input file. His code is given in his message here:
> http://www.mechanik.tu-darmstadt.de/GMT-Help/Archiv/11021.html.

Thanks for the helpful information. I do have a question however.

>
> The output from his shp2gmt version looks like this for a tiny extract
> from a state map of the US:
>
> > Polygon #0 nVertices=17 nParts=1 cat=2 AREA=0 PERIMETER=0.224 STATESP020=3 STATE="Alaska" STATE_FIPS="02"
> -2.86129e+06 6.01613e+06
> -2.86152e+06 6.01607e+06
> -2.86163e+06 6.01615e+06
> -2.86178e+06 6.01631e+06
> -2.86248e+06 6.01724e+06
> -2.86327e+06 6.01807e+06
> etc.
>

OK. Got that code and that is working as reported.

> What I can then do is write an awk script to substitute in the correct
> -G and -W plotting parameters depending on the attributes. For instance
> running this bit of awk over the above input file
>
> {if ($0 ~ /Alaska/) printf "> -G175 -Wblack\n";
> else print;}

You lost me here. In what context am I going to process the file with
awk, in the psxy command line or before hand?

Please provide a slightly more fleshed out example for me to follow.

Thanks in advance.

T

--

9

On Wed, 26 Apr 2006 11:30:59 -0700 (PDT)
Allan Hollander <adh@ice.ucdavis.edu> wrote:

Trevor,

I'm sorry that was a bit terse. Basically, both the use of shp2gmt and the
awk script are preprocessing steps before you run psxy. The processing
flow is like this:

1) Run the shp2gmt script over the shapefile to produce a
preliminary input file for psxy. In a bash shell, this would be something
like:

% shp2gmt usstates > usstates0.xy

2) Run the awk script over the preliminary input file to produce the final
version of the input file for psxy which is tagged for polygon-by-polygon
symbology. Assuming the awk script below is in a file named alaska.awk,
the command line here would be:

% awk -f alaska.awk usstates0.xy > usstates.xy

3) You can then plot this input file using psxy, as follows:

% psxy usstates.xy [other GMT parameters go here...] >> outputplot.ps

By the way, there's nothing magical about using awk -- most any scripting
language can be used for this sort of processing (a lot of people would
use perl, for instance).

Hope this helps,

Hmm...

OK, now I understand. I think however that probably for my needs having
the ability to select objects by attribute values using SQL is
important so I will probably load my base information into PostGIS and
export my themes as different layers. I will keep this option in mind
however because this way would be faster than going PostGIS to shape
and then shape to xy.

Thanks again.

T
--
Trevor Wiens
twiens@interbaun.com

The significant problems that we face cannot be solved at the same
level of thinking we were at when we created them.
(Albert Einstein)