search

ThomasLecocq / msnoise-tomo

The Ambient Noise Tomography Plugin for MSNoise
http://www.msnoise.org
European Union Public License 1.1
21 stars 15 forks source link

observed travel time calculation #22

Open dylanmikesell opened 4 years ago

dylanmikesell commented 4 years ago

I think there is an error in the t_obs calculation in ANSWT.py. This value is currently computed as

dist=(np.hypot(x1-x2, y1-y2)) t_obs = dist/Vg1

where dist is being computed from x1,x2,y1,y2, which are in units of degrees. I spoke with Aurelien today and he suggested using the rows in the G-matrix. This approach worked, but only after converting x1,x2,y1,y2 to a local cartesian grid. I tried many different things prior to this, and I could never get a good solution using units of degrees. Right now I am working on a grid that is 200m on either side.

I propose this solution.

# Compute the observed times in whatever units the grid is. dist = GG.toarray().sum(axis=1) # distance in units of the grid t_obs = dist / Vg1 # [??] arrival time data

But as I said, this has only worked for cartesian coordinates, and I think it has to do with src/mk_MatPaths.c. This is again just using the hypotenuse to compute distances in each grid cell and build the G-matrix.

Regardless of units, t_obs and t_calc should match in the way they are computed and right now they are not.

I also propose that we convert station coordinates to a local cartesian system inside of ANSWT.py and use a local grid. I can work on this. We can write the output back to actual coordinates, but I can't get the code to work right now without converting units...I even tried tracking down which units are causing the problem because Aurelien thought G is normalized and distance units would not matter. But I cannot find the problem yet.

Thoughts?

ThomasLecocq commented 4 years ago

Hi Dylan,

Indeed, in case the coordinates are DEG it's not good (except close to the equator)... It's been a long long time I haven't looked at the code. If you have a solution, please indeed PR it. Actually there are already "interstation distances" calculated elsewhere in the code (and in the prepare_ccf file) where the real distance is computed (using msnoise.api.get_interstation_distance)

thanks for all this

dylanmikesell commented 4 years ago

I have moved from a grid in units of degree to a grid in units of kilometers. This uses pyproj to convert to UTM coordinates and then kilometers from meters. Now units of data match G matrix and strange issues related to model values and data values appear to be gone. Commit is 3285c04. Will do a pull request soon.