Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

NLLOC lat/lon conversion #72

Open
saeedsltm opened this issue Jan 12, 2018 · 14 comments
Open

NLLOC lat/lon conversion #72

saeedsltm opened this issue Jan 12, 2018 · 14 comments
Labels

Comments

@saeedsltm
Copy link

Hi,
Could you please help me how to set gk2lonlat function correctly when using our own coordinate system in NLLOC?

Regards

@megies
Copy link
Owner

megies commented Jan 12, 2018

Well.. you should just replace the function such that it converts from the coordinates you use in NonLinLoc (this is assuming that in NonLinLoc you use TRANS NONE, I think) to WGS84. If you have a coordinate system that has an EPSG code you can do it pretty much the same as in the original gk2lonlat function.

@saeedsltm
Copy link
Author

Thanks, OK but when we use something like this "TRANS LAMBERT WGS-84 35.50 50.50 32.0 38.0 0.0" which is already WGS, then what would be gk2lonlat? In this case our reference is WGS84!

@megies
Copy link
Owner

megies commented Jan 12, 2018

OK, so currently, I guess you should just put

def gk2lonlat(x, y):
    return x, y

@saeedsltm
Copy link
Author

saeedsltm commented Jan 12, 2018 via email

@megies
Copy link
Owner

megies commented Jan 12, 2018

Attachments to email replies are not retained by Github..

@saeedsltm
Copy link
Author

sorry, by the way it did not work. i think the problem is i need to convert x/y from Cartesian coordinate system (source and station relative to map center in TRANS line of NLL control file) to lon/lat. So the output coordinate system is WGS-84 but i dont have any idea about the input coordinate system.

@saeedsltm
Copy link
Author

saeedsltm commented Jan 15, 2018

Finally i found a solution and here is the change list:
1- write new gk2lonlat function in util.py like:
def gk2lonlat(lon0,lat0,x,y): lon0*=ones_like(x) lat0*=ones_like(y) wgs84_geod = Geod(ellps='WGS84') new_lon, new_lat, _ = wgs84_geod.fwd(lon0,lat0,rad2deg(arctan2(y,x)),sqrt(x**2+y**2)*1e3) return new_lon, new_lat
2- and change the following line:
data[0], data[1] = gk2lonlat(data[0], data[1])
to:
data[0], data[1] = gk2lonlat(lon0, lat0, data[0], data[1])
3- add the following lines to obspyk.py's line 2191:
# read reference coordinate with open(files['summary']) as f: for l in f: if "TRANSFORM" in l: lat0 = float(l.split()[5]) lon0 = float(l.split()[7])
4- and change the following lines:
4-1) line:
lon, lat = gk2lonlat(x,y)
to:
lon, lat = gk2lonlat(lon0,lat0,x,y)
4-2) line:
data = readNLLocScatter(PROGRAMS['nlloc']['files']['scatter'], self.widgets.qPlainTextEdit_stderr)
to:
data = readNLLocScatter(lon0, lat0, PROGRAMS['nlloc']['files']['scatter'], self.widgets.qPlainTextEdit_stderr)
Thanks.

@megies
Copy link
Owner

megies commented Jan 15, 2018

Ah, I think I understand what the problem is now.. basically in your case we need to parse the GEOGRAPHIC line in the nonlinloc hypocenter-phase output file, instead of the HYPOCENTER line that I usually use because I have to convert coordinates manually after the nonlinloc run (all of our models are in local XYZ coordinate systems).

Compare two flavors of the hypocenter output file:

I'll try to come up with a simple switch for the configuration file, so that this can be done without hacking the source code in the future.. maybe even later today or this week

@megies megies added the task label Jan 15, 2018
@megies
Copy link
Owner

megies commented Jan 15, 2018

Oh, just saw your python repos.. nice! Need to give PyVelest and PyGAVel a try at some point.. 🚀 :-)

@saeedsltm
Copy link
Author

Yes, actually i made a change to use GEOGRAPHIC but it didn't resolve my problem because stations are also needed to change their positions.
About the repos, they may need some minor changes i think ;).

@megies
Copy link
Owner

megies commented Jan 15, 2018

but it didn't resolve my problem because stations are also needed to change their positions.

Hmm.. I don't understand, can you clarify? I don't see how station coordinates are ever handled in something else than lon/lat..

@saeedsltm
Copy link
Author

saeedsltm commented Jan 15, 2018 via email

@megies
Copy link
Owner

megies commented Jan 16, 2018

Ok, so what I would do is add a config option for nlloc that can either be..

  • None which means that no transform will be done and GEOGRAPHIC line will be used (can you please confirm that the scatter is also in lat/lon in your case?)
  • int, in which case this will be used as an EPSG code (with an additional option for a scaling factor to apply -- I usually have kilometers in nlloc even though I use a meter-based EPSG coordinate system)
  • str, in which case this will be interpreted as a function that is imported and used for coordinate conversion

The latter two cases would then use HYPOCENTER line and the first option the GEOGRAPHIC line.

@saeedsltm
Copy link
Author

scatters are in x,y relative to geographical center in introduced in TRANS line. So they need to be converted according to reference coordinate system.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants