please add DIST to SAC headers

[alexhutko]

STAL/STLO and EVLA/EVLO are in the SAC headers, but DIST is not, which many users may find useful.

[krischer]

No problem. Something to keep in mind: everything within Instaseis happens in a perfectly spherical earth so the distance would also be calculated on a sphere. I think this is acceptable if it is documented.

[sceylan]

I think that's a very good idea. And I think it should be a read-only parameter. Probably you would also consider this, but I just wanted to point it out to be sure.

Second, I would introduce ellipticity as a parameter at a method interface (if possible). SAC does use spherical planet for distance calculation; therefore tomography people like myself often use their own codes for distance calculations. I think this would be a nice utility to have.

[krischer]

@sceylan FYI: for now this would only be happening for some of the server routes: http://instaseis.net/server.html

Instaseis currently does otherwise not bother with I/O at all and delegates that responsibility to ObsPy. We could totally add a trace.stats.sac header to the objects returned from Instaseis. In that case these headers would be used by ObsPy if written out in the SAC format. Can you please open a separate issue for that?

Another question: If SAC uses a spherical Earth for the distance calculation - are the coordinates spherical as well or defined on the WGS84 ellipsoid?

[alexhutko]

SAC stores them as geographic and converts to geocentric. https://ds.iris.edu/files/sac-manual/commands/sss.com/traveltime.html

[krischer]

That's for theoretical travel times which are always calculated in a spherical earth due to using the tau-p method. In that sense it is kind of similar to instaseis.

I assume you are fine with a spherical/geocentric distance?

[alexhutko]

IRIS assumes geographic latitudes so as long as they get conveted to geocentric, that's fine. I think this is what you are saying.

[CTrabant]

According to the notes in @sac-101.6a/src/ucf/distaz.c@:

@note Calculations are based upon the reference spheroid of 1968 and are defined by the major radius (RAD) and the flattening (FL).

So not even WGS84. I cannot find a modern reference for a 1968 ellipsoid, looks like same radius at the SOUTH AMERICAN 1969 ellipsoid, but different flattening ( http://www.arsitech.com/mapping/geodetic_datum/). Anyway. It's not spherical, the Earth it not a sphere after all.

[krischer]

Done in e5d22b9bb130bc2d4d1aff2cdbf705f5ef3a917c. Please test it.

Just to clarify: Syngine currently takes coordinates in WGS84 as input values. Your middleware converts these to geocentric coordinates which get passed to the Instaseis server. Instaseis always treats the Earth as a sphere. Thus the distance is also calculated on a sphere as that is the only always correct thing.

Now in your middleware you could reset the dist header field to the ellipsoidal difference between source and receiver. That would be more correct in terms of physical distance but not correct in terms of the distance instaseis sees.

I propose to keep the spherical distance - it just has to be properly documented.

[CTrabant]

Providing geodetic coordinates with a distance on a sphere feels like an unexplainable mismatch, users will make assumptions so documentation of quirks is non-ideal.

I ended up doing the simple calculations in the IRIS middleware. Now we are populating: DIST, AZ, BAZ and GCARC headers using the geodetic coordinates provided by the user and the WGS84 ellipsoid.

Ideally the middleware would not need to do this, maybe in can be incorporated into Instaseis in the future.

[krischer]

You are totally right - my bad. I forgot that SAC header fields are always in WGS84 and instaseis already converts them to the correct coordinates... Fixed in 0461b42927c1d88e6c583f3942d8ccc9ebfc12c9

If you use this in syngine, please also install geographiclib:

$ pip install geographiclib

ObsPy has a fallback implementation but the geographiclib one is faster and probably also more accurate for the edge cases. ObsPy will use geographiclib if it is available.

I also added the required additional header fields. I'm not sure the GCARC one if fully correct. This particular one uses geocentric coordinates to get the greatcircle arc distance which I think is correct in this case.

[CTrabant]

Nice. Pulled update and it's running at the DMC.