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

Q&A
22

First note that the patches are not aligned N-S, so they would have to be build from multiple little cells that build up to have that appearance.

Second, the fault planes are not horizontal so a grid would have to somehow represent also the fault dip and orientation.

And third, at least the data that I know, produced by agencies type USGS when inverting earthquakes waveforms, produce vector data, not grids.

23

I’m probably not a big help to you. The cases I dealt with was in ingesting seismic inversion solutions and use them to calculate Okada deformations for tsunami modeling. I do that in Mirone where I also produce maps of the horizontal projection of those fault patches (similar to what you seem to be trying to obtain) but I did not calculate the patches. “Just” read those complicated formats.

I haven’t visited this for a long time but code that reads and plots “Sub Fault Format” is in

24

Thanks a lot for explaining the above.

Yes sure, I noticed very well the patches aren’t aligned N-S, but most of them look regular on the map. Can a grid be defined on an oblique projection or is it just the map view is deceitful?

If such an irregular (with regard to N-E) grid can be defined, other parameters can be specified on the same grid, the dip and the orientation. This may again be a stupid idea, useless or irrelevant as I don’t understand the data and much of the context. If it is so don’t hesitate pointing this out (and explanations are of course much appreciated if provided). I’m just a sediment/benthic biogeochemist who happened interested about GMT, GIS and related subjects.

Can a grid be defined on an oblique projection at all, regardless the topic’s data?

Thanks.

25

The CSI package has several Python data structures (objects) that contain the information about the patches, including the list of corners for each patch. It calculates the corners when it creates the patches in the first place, and then keeps that information.

26

Understood. Manual calculation is not suitable. Since I have been conducting all the analysis in MATLAB, I need to examine the (CSI) package. Thank you for guiding me in the right direction, as I was unable to pinpoint where the issue was.

27

I need to understand the code and maybe it will be a help to me.Thank you for sharing.

28

Hi Eric,
Sorry for bothering you again, but I’ve written a code to calculate the corners of each patch. Even after using the -Z option, the figure I’m getting remains the same. I am attaching the file.
patch_polygons.txt (118.4 KB)
When running this in GMT, I’m still getting gaps, but the pattern it shows is correct. However, it’s not continuous.

coupling_patches
coupling_patches1959×1409 157 KB

29

What plotting command are you using? You should not use the -Sj option with the polygon files. That option is for point files.

If you are using the polygon option, then the problem is the way that you generated the polygons for each patch. If you take the maximum and minimum of the latitude and longitude for the patch, then you will get rectangles that are like the ones you show. You need to figure out the patch corners that were rectangles in the model fault coordinate system but are rotated in the lat-long coordinate system.

You should go back to the program that generated your model fault patches in the first place and use that to calculate the patch corners. It will be extremely difficult to get from the patch centroids.

30

#!/usr/bin/env bash

gmt begin coupling_patches png

# Set the region and projection
gmt basemap -R72/80/30/35 -JM6i -Ba2f1 -BWSen
gmt plot patch_polygons.gmt -R72/80/30/35 -JM6i -Ccombined.cpt -L -W0.01c
# Add a color bar for coupling values
#gmt colorbar -Dx6i/0.3i+w5i/0.25i+h -Ccombined.cpt -Bxa0.2f0.1+l"Coupling"

gmt end show

This the script I am using

31

I attach one of the Python programs that I used to generate the patches for CSI. In this case, I am making a model fault on the Mexican subduction zone so I used the Slab 1 database to get the depths of the slab surface and the location of the trench. Then I create a rectangular model fault that has the same strike as the trench. This is a step that you need to do for your model fault to get the patch corners.

MakeFaultPlane-Ometepec-v1.txt (5.8 KB)

32

No, rotated grids are not supported in GMT.

33

Not sure I understand this answer (or my question was not properly formulated)

We know there are e.g. UTM grids, that is, Cartesian grids with UTM assigned as a coordinate system. Can the same not be done, just using e.g. oblique Mercator, in principle?

I imagine there can be an issue defining the corresponding Cartesian coordinate system, but e.g. projected coords for oblique Mercator should be “the same” as for the non-oblique, regardless the oblique angle?

34

Sorry for the short answer.
Yes, I don’t see a reason why we cannot create grids in coordinates that have been projected in Oblique Mercator and thus obtain a Cartesian grid in that referencing system . But that grid would not have stored the information grid is oblique and the oblique angle. It is in this sense that I say that we don’t support rotated grids. ~30 years ago I used to work with grids that supported this notion in order to process/map magnetic grid acquired along the the magnetic declination but have not seen them ever since. Without this information and software support one cannot use those grids to mix, for example, with geographic data, e.g. plotting coastlines, etc…

35

Thanks again for the explanations!