# Histograms of .grd file

**URL:** <https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941>\
**Category:** Q&A\
**Created:** [October 12, 2020, 3:11am UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941 "2020-10-12T03:11:41Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![KRL](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/krl/32/583_2.png) [@KRL](https://forum.generic-mapping-tools.org/u/KRL)\
**Post date:** [October 12, 2020, 3:11am UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/1 "2020-10-12T03:11:41Z")

</div>

Dear GMT Comunity,

I was wondering if you could help me out with this I have a .grd data set which I need to extract some statistics and make a Classification (supervised) from the .grd file. Is this possible to do with GMT?

Hope not to take to much of your time,

Kind regards,

K

---

<div class="post-metadata">

**Author:** ![Esteban82](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/esteban82/32/54_2.png) [@Esteban82](https://forum.generic-mapping-tools.org/u/Esteban82)\
**Post date:** [October 12, 2020, 3:20am UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/2 "2020-10-12T03:20:28Z")

</div>

If by “supervised clasification” you mean human-guided, then I would think that It would better to use a GIS. What exactly do you want to do?

---

<div class="post-metadata">

**Author:** ![KRL](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/krl/32/583_2.png) [@KRL](https://forum.generic-mapping-tools.org/u/KRL)\
**Post date:** [October 12, 2020, 3:36am UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/3 "2020-10-12T03:36:10Z")

</div>

I would like to make 5 groups out of a continuum dataset of Velocities then I would established ranges of velocities let´s say the continuum is from 300 m/s to 900 m/s then I would like to make 1 class from 300 to 400 m/s , 400 to 500 m/s, 500 to 600 m/s, 600 to 700 m/s and 700 to 900 m/s. based on the mean value and the dispersion fo the dataset which is interpolated with near neighbour interpolation.

of for example let’s say I would like to subtract the average velocity value and just see the max and min values i think this could be achived with grdmath not so sure really because I am still learning the function

Hope not to take too much of your time.

K

---

<div class="post-metadata">

**Author:** ![KRL](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/krl/32/583_2.png) [@KRL](https://forum.generic-mapping-tools.org/u/KRL)\
**Post date:** [October 12, 2020, 3:46am UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/4 "2020-10-12T03:46:28Z")

</div>

Also it would be great if I could see the histograms of the studied dataset. this would defenely would be great

---

<div class="post-metadata">

**Author:** ![Esteban82](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/esteban82/32/54_2.png) [@Esteban82](https://forum.generic-mapping-tools.org/u/Esteban82)\
**Post date:** [October 12, 2020, 4:01am UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/5 "2020-10-12T04:01:51Z")

</div>

The histrograms are easy. You could create a table with grd2xyz and them used (ps)histrogram to make the plot.  
For the first part, if I undesertand correctly, I think it can be done with GMT. With gmt select -Z you could filter the grd based on the velocities values, then calculated the mean and dispersion of each sub-set with grdinfo.

---

<div class="post-metadata">

**Author:** ![KRL](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/krl/32/583_2.png) [@KRL](https://forum.generic-mapping-tools.org/u/KRL)\
**Post date:** [October 12, 2020, 1:27pm UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/6 "2020-10-12T13:27:49Z")

</div>

Thank you for your advice I will give it a try ! and be back!  
Thanks !

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [October 12, 2020, 8:26pm UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/7 "2020-10-12T20:26:38Z")

</div>

An example on how the Julia wrapper makes our life easier

```auto
using GMT

# A test grid
G = grdmath("-R300/900/0/600 -I10 X");

G[G .> 700] .= 5;
G[G .> 600] .= 4;
G[G .> 500] .= 3;
G[G .> 400] .= 2;
G[G .>= 300] .= 1;

# Need to update z_min & z_max
G = grdedit(G);

imshow(G, view=(200,30), fmt=:png)

```

 ![GMTjl_tmp](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/1X/5c099583de334a61df18b20dfe4f0ef287112e7d.png)

---

<div class="post-metadata">

**Author:** ![pwessel](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/pwessel/32/30_2.png) [@pwessel](https://forum.generic-mapping-tools.org/u/pwessel)\
**Post date:** [October 12, 2020, 8:59pm UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/8 "2020-10-12T20:59:38Z")

</div>

It is obviously true that working from Julia, Python, or Matlab adds tremendous extra computational options to such tasks. But that does not mean this is hard in GMT proper, e.g.,

```
gmt grdmath -R300/900/0/600 -I10 X 300 SUB 600 DIV 5 MUL CEIL = bin.nc
gmt grdview t.nc -JX4i -JZ3i -Qi100 -p200/35 -B -Bz -png map

```

 ![map](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/1X/9fde6d9525857e8eeac33e553de5f26829c4c962.png)

---

<div class="post-metadata">

**Author:** ![KRL](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/krl/32/583_2.png) [@KRL](https://forum.generic-mapping-tools.org/u/KRL)\
**Post date:** [October 15, 2020, 2:17am UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/9 "2020-10-15T02:17:47Z")

</div>

Dear all, Thank you for your many devices. I have managed to get this far However still the Velocities are in a continuum form and is best to get ranges of velocities rather than a continuum.

 ![Cumana_Vs30_](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/1X/2470af7b003c280a766ac6927a235698f440f5d7.png)  
and this is the classification I would like to implement on the map to separate velocities y 4 classes with ranges AB \>650; C = 400 to 650 m/s ; D =200 - 400 m/s; D \>200 m/s  
 ![Screen Shot 2020-10-14 at 22.11.51](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/1X/59118c7259ebf5f4c6b1ab355e1911268ad19be8.png)

---

<div class="post-metadata">

**Author:** ![Esteban82](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/esteban82/32/54_2.png) [@Esteban82](https://forum.generic-mapping-tools.org/u/Esteban82)\
**Post date:** [October 15, 2020, 2:43am UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/10 "2020-10-15T02:43:21Z")

</div>

Not sure what would like to do but with grdclip you could classified the grid based on your ranges. Maybe also grdmask to create masks of your selected areas (based on your ranges) to perform math operationes.

---

<div class="post-metadata">

**Author:** ![KRL](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/krl/32/583_2.png) [@KRL](https://forum.generic-mapping-tools.org/u/KRL)\
**Post date:** [October 17, 2020, 4:08am UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/11 "2020-10-17T04:08:51Z")

</div>

Hi, well this is kind of the idea behind, based on the classification trying to get the polygons and then get the Discretized each area. this was made with GMT but I am still getting around with the script to just the idea with my data

 ![Screen Shot 2020-10-17 at 00.02.08](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/1X/55835980bf1e596dbfc84e044e822c3b0ed3536f.jpeg)

#!/bin/sh  
proj=m1:125000 #5/7  
region=139.875/140.11/35.79/35.97 #33=\>23.5, 25=\>18  
gmtset GRID\_PEN\_PRIMARY 0.25p,0.25\_1:0  
gmtset BASEMAP\_TYPE plain  
gmtset ANNOT\_FONT\_SIZE\_PRIMARY 12  
gmtset FRAME\_PEN 1.0p  
gmtset ANNOT\_FONT\_PRIMARY 1  
gmtset GRID\_PEN\_PRIMARY 0.25p,0.25\_1:0  
gmtset BASEMAP\_TYPE plain  
gmtset ANNOT\_FONT\_SIZE\_PRIMARY 15  
gmtset PLOT\_DEGREE\_FORMAT +D  
gmtset TICK\_LENGTH -0.05c  
gmtset ANNOT\_OFFSET\_PRIMARY 0.15c  
gmtset GRID\_PEN\_PRIMARY 0.25p,0.25\_3:0  
echo 1 1 | psxy -Jx1 -R0/2/0/2 -Ss90 -G219/219/219 -K -N -Y1 -X2\>test.eps  
echo 1 0 0 121 100 0 0 121 \> tmp.cpt  
makecpt -Cpanoply -Q -T2/4.2/0.16 \>\>tmp.cpt  
#makecpt -Cpanoply -T1/20000/1000 \>tmp.cpt  
echo B 255 255 255\>\>tmp.cpt  
echo N 0 99 0 \>\>tmp.cpt  
#echo B 0 0 99 \>\>tmp.cpt  
#echo F 99 0 0 \>\>tmp.cpt  
#pscoast -J$proj -R$region -C200/200/255 -K -O -Df -Ia\>\>test.eps  
awk -f jinko.awk kashiwajinko.gmt \> tmp1  
awk -f jinko.awk nagareyamajinko.gmt \>\> tmp1  
awk -f jinko.awk abikojinko.gmt \>\> tmp1  
psxy -J$proj -R$region tmp1 -Ctmp.cpt -M -L -K -A -O -W1,black\>\>test.eps

psxy -J$proj -R$region roadchiba.gmt -M -K -W1,99/99/99 -O \>\>test.eps  
psxy -J -R stations1995.gmt -M -L -A -W6,0/0/0 -O -K -Sc0.2 \>\>test.eps  
psxy -J -R tx.dat -W6,0/0/0 -O -K -Sc0.2 \>\>test.eps  
psxy -J -R railways19552.gmt -M -L -A -W6,49/49/49,… -O -K\>\>test.eps  
psxy -J$proj -R$region lakes22.gmt -M -K -G169/169/169 -O -W1,black\>\>test.eps  
gmtset ANNOT\_FONT\_SIZE\_PRIMARY 18  
gmtset ANNOT\_OFFSET\_PRIMARY 0.28c  
psscale -D21.4/9.7/17/1 -Ctmp.cpt -K -O -L -E0.9\>\>test.eps  
gmtset ANNOT\_FONT\_PRIMARY 0  
gmtset GRID\_PEN\_PRIMARY 0.25p,0.25\_3:0  
gmtset ANNOT\_FONT\_SIZE\_PRIMARY 16  
echo 140.195 35.76 15 0 1 “MC” “[MUs/km^2]” | pstext -J -R -K -O -N -BSWneg0.02a0.04/g0.02a0.02\>\>test.eps  
cp test.eps test2.eps  
awk ‘{if($1==“%%BoundingBox:” || $1==“%%HiResBoundingBox:”){print $1, $2,$3,595,750}else{print $0}}’ \< test2.eps \> test.eps  
convert test.eps test.png  
convert -rotate 90 test.png test.png

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [October 17, 2020, 11:39am UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/12 "2020-10-17T11:39:19Z")

</div>

I think you have to make an effort to explain better what you are trying to achieve. On my side, I’m lost. Your last example has nothing to do with the Vs30 image question.

---

<div class="post-metadata">

**Author:** ![KRL](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/krl/32/583_2.png) [@KRL](https://forum.generic-mapping-tools.org/u/KRL)\
**Post date:** [October 17, 2020, 12:19pm UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/13 "2020-10-17T12:19:07Z")

</div>

Oh okay sorry, I would try to explain myself better, I am trying to make a classification.

My dataset is a .grd file that integrates measurement Velocities (Vs30: Shear wave velocity of the first 30 meters) that were acquired in the field and a topographic model that estimates the Vs30 from the Slope of the relief [https://earthquake.usgs.gov/data/vs30/](https://earthquake.usgs.gov/data/vs30/), I have to integrate both datasets because field measurements have much more resolution than only the topographic model and Also the topographic model itself has more coverage in extension but lower resolution, I have achieved this task, with the Surface command. The generation of an integrated .grd that would include field measurement and the topographic model.

The next step is to study the velocities which I did with the Histogram command and establish ranges that would correspond to classes (from the image above the first one I have published ) I have built 4 classes using the integrated grid of Vs30 (field measurements + topographic approx)

Now I have the ranges for the classes and I need to separate those classes and give One Color to each class so my map represent discrete areas of the Velocities rather than a continuum variation, this means that my map should have 4 colours one for each class

and with these areas or this new grid with 4 classes, I will build polygons.

I hope I have made my best explaining what I am trying to achieve. I am sorry if I confuse you guys. and Thank you for your support.

K

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [October 17, 2020, 1:02pm UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/14 "2020-10-17T13:02:37Z")

</div>

OK, so what is wrong with the solution I gave you?  
Paul’s solution depends on using and expression but you can split it also in multiple grdmath commands.  
Once you have the grid with the classes, `grdcontour` can be used to output the classes polygons.

BTW, if you can build GMT, the `shakemap` branch let you can compute the Vs30 velocities much easier than that USGS recipe, but it only has a online help.

---

<div class="post-metadata">

**Author:** ![KRL](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/krl/32/583_2.png) [@KRL](https://forum.generic-mapping-tools.org/u/KRL)\
**Post date:** [October 17, 2020, 2:32pm UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/15 "2020-10-17T14:32:18Z")

</div>

[quote=“Joaquim, post:14, topic:941”]

I don´t think there is a problem with your solution. The situation is that I don´t know how to use Julia at this moment I have downloaded Julia however GMT wrapper don´t know how to install it and get it working. For the moment I have instal GMT 6 and learning to get things with it  
K

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [October 17, 2020, 2:48pm UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/16 "2020-10-17T14:48:07Z")

</div>

> [@KRL](#):
>
> I have downloaded Julia however GMT wrapper don´t know how to install it

As explained in the GMT.jl page. At the julia command line

 ![image](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/1X/981eb77baa8a1b2baf8da10a4dfdcfe8f1e6279c.png)

Ofc you don’t have to use the Juia solution (although it’s infinitely simpler). You can use `grdmath` too.

---

<div class="post-metadata">

**Author:** ![KRL](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/krl/32/583_2.png) [@KRL](https://forum.generic-mapping-tools.org/u/KRL)\
**Post date:** [October 17, 2020, 3:56pm UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/17 "2020-10-17T15:56:02Z")

</div>

okey I got

 ![Screen Shot 2020-10-17 at 11.54.50](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/1X/1092a4e07d823b9cddf0def590e85a8461cab7a2.png)

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [October 17, 2020, 4:29pm UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/18 "2020-10-17T16:29:42Z")

</div>

OK, so now, in my example, just replace

```auto
G = grdmath("-R300/900/0/600 -I10 X");

```

by

```auto
G = gmtread("/path/to/vs30/grid");

```

at the end, if you want to save the classes grid, do

```auto
gmtwrite("/path/to/classes_grid/classe.grd", G)

```

---

<div class="post-metadata">

**Author:** ![KRL](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/krl/32/583_2.png) [@KRL](https://forum.generic-mapping-tools.org/u/KRL)\
**Post date:** [October 18, 2020, 9:54pm UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/19 "2020-10-18T21:54:39Z")

</div>

Hi,  
I am trying to run the first example you game me so I can understand how GMT julia works or I can make it run and I am getting this ERROR. any suggestions on what should I do?

 ![Screen Shot 2020-10-18 at 17.50.45](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/1X/f01a6f852edbf4b6bdb847c4596c11f4c2a6f168.png)

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [October 18, 2020, 10:17pm UTC](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941/20 "2020-10-18T22:17:19Z")

</div>

Did you install GMT by using the GMT bundle? It wont work in Julia. See [https://github.com/GenericMappingTools/GMT.jl/issues/462](https://github.com/GenericMappingTools/GMT.jl/issues/462)  
MacOS screws the the dependencies

[Next page](https://forum.generic-mapping-tools.org/t/histograms-of-grd-file/941.md?page=2)