How to use -Sj option in gmt plot to plot coupling value of a fault?

AdiMoh08 · August 30, 2024, 2:32pm

Hi GMT team,

I am using the -Sj option in gmt plot to plot interseismic coupling.But somewhere I am confused regarding this option.I dont know what mistake I am doing but I am not getting the results the way I want.

#!/bin/bash
gmt begin plot_rectangles png
gmt basemap -R72/78/31/36 -JM6i -Baf -BWSen
gmt makecpt -Chot -Ic -T-0.3/1/0.01 -H > coupling.cpt
'{print $2, $1, $3}' trialcoupling.txt | gmt plot -Sj314/1/0.5 -Ccoupling.cpt -W0.5p,black
#gmt colorbar -Dx8c/0.5c+w12c/0.5c+h -Baf -By+l"Coupling"
gmt end show

Please find the data file below.Do i need to add anything more in the data file to get the plot?
trialcoupling.txt (52.1 KB)

AdiMoh08 · August 30, 2024, 2:34pm

I dont know why the output is coming like this. I have 1034 values in the table for each fault patch but I think its not taking all the values.
Please help!!

Joaquim · August 30, 2024, 3:09pm

Maybe you want to specify a -R that encloses all your data?

AdiMoh08 · August 30, 2024, 5:21pm

Hi Joaquim,
Even after using -R option , I am getting the same figure.
Somewhat I am trying to get this type of plot (though its generated in matlab) in GMT . I have the centroid long lat value of each element and its corresponding coupling value.
But am not able plot it using -Sj option.How to do it then?

EJFielding · August 30, 2024, 5:45pm

@AdiMoh08 Your “trialcoupling.txt” file only has three columns, which I assume are the latitude and longitude of the center of your fault patch and the coupling value. I would suggest trying the -Sc plotting option first to plot circles at each point to make sure it is reading the file correctly. I think your plotting line is missing the command “awk” at the beginning of the line.

Joaquim · August 30, 2024, 6:21pm

Don’t need this (that intends to swap the first two columns). Just use -:
But the result is what you got (except that the limits should be others). Sorry, but that’s what you have in the file.

AdiMoh08 · August 30, 2024, 6:46pm

Hi
awk i have used in my script. Sorry, i missed it while copying .yes those are the lat lon values of center of each fault patch and 3rd one is the coupling. I guess there is some error in my input file which I am not even able to pick.Could you suggest me how would I plot these fault patches and the variation of coupling values along them?
And I am confused in the -Sj option regarding what it exactly does??

yvonnefroehlich · August 30, 2024, 6:48pm

You can flip the first two columns (lat lon → lon lat) of your input file by using -:i.

I feel the issue is related to your input data: When plotting a very small subregion, there is this output figure, which appears a bit strange to me.

gmt begin subregion_squares png
gmt basemap -R77.3/77.5/32.7/32.9 -JM5c -Baf -BWSen
gmt plot trialcoupling.txt -:i -Ss0.1c -Chot 
gmt end show

AdiMoh08 · August 30, 2024, 6:50pm

Yes, I guess I am missing something in my input option but at the same time I am searching that how do i represent these fault patches showing coupling values in GMT , as shown in the above figure which I have attached just now.

AdiMoh08 · August 30, 2024, 6:54pm

yup, somewhere the input file is missing few more columns i guess which even i am trying to figure out.Not sure if i need the four coordinates of the edges of each patch to plot the coupling value or only centroid lat lon will do.

AdiMoh08 · August 30, 2024, 9:38pm

I modified my input file now and this is the image I am getting now using -Sj option.Why are the squares overlapping? How to correct it?

This is my new input file.
patch_info.txt (55.9 KB)

SCRIPT:

#!/usr/bin/env sh
gmt begin coupling_map png

# Set the map region and projection
gmt basemap -R72/85/25/35 -JM10i -Ba2f2

# Add coastlines for reference
gmt coast -R72/85/25/35 -JM10i -W1/0.5p -N1/0.8p -Df

# Create a color palette for the coupling values

gmt makecpt -Crainbow -T0/1 -Ic > coupling.cpt

# Plot the coupling values at the respective lat/lon positions

awk ‘{print $2, $1, $5}’ patch_info.txt | gmt plot -Sj45/1c/1c -Ccoupling.cpt -W0.5p,black

# Add a color scale bar
#gmt colorbar -Ccoupling.cpt -Dx8c/0.5c+w10c/0.5c+h -Ba0.2f0.1+l"Coupling"

gmt end show

Is it correct?

Joaquim · August 30, 2024, 9:52pm

Because they are too big?

AdiMoh08 · August 31, 2024, 2:53pm

