Problem with oblique mercator projection

[cossatot]

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

[jploveless]

Thank you for pointing this out, Richard. The oblique Mercator routine changed when we adapted Blocks to work on global models, and one routine that is automatically done is that any segment that spans 0º longitude is split into two segments, with the split endpoint at the prime meridian. I think you're right that faultobliquemercator.m may not work as a general purpose projection tool, and I recall some issues in trying to use Matlab's Mapping Toolbox and the USGS projection reference in handling these special cases of spanning the prime meridian. Given the other routines within GetElasticPartials.m, the elastic dislocation solution is output correctly, even if the projected coordinates look funky in the intermediate steps, although there are still issues with stations that lie at exactly 0º longitude (see correction below).

My test, run after generating the variables from the code you posted, is this:

s.lon = s.lon+1e-8; % Perturb station longitudes by ~1 mm
G = GetElasticPartials(f, s);
figure
plot([f.lon1 f.lon2], [f.lat1 f.lat2], 'r')
hold on
quiver(s.lon, s.lat, G(1:3:end, 1), G(2:3:end, 1)); % Just using the first column of the full partials matrix, corresponding to unit strike slip
title('Okada strike-slip via GetElasticPartials.m');
xlabel('deg lon');
ylabel('deg lat');

and it produces:

Hopefully this addresses the issue, although it's a bit unsatisfactory that the projection routine doesn't work as a standalone tool. Seems like moving to proj4 could be a good solution!

Thanks, Jack

[cossatot]

Thanks Jack, I'm glad that the issue is handled properly by the code.

It's always frustrating when you're simplest-possible test case is actually a weird edge-case...

Richard