Ladies and Gents,
I need to find all points where the terrain penetrates the surface of a cone placed on its tip on the reference point.
In this vertical cut mockup the terrain is grey, the blue star is my reference point (lon, lat, elev), the light yellow area is the cone with a given slope (1:20) and the orange areas are the ones I’m interested in.
How would you approach this problem? My experimentation so far was not very successful. I appreciate your nudging me in the right direction.
Are your data gridded or arbitrarily spaced?
That’s a line-of-sight problem. We don’t have anything (that I know) in GMT for that. GDAL actually has a program that in principle could be called (or added to) grdgdal but when I once looked at it I got a bit disappointed with possibilities. Don’t remember more now.
I would imagine that you could create an artificial grid for the cone surface, difference the two grids, and NaN the values <0 (based on the interpretation that you’re looking for the points in orange), all using grdmath. If you want an ASCII file with the points above the cone rather than a grid, you could export with grd2xyz and -s.
With no guarantees, because I haven’t tested or checked this, these are generally the operations that I was thinking of. Feel free to let me know the reasons for this failing, in the case that it doesn’t work.
gmt grdmath -R$R -I30s $lon $lat SDIST 1000 MUL 20 DIV $elev ADD data.nc SUB 0 GT 0 NAN = cone_pts.nc
Hi @maxrjones unfortunately I didn’t get lucky with your suggestion. Here is what I tried based on your idea:
gmt grdcut -R7/10/49/51 @earth_relief_30s -Grelief_fra.nc
No questions here so far.
gmt grdmath -R7/10/49/51 -I30s 008:34:13.64 50:01:59.90 SDIST = dist_from_fra.nc
I’m unsure what the unit of my calculated values is. This should give me the distance to my supplied coordinate in 30s increments or doesn’t it (-I30s)? If it is in half arc minute increments and I want Meters I need to multiply it by 3704 (1’ = 1 NM = 1852 m) and divide it by 20 to get my slope.
gmt grdmath dist_from_fra.nc 3704 MUL 20 DIV = slope.nc
I hope so far my attempt makes sense? Here I get caught by my inadequate knowledge of grdmath:
gmt gmtmath relief_fra.nc slope.nc GT = terrain.nc
I want all values which are above slope.nc and the rest gets ignored. Instead I get a segfault.
gmtmath [ERROR]: NetCDF variable lat has different dimension (240) from others (360)
gmtmath [WARNING]: Long input record (4131 bytes) was truncated to first 4094 bytes!
ERROR: Caught signal number 11 (Segmentation fault: 11) at
0 libgmt.6.dylib 0x00000001028a167e GMT_gmtmath + 8718
1 ??? 0x0000000000000000 0x0 + 0
Stack backtrace:
0 libgmt.6.dylib 0x0000000102686a02 sig_handler + 578
1 libsystem_platform.dylib 0x00007fff67ccb5fd _sigtramp + 29
2 ??? 0x0000000000000000 0x0 + 0
3 libgmt.6.dylib 0x00000001026aa52b GMT_Call_Module + 2331
4 gmt 0x000000010267b998 main + 2568
5 libdyld.dylib 0x00007fff67ad2cc9 start + 1
6 ??? 0x0000000000000007 0x0 + 7
I’m pretty sure I made some mistakes along the way which lead to this mess but I need you kindly pointing me towards it. Thank you for your help!
Hi @KristofKoch, the SDIST always returns distances in km:
SDIST 2 1 Spherical (Great circle geodesic) distance (in km) between nodes and stack (A, B)
So if you want meter then 1000 MUL needs to be added.
Also, you cannot use gmtmath (expects tables) to read grids hence it crashes. So stick to grdmath until you have a true/false (1/0) decision for each node, then convert false to NaN, possibly multiply that with your grid if you want the z-values, and then save. Then use grd2xyz -s to skip the NaN nodes to yield the final table.
Hi @pwessel thank you for pointing my errors out – those who can read the docs have a clear advantage. That gmtmath thing was an autocompletion SNAFU. I ran into another snag: @earth_relief_30s comes with pixel node registration … switching it to @earth_relief_30s_g did the trick.
Thanks again for your help!
Your questions and problems are so intriguing and interesting.
If you feel like it, please share some of your results (if it results in e.g. a fancy map). I would love to see what you’re up to!