Hey I am wondering if anyone could look at my code and verify that I am heading in the right direction. Here is a summary of my project thus far. I have a contour plot of elevation points of a known latitude and longitude section. Next I made contour profiles of desired locations on the contour map, explained here:
http://www.idlcoyote.com/ip_tips/image_profile.html
Now I am trying to find the distance between longitude and latitude points using the vincenty formula form here: http://en.wikipedia.org/wiki/Great-circle_distance. I have completed this and I am getting logical answers I am just wondering if it is accurate or a correct.
Pro ariel
openr, lun, 'arial.txt',/get_lun
data=dindgen(824,914)
readf, lun, data
close,lun
; CONTOURING ARIEL
s=Size(data,/dimensions)
lons=scale_vector(dindgen(s[0]), 270, 360)
lats=scale_vector(dindgen(s[1]), 0, -90)
nlevels = 850
step = (Max(data) - Min(data)) / nlevels
levels = IndGen(nlevels) * step + Min(data)
device, decomposed=0
loadct, 1
window, 1, retain=2
;contour, data, lons, lats, background=cgcolor('white'), color=cgcolor('black'), xstyle=1, ystyle=1, title='Ariel', ytitle='Latitude', xtitle='Longitude', levels=levels
; obtaining the x1,x2,y1,y2 values needed for interpolation and contour profile
mousebutton=!mouse.button & !mouse.button=0
print,'right click to quit'
while !mouse.button ne 4 do begin
wset,1
print,'click on the starting point.'
cursor,x1,y1, 3,/data
if x1 gt s(0) || y1 gt s(1) then begin
print,'first point is off the image.'
goto,skip
endif
print,'click on the ending point.'
cursor,x2,y2,3,/data
if x2 gt s(0) || y2 gt s(1) then begin
print,'second point is off the image.'
goto,skip
endif
wset,1
plots,[x1,x2],[y1,y2], color=150
;MAKING IMAGE PROFILES
window, 2, retain=2
;npoints= ABS(x2-x1+1) > ABS(y2-y1+1)
npoints=150
xloc= x1+(x2-x1)*findgen(npoints)/(npoints-1)
yloc= y1+(y2-y1)*findgen(npoints)/(npoints-1)
profile=interpolate(data, xloc, yloc)
plot, profile
;CALCULATING THE DISTANCE
lond1=lons[x1]
lond2=lons[x2]
latd1=lats[y1]
latd2=lats[y2]
lon1=lond1*!DTOR ; converting to radians
lon2=lond2*!DTOR
lat1=latd1*!DTOR
lat2=latd2*!DTOR
dlon=abs(lon2-lon1)
dlat=ABS(lat2-lat1)
r=578.9
; vincenty formula
a=sqrt((cos(lat2)*sin(dlon))^2 + (cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dlon))^2)/(sin( lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(dlat))
c= atan(a)
d=r*c
print, d
wset,1
skip:
endwhile
; user quites
!mouse.button=mousebutton
stop
END
Also I am wondering if anyone has any advice on how use my new distance value and plot it against my contour profile. Right now my x axis for my contour profile is just the number of npoints I used, which is 150 but I would like it to represent the distance.
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
