jploveless / blocks Goto Github PK

I processed around 150 SLC images for estimating the crustal deformation (only descending orbit available) using PS-InSAR technique. How to give this input to Blocks?

Hi,
First of all thank you so much for keeping this wonderful software in open-source. I have doubt regarding intraplate slip rate estimation using GPS network. Is it possible, if it is please kindly provide some suggestions to calculate the slip rate in intraplate regions using your software.

While running the sample file in Block software, I'm getting below error. Please check it once.

Blocks('/home/rakesh/Documents/Blocks-master/Kachchh1/command')
Parsing input data...Error using fgets
Invalid file identifier. Use fopen to generate a valid file identifier.

Error in fgetl (line 32)
[tline,lt] = fgets(fid);

Error in ReadCommand (line 57)
line = fgetl(infile);

Error in Blocks (line 5)
Command = ReadCommand(commandFile); % Read command file

Can InSAR data be used in the Blocks?

Dear members,
I already mentioned that I'm working on Intraplate region. I want to separate the Intraplate region to two different blocks because of GPS velocities (Northward stations showing southward motion and vise versa; attached figure). Then based on GPS velocities I wanna calculate the two different Euler poles for two blocks. Is it correct procedure to calculate the slip rate of different faults presented in two different blocks?

Your kind suggestions will enhance my knowledge.

(Dumka et al., 2023; GPS and InSAR derived evidences of intra‐basin stress and strike‐slip tectonics in the vicinity of 2001 (M7.7) earthquake, Kachchh, western India)

Hello,

When I was trying to learn from the examples you gave, the function “StartsWith ” seems to be missing. Is there a problem with my use? I don't know what happened. And I can't find this function in all the files. Please tell me what should I do.

Thanks,
sunset

Dear members,
I used the Blocks for an interseismic coupling inversion of a strike slip fault. The result of strike slip rate from relative block motion in Mod.patch file is shown in figure below. It is noteworthy that certain fault segments exhibit positive slip rates while others display negative values. I found that this variation in sign is related to the south-dipping or north-dipping of faults. However, the definition of left lateral or right lateral strike slip should be independent of the fault dip. So I'm confused about such reversal of strike slip direction.

Another question is about the dip slip rate from relative block motion in Mod.patch file. The inversion result of dip slip is unrealistically large, especially in the element whose dip angle is close to 90°. Is there any method to constrain the dip slip rate in the patch inversion?

Thanks,
Jeff

Hi,
I am using the most recent version of Blocks from this repository, and it seems like the BlocksForward problem is not working. I get this error.

In an assignment A(:) = B, the number of elements in A and B must be the same.

Error in BlockLabel (line 218)
stl(stl == 0) = ext;

Error in BlocksForward (line 41)
[Segment, Block, Station] = BlockLabel(Segment, Block, Station);

Error in visualize_fwd_blocks (line 11)
V=BlocksForward(EE,NN,outdir);

However, I checked with the version of Blocks in brendanjmeade's repository and that seems to be working fine. Can you help me figure out what is wrong?

Thanks,
Rishav

Hi,

I'm using some functions from Blocks to benchmark/test a project I'm working on. Thanks for making the code available and FOSS for such as great project.

One of the issues that I ran into was that the projection done in faultobliquemercator.m does not seem to be correct, which causes the displacements from the Okada dislocations to be funky.

I have checked the projection against the Proj4 library and do not run into this issue. I've also checked the code against the USGS reference and it looks correct, so I am not sure where the issue is. Nonetheless for my project I am just using Proj4 (I'm not working in Matlab), so it doesn't impact me but I thought I would let yall know.

To test, I'm using a script, which will be below, that takes a vertical strike-slip fault embedded in a grid of points in lon-lat coords, projects things, and does the Okada dislocation for unit strike slip. The results are here:

Here are the points (blue) and fault (red) in geographic coordinates:

Here is everything in the projection from faultobliquemerc.m:

Here is the resulting Okada solution:

If I use the projected coordinates from Proj4 (in Python, script to follow) and then bring these back into the matlab code, this is what I get for the projected geometries:

And the Okada results look a lot more like what we'd expect:

Here is the matlab code:

clear;

f = [-1. 0. 1. 0.];

[slon, slat] = meshgrid(-8.:0.5:8., -8.25:0.5:8.25);

slon_ = reshape(slon, 1, []).';
slat_ = reshape(slat, 1, []).';