reduced the size now they look like this.How to make them a continuous image showing transition?

EJFielding · August 31, 2024, 4:40pm

If you want to get a plot of the exact rectangles of your model fault colored by the coupling value, you will need to make a file that contains the corners of each patch. As you found, the “-Sj” option plots rectangles at the locations of the centroids. You might be able get those rectangles to approximately connect, but it won’t be a continuous image.

I make plots like this by generating a GMT polygon output file with the slip or coupling on each patch specified with “-Z” lines. This is a part of my slip file:

> -Z0.0 
-98.28774025626007 15.962808020882274 11.07162015337847 
-98.2435904216405 15.948026593547882 11.07162015337847 
-98.22863663544821 15.989809844124903 12.108099529388566 
-98.27279425830342 16.004594623892128 12.108099529388566 
> -Z0.0 
-98.33189691244978 15.977580637505 11.07162015337847 
-98.28774025626007 15.962808020882274 11.07162015337847 
-98.27279425830342 16.004594623892128 12.108099529388566 
-98.31695873048562 16.0193705739944 12.108099529388566 
> -Z0.04733887167672173  
-98.37606037218372 15.992344416809255 11.07162015337847 
-98.33189691244978 15.977580637505 11.07162015337847 
-98.31695873048562 16.0193705739944 12.108099529388566 
-98.36113003396213 16.03413766780622 12.108099529388566 
> -Z0.30485786203874843 
-98.4202306174012 16.007099332205545 11.07162015337847 
-98.37606037218372 15.992344416809255 11.07162015337847 
-98.36113003396213 16.03413766780622 12.108099529388566 
-98.40530815066529 16.048895878719456 12.108099529388566

This is the GMT plotting command (classic mode, so you need to convert to modern mode) that I use:

# fault slip rectangles colored by their -Z values
gmt psxy $file.dat $LLproj $LLreg -C$cpt -L -W0.01c -O -K >> $out

This is an example of the plot:

Joaquim · August 31, 2024, 6:36pm

How do you do that? I do it as well on Mirone (Matlab) for the .fsp and subfault formats but the code is not exactly simple. I wouldn’t mind adding this to GMT.jl

EJFielding · August 31, 2024, 8:00pm

I use the Classic Slip Inversion (CSI) package. It is on GitHub at:

I am afraid it is all in Python, not Julia.

The CSI package has long had a function to write out the patches with their slip values in the GMT polygon format. It includes long, lat, depth or UTM coordinates, depending on a flag. This allows perspective views in GMT.

Joaquim · August 31, 2024, 8:12pm

Thanks. I may have to revisit my old Matlab code to see how hard it is to port it.

mkononets · August 31, 2024, 9:18pm

Maybe this is a stupid idea, but can these data not be gridded? Look very regular. I’ve just got curious about it, so explanations are very welcome.

AdiMoh08 · August 31, 2024, 10:44pm

I tried the same earlier but I dont know why I am not getting the image you got.I am using basic maths to calculate the edges from the centroid of a rectangle but still not getting results like this.

#!/usr/bin/env bash

gmt begin coupling_patches png
gmt basemap -R72/85/25/35 -JM6i -Ba2f1 -BWSen
awk '{
lon_min = $1 - $4/2; where $1=long,$2=lat,
$4=width=0.1 degrees,$5=height=0.05 deg
lon_max = $1 + $4/2;

     lat_min = $2 - $5/2;  
 
     lat_max = $2 + $5/2;  
  
  print lon_min, lat_min;
  print lon_max, lat_min;
  print lon_max, lat_max;
  print lon_min, lat_max;
  print lon_min, lat_min;
  print "> -Z" $3;

}’ patch_uk.txt > patch1.txt
gmt plot patch1.txt -R72/85/25/35 -JM6i -Ccombined.cpt -L -W0.01c
gmt colorbar -Dx6i/0.3i+w5i/0.25i+h -Ccombined.cpt -Bxa0.2f0.1+l"Coupling"
gmt end show
Do you calculate the edges of these patches like this?
I am curious because somewhere I am close but may be making some mistakes.

AdiMoh08 · August 31, 2024, 10:53pm

Hi Joaquim,
I’m currently using MATLAB code to calculate the slip and coupling of subfaults. However, I’m facing an issue with merging neighboring regions that have different convergence rates. The code I’m using can only estimate slip and coupling for one region at a time. If the same fault has different convergence rates, I need to estimate them separately and then attempt to merge them. Since I’m unable to write a code that can merge all the individual blocks with different convergence rates together, I’m looking for a solution from GMT.

Could you suggest me any solution that I can try in matlab?