Instantly share code, notes, and snippets.

Embed
What would you like to do?
set terminal pngcairo font "arial,8" size 640,350
set output 'simple_distance.png'
set multiplot layout 1, 2
r(d) = d*pi/180
dist_hav(x,y,x0,y0) = 2*asin(sqrt(sin(r(y-y0)/2)**2 + cos(r(y0))*cos(r(y))*sin(r(x-x0)/2)**2))*6371
dist_sim(x,y,x0,y0) = sqrt(r(y-y0)**2 + (cos(r(y0))*r(x-x0))**2)*6371
error(x,y,x0,y0) = 100*abs(dist_hav(x,y,x0,y0)-dist_sim(x,y,x0,y0))/dist_hav(x,y,x0,y0)
set view map
set palette cubehelix negative
set contour
set isosamples 200, 200
set cntrparam level discrete 5
set title "Relative error at 89°N"
set xrange [-.05/cos(r(89)): .05/cos(r(89))]
set yrange [89-.05:89+.05]
splot error(x,y,0,89) notitle with pm3d, \
dist_hav(x,y,0,89) nosurface notitle
set title "Relative error at 65°N"
set xrange [-.05/cos(r(65)):.05/cos(r(65))]
set yrange [64.95:65.05]
splot error(x,y,0,65) with pm3d notitle, \
dist_hav(x,y,0,65) nosurface notitle
unset multiplot
@Tobii42

This comment has been minimized.

Tobii42 commented Nov 15, 2018

This is awesome! You can become even more accurate if you use the formula for equirectangular distance calculation from https://www.movable-type.co.uk/scripts/latlong.html :

earth_radius = 6371
dist_eqi_x(lng0, lat0, lng1, lat1) = r(lng1 - lng0) * cos(r(lat1 + lat0)/2)
dist_eqi_y(lat0, lat1) = r(lat1 - lat0)
dist_eqi(lng0, lat0, lng1, lat1) = earth_radius * sqrt(dist_eqi_x(lng0, lat0, lng1, lat1)**2 + dist_eqi_y(lat0, lat1)**2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment