search

GeoscienceAustralia / eqrm

Automatically exported from code.google.com/p/eqrm
Other
5 stars 4 forks source link

Rupture distance calculation is incorrect #30

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago 
What steps will reproduce the problem?
1. Set DISTANCE_LIMIT = 0.000001 in distance functions. (I will check this in 
soon)
2. run test_all.py

What do you see instead?
A distance unit test fails (that I just added)

Please use labels and text to provide additional information.

This error is due to the site locations being set to a local co-ord system with 
the trace start as the origin.  The function rupture_xy (in 
distance_functions.py) is expecting the mid point of the
    rupture or the mid point of the rupture trace as the origin (The exact origin used needs to be determined.)

So a solution to this could be pass the rupture_centroid_lat/long to the 
distance functions.

More unit tests need to be added though.  I added a simple one, but a more 
complex one with dip, depth and azimuth is needed.

Original issue reported on code.google.com by duncan.g...@gmail.com on 9 Mar 2012 at 5:37

GoogleCodeExporter commented 9 years ago

Original comment by b...@girorosso.com on 14 Mar 2012 at 5:01

GoogleCodeExporter commented 9 years ago 
distance_functions.Rupture_xy is based on definition in 
Documentation/Distance2.doc

The code looks pretty accurate to what is documented.

There are two issues here. The first one is explained here:

Key pseudocode for Rupture_xy
1. x = abs(x)
2. y = abs(y)
3. l = lengths / 2
4. w = widths / 2
5. x = where(x < l, l, x)
6. y = where(y < w, w, y)
7. x = x - l
8. y = y - width
9. distance = sqrt(x^2 + y^2 + z^2)

The idea behind lines 5 and 6 is that if the x and y coordinates are inside the 
rupture plane then set the values to the edge.

This falls over, however when x and y are not within the rupture plane, but 
abs(x) and abs(y) are. This is shown in the failing test 
test_distance_functions.test_Rupture. 

e.g.

site 7 -> x = -0.1, y = -0.1
and
length = 0.4, width = 10.0 and dip = 90.0

1. dip = 90 degrees, width is projected to 0 on the 2D plane, so abs(y) > width
2. abs(x) is less than length / 2 and so it is set to this (0.2)

This means that the distance for site 7 would be calculated as the same as site 
4 as the coordinates are adjusted to be the same (abs(x) = 0.2, abs(y) = -0.1)

The assertion message confirms this:

AssertionError: Expected Rx=
[[ 35.14084638]
 [ 33.33753404]
 [ 35.14084638]
 [ 11.11251135]
 [  0.        ]
 [ 11.11251135]
 [ 11.11251135]
 [  0.        ]
 [ 11.11251135]]
got
[[  3.51407126e+01]
 [  3.33375340e+01]
 [  3.51407126e+01]
 [  1.11124436e+01]
 [  1.00000000e-06]
 [  1.11124436e+01]
 [  1.11124944e+01]
 [  1.00000000e-06]
 [  1.11124944e+01]]

Original comment by b...@girorosso.com on 15 Mar 2012 at 3:35

GoogleCodeExporter commented 9 years ago 
The second issue:

Section 4.2 of the EQRM manual states:

Fault distance (Rrup): Closest distance between rupture plane and site of
interest.

This is not what distance_functions.Rupture_xy implements. The relevant lines 
in the comment above:

3. l = lengths / 2
4. w = widths / 2
7. x = x - l
8. y = y - width

This reads as calculate the shortest distance to the edge of the plane where l 
and w are the x and y coords of the rupture centroid.

test_distance_functions.test_Rupture tests the shortest distance to the rupture 
plane i.e.

sites 
1 - 3: distance to end
4 - 6: distance to centroid
7 - 9: distance to start

Original comment by b...@girorosso.com on 15 Mar 2012 at 3:43

GoogleCodeExporter commented 9 years ago 
Duncan,

Let me know how you want to proceed given the two issues highlighted above.

Cheers,
Ben

Original comment by b...@girorosso.com on 15 Mar 2012 at 3:44

GoogleCodeExporter commented 9 years ago 
After discussion with Duncan, this distance function will be removed and a 
distance calculation provided by Trevor Allen will be used.

Original comment by b...@girorosso.com on 15 Mar 2012 at 11:14

GoogleCodeExporter commented 9 years ago 
I've turned off the test_Rupture test, so test_all is clean for the install 
Robbie is doing.

Original comment by duncan.g...@gmail.com on 16 Mar 2012 at 6:02

GoogleCodeExporter commented 9 years ago 
Kaklamanos et al. (2011) contains a well described methodology for calculating 
Rrup. 

I have implemented this and it passes the vertical strike test added by Duncan 
with the following changes:

- removing the distance limit in the Joyner-Boore distance function
- converting the x,y values of the sites in the test to lat,lon based on the 
trace start lat,lon
- passing through rupture centroid x,y to the distance function in the test
- removing the inconsistent deg per km multiplications

Further tests to be added:
- non-vertical strike tests
- Joyner-Boore distance tests (Kaklamanos et al. uses Rjb and there are 
currently no tests for this)

Original comment by b...@girorosso.com on 20 Mar 2012 at 4:08

GoogleCodeExporter commented 9 years ago 
Revision 1008 includes the implementation of Rrup from Kaklamanos et al. 
(2011), along with unit tests and an Excel worksheet used to produce test data.

Work to be completed:
- Unit tests for Joyner-Boore
- Refactor references to the existing Rrup function to use the new one

Original comment by b...@girorosso.com on 21 Mar 2012 at 3:12

GoogleCodeExporter commented 9 years ago 
Revision 1009 adds
- unit tests for Joyner-Boore
- add depth to top of rupture plane to the distance function interface and 
ensures Kaklamanos_Rupture uses this rather than hypocentral depth

Original comment by b...@girorosso.com on 21 Mar 2012 at 5:16

GoogleCodeExporter commented 9 years ago 
Revision 1010 obsoletes the existing Rupture distance function and uses 
Kaklamanos_Rupture for this calculation. 

Scenarios and long scenario standard results updated to reflect this change.

Original comment by b...@girorosso.com on 21 Mar 2012 at 11:50

GoogleCodeExporter commented 9 years ago 
Reopening and assigning to me as David B has found this implementation does not 
resolve the issue highlighted, where for certain sites the calculation 
overestimates distance.

Original comment by b...@girorosso.com on 22 Mar 2012 at 5:56

GoogleCodeExporter commented 9 years ago 
The issue highlighted by David B was as follows
- length, width and depth are specified as integers
- the calculation of length and width of rupture centroid used in the x and y 
values also uses integer values:
l = lengths/2
w = cos_dip*widths/2
- in this case length was 25 where length/2 = 12 not 12.5

The fix is ensures these are calculates as floating point values

l = lengths/2.0
w = cos_dip*widths/2.0

This is checked in in revision 1015

Original comment by b...@girorosso.com on 22 Mar 2012 at 11:21