Hi all,
Not directly GMT related, but with the bright minds here, maybe I can get some help.
I have a nice grid defined by 1degx1deg in the ED50 datum. If I were to transform this grid to WGS84, the coordinates would become very ugly; e.g the point 21E,68N (ED50) → 21.9923E,67.9932N (WGS84). I want to avoid this.
One possible solution would be to cheat - just ‘move the grid’; instead of treating the grid as ED50, lie, and ‘say’ it’s in WGS84. This would shift the grid (by a relatively small amount, ~100 m or so depending on where you are), and we do not clutter up the nice coordinates. Any points/polygons defined on the ED50 grid, may be transformed to WGS84 without any issues, so this is not a big deal.
My question is, how much will the grid move, and how can I best visualize this?
I’ve attached two screenshots; I have checked how far the grid moves in two places of the area. This is just to illustrate that the amount of shift depends on where you are. The red lines are the ED50 grid and the black lines are the WGS84 grid.
Any thoughts, tips or tricks are very welcome! If this is a terrible solution, I would like to know!
Presumably over your small 1x1 degree area, the bulk of the shift is a constant, plus deviations. Maybe create two grids via grdmath: One with longitude and one with latitudes, then dump to tables and form lon,lat table you run via mapproject to do the datum shift, then reformat back to grids and subtract the original grid. Now you have delta lon and delta lat on a grid. You can scale these to meters assuming a local flat Earth, and now you can plot these grids to see how the distortion/shift works. Maybe take out the mean lon and mean lat shift first and look at the deviations. Perhaps just plotting the grid as is (as you say) but shift it by teh mean offset is good enough for government work.
Thanks for the suggestion Paul, I’ll follow your bread crumbs.
Bonus question: Why is the output of gmt mapproject -Qd sent to stderr?
With external calls itwould be bad if mapproject returned taht freeform text to the caller I think. So we implemented it as a verbose thing.
Alright, I now have delta lon and delta lat on a grid.
Could you please advise, how I should proceed to properly scale the deltadegrees to meter? As a rough test I converted the deltadegrees to meter by multiplying with 60*1852. Then I applied the distance formula. This seems to give a too big shift, compated to the figures I attached in my first post.
I didnt quite describe the area I’m dealing with good enough; nevermind the grid that I spoke about. The area is approximately between 2E - 36E and 55N to 74N. So it’s quite large.
(Within this area, there are rows and columns [i.e. graticules] of 1degx1deg defining a grid)
Maybe it would be simpler if I worked with projected coordinates(?)
I don’t understand well. Do you want to convert Spherical Degrees to Kilometers?
Maybe you need to use:
gmt grdmath deg2km
I was thinking you would simply scale your delta_lat grid by 111113 to get meters and for delta_lon you need to scale by Y COSD 111113 MUL, with 111113 a number from memory for meters in a degree at Equator (or along parallel), but you can possibly use DEG2KM instead in both places.
Thank you both for great ideas. I’ll work some more on this. Forgot there was a deg2km.
Paul, I’m very content that I (to a degree [pun intended]) tried what you described. Just forgot the COSD factor for the longitude.
Thanks a bunch to Paul and Esteban for helping me out.
For anyone interested, this is what the final script looks like:
#!/bin/bash
rm 2a-ed50.txt 2b-konv_fra_ed50_til_wgs84.txt &>/dev/null
for ((x=-14;x<=38;x++));
do
for ((y=56;y<=85;y++));
do
printf "%i %i\n" $x $y >> 2a-ed50.txt
#coordinate conversion between datums from (64) European 1950 to (219) WGS 1984
printf "%i %i\n" $x $y | gmt mapproject -T64/219 >> 2b-konv_fra_ed50_til_wgs84.txt
done
done
#make deltadeg files
paste 2a-ed50.txt 2b-konv_fra_ed50_til_wgs84.txt | awk '{print $1-$3}' > 3a-delta_long.txt
paste 2a-ed50.txt 2b-konv_fra_ed50_til_wgs84.txt | awk '{print $2-$4}' > 3b-delta_lat.txt
#make xyz of deltasets
paste 2a-ed50.txt 3a-delta_long.txt > 4a-grid_deltalong.txt
paste 2a-ed50.txt 3b-delta_lat.txt > 4b-grid_deltalat.txt
#grid deltasets
gmt xyz2grd -G5a-grid_deltalong.nc $(gmt info -I1 4a-grid_deltalong.txt) -I1 4a-grid_deltalong.txt
gmt xyz2grd -G5b-grid_deltalat.nc $(gmt info -I1 4b-grid_deltalat.txt) -I1 4b-grid_deltalat.txt
#convert to meter
gmt grdmath Y COSD 111120 MUL 5a-grid_deltalong.nc MUL = 6a-grid_deltalong_METER.nc
gmt grdmath 111120 5b-grid_deltalat.nc MUL = 6b-grid_deltalat_METER.nc
#I tried this first, but this gave too big numbers (about twice as big, I think). due to no long scaling?
#gmt grdmath -fg 5a-grid_deltalong.nc DEG2KM 1000 MUL = 6a-grid_deltalong_METER.nc
#gmt grdmath -fg 5b-grid_deltalat.nc DEG2KM 1000 MUL = 6b-grid_deltalat_METER.nc
#use pythagoras to calculate total shift d=[(deltax_meter)^2+(deltay_meter)^2]^0.5
gmt grdmath 6a-grid_deltalong_METER.nc 2 POW 6b-grid_deltalat_METER.nc 2 POW ADD 0.5 POW = 7-tot-disp.nc
#clean up
rm gmt.history
And now a bit of some interesting observations.
In addition to mapproject, I also did a test with proj's cs2cs. This uses 7 parameters to do the datum transform, instead of 3 (mapproject).
Take a look at the time it takes to complete for mapproject and cs2cs:
With mapproject
#!/bin/bash
for ((x=-14;x<=38;x++));
do
for ((y=56;y<=85;y++));
do
printf "%i %i\n" $x $y | gmt mapproject -T64/219
done
done
real 1m4.507s
user 0m55.031s
sys 0m9.737s
With cs2cs:
#!/bin/bash
for ((x=-14;x<=38;x++));
do
for ((y=56;y<=85;y++));
do
printf "%i %i\n" $x $y | cs2cs -f "%.2f" +proj=longlat +ellps=intl +towgs84=116.641,56.931,110.559,-4.327,-4.464,4.444,3.520 +no_defs +to +proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs
done
done
real 0m1.153s
user 0m1.187s
sys 0m0.321s
cs2cs does this operation about 60x faster than mapproject.
1 Like
Yes, GMT has a startup cost and you are calling it again again with one data point. I would rewrite your script to build the input file first, then call mapproject once. While cs2cs will still be faster I bet it wont be as bad as this worst case scenario!
Programming (or scripting…) is fascinating. Only limited by imagination. You can do brilliant things - and - extremely stupid things. My script represent the latter.
Thanks Paul. Again.
#!/bin/bash
rm temp &> /dev/null
for ((x=-14;x<=38;x++));
do
for ((y=56;y<=85;y++));
do
printf "%i %i\n" $x $y >> temp
done
done
gmt mapproject -T64/219 temp
real 0m0.205s
user 0m0.134s
sys 0m0.041s
#!/bin/bash
rm temp &>/dev/null
for ((x=-14;x<=38;x++));
do
for ((y=56;y<=85;y++));
do
printf "%i %i\n" $x $y >> temp
done
done
cs2cs -f "%.2f" +proj=longlat +ellps=intl +towgs84=116.641,56.931,110.559,-4.327,-4.464,4.444,3.520 +no_defs +to +proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs temp
real 0m0.078s
user 0m0.047s
sys 0m0.028s
