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
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?
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
Also it would be great if I could see the histograms of the studied dataset. this would defenely would be great
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.
Thank you for your advice I will give it a try ! and be back!
Thanks !
An example on how the Julia wrapper makes our life easier
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)
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
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.
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.
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
#!/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
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.
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/, 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
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.
[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
As explained in the GMT.jl page. At the julia command line
Ofc you donât have to use the Juia solution (although itâs infinitely simpler). You can use grdmath too.
OK, so now, in my example, just replace
G = grdmath("-R300/900/0/600 -I10 X");
by
G = gmtread("/path/to/vs30/grid");
at the end, if you want to save the classes grid, do
gmtwrite("/path/to/classes_grid/classe.grd", G)
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?
Did you install GMT by using the GMT bundle? It wont work in Julia. See https://github.com/GenericMappingTools/GMT.jl/issues/462
MacOS screws the the dependencies