s = [slon_ slat_];

nu = 0.25;

[f, s] = ProjectSegCoords(f,s);

f.dip = 90.;
f.lDep = 20.;
f.bDep = 0.;

[strike, L, W, ofx, ofy, ofxe, ofye, ...
               tfx, tfy, tfxe, tfye]        = fault_params_to_okada_form(f.px1, f.py1, f.px2, f.py2, deg_to_rad(f.dip), f.lDep, f.bDep);
           
 [ves, vns, vus,...
  ved, vnd, vud,...
  vet, vnt, vut] = okada_partials(ofx, ofy, strike, f.lDep, deg_to_rad(f.dip), L, W, 1., 0., 0., s.fpx, s.fpy, nu);
 

figure
plot(s.lon, s.lat, '.')
hold on
plot([f.lon1 f.lon2], [f.lat1 f.lat2], 'r')
title('Coords in Lon-Lat');
xlabel('deg lon');
ylabel('deg lat');

figure
plot(s.fpx, s.fpy, '.')
hold on
plot([f.px1 f.px2], [f.py1 f.py2], 'r')
title('projected coords from faultobliquemerc.m');
xlabel('x (m)');
ylabel('y (m)');

figure
plot([f.lon1 f.lon2], [f.lat1 f.lat2], 'r')
hold on
quiver(slon_, slat_, ves, vns)
title('Okada strike-slip via faultobliquemerc.m');
xlabel('deg lon');
ylabel('deg lat');

%%%%% using external coords from Proj4

proj_coords = csvread('obl_coords.csv');
proj_fault_coords = csvread('obl_fault_coords.csv');
ps = s;

ps.fpx = proj_coords(:,1);
ps.fpy = proj_coords(:,2);

pf = f;
pf.px1 = proj_fault_coords(1,1);
pf.px2 = proj_fault_coords(2,1);
pf.py1 = proj_fault_coords(1,2);
pf.py2 = proj_fault_coords(2,2);

[pstrike, pL, pW, pofx, pofy, pofxe, pofye, ...
 ptfx, ptfy, ptfxe, ptfye]        = fault_params_to_okada_form(pf.px1, pf.py1, pf.px2, pf.py2, deg_to_rad(pf.dip), pf.lDep, pf.bDep);
           
[pves, pvns, pvus,...
  pved, pvnd, pvud,...
  pvet, pvnt, pvut] = okada_partials(pofx, pofy, pstrike, pf.lDep, deg_to_rad(pf.dip), pL, pW, 1., 0., 0., ps.fpx, ps.fpy, nu);

figure
plot(ps.fpx, ps.fpy, '.')
hold on
plot([pf.px1 pf.px2], [f.py1 f.py2], 'r')
title('projected coords from Proj4')
xlabel('x (m)');
ylabel('y (m)');

figure
plot([pf.lon1 pf.lon2], [pf.lat1 pf.lat2], 'r')
hold on
quiver(ps.lon, ps.lat, pves, pvns)
title('Okada strike-slip via Proj4');
xlabel('deg lon');
ylabel('deg lat');

And here is the Python code to do the projection:

import numpy as np
from pyproj import Proj, transform


lon1 = -1.
lon2 = 1.
lat1 = 0.00001
lat2 = -0.00001

lons, lats = np.meshgrid(np.arange(-8., 8.5, 0.5), np.arange(-8.25, 8.75, 0.5),
                         indexing='ij')

lons, lats = lons.ravel(), lats.ravel()

wgs84 = Proj(init='epsg:4326')

fm_str = "+proj=omerc +lat_1={lat1} +lon_1={lon1} +lat_2={lat2} +lon_2={lon2}".format(
    lat1=lat1, lon1=lon1, lat2=lat2, lon2=lon2)

obl_m = Proj(fm_str)

x, y = transform(wgs84, obl_m, lons, lats)

fx, fy = transform(wgs84, obl_m, [lon1, lon2], [lat1, lat2])

np.savetxt('obl_coords.csv', np.vstack([x,y]).T, delimiter=',')
np.savetxt('obl_fault_coords.csv', np.vstack([fx,fy]).T, delimiter=',')

There is a very small difference in the lats of the fault because Proj4 errors out with a purely 90-striking fault.

I've also tested in some different geometries and the matlab code still does not produce the expected results.

Maybe I'm doing something wrong, though!

Thanks,
Richard