#!/usr/local/bin/davinci -f
###############################################################################
#
# Davinci Program to Compute corrections due to plane motions for
# AVIRIS (but priniples will apply to most scanners.
#
# Roger N. Clark, U.S. Geological Survey
# MS 964, Box 25046 Federal Center
# denver, CO 80225
# rclark@speclab.cr.usgs.gov
#
#
# REVISIONS
# $Log: cmd.davinci.avirisgeom.2,v $
# Revision 1.2 97/10/21 15:38:20 15:38:20 raymond (Ray Kokaly)
# Changed code to compute ranges correctly for vhead (yaw calc.
#
# Revision 1.1 97/10/21 11:24:20 11:24:20 raymond (Ray Kokaly)
# Initial revision
#
# Revision 11/4/97 (version 4) includes along track velocity
# variations and roll position at start of data acquisition.
# - Roger N. Clark 11/4/97
#
# Revision 11/10/97 (version 5) includes do some computations in
# double precision, write output to specpr files.
# add user titles to plots, command line specpr file name, grid size
# and mean topographic elevation.
# - Roger N. Clark 11/10/97
#
# Revision 11/11/97 (version 6) complete version 5 additions
# add bow computation in case flight is not straight
# - Roger N. Clark 11/11/97
#
# add equatorial and polar radii, computations for cross track scan velocity
# and aircraft elevation
# - Roger N. Clark 11/12/97
#
# version 8: refine and speed up yaw computation, fix small things
# - Roger N. Clark 11/13/97
#
###############################################################################
verbose=0 # be silent
#
# READ COMMAND LINE OPTIONS
#
args=$ARGC # number of command line arguments
echo("command line arguments = " + args)
if (args < 3) {
printf("%s\n","ERROR: insufficient input")
Usage = $0 + " Engineering_File Navigation_File title_for_plots_and_specpr"
Usage = Usage + " [-plot] [-specpr file_name]"
Usage = Usage + " [-grid grid_size_in_meters]"
Usage = Usage + " [-topo height_in_meters]"
printf("%s\n","Usage: " + Usage)
printf("%s\n","exit 1")
}
if (args < 3) { # bug. this must be separate or
# last if statement exists before printing.
#exit (1) # no exit because davinci exits always when it sees it
}
file11 = $1 # Engineering File (required)
file12 = $2 # Navigation File (required)
title1 = $3 # title for plots and specpr files (required)
do_plot = 0 # default: no plots
specprout = 0 # default: no specpr output
specpr_file="specpr_file.sp"
grid = 16.0 # default grid is 16 meters
topo = 1500. # default mean topographic height in meters
# really need digital elev model for each pixel
gvcal = 1.000 # velocity calibration (empirically derived)
# if not equal 1.0, then the ground velicity measurement of
# the aircraft is not perfectly calibrated.
# note: might be 0.993
exit_flag = 0 # flag to exit when > 0
if (args > 3) { # find other options
for (i = 4; i <= args; i = i+1) {
#echo("DEBUG: $ARGV[" + i + "] = " + $ARGV[i])
if($ARGV[i] == "-plot") { # do plots
do_plot=1
printf("%s\n","enabling plots")
}
if($ARGV[i] == "-specpr") {
specprout = 1
if ( args >= i+1) { # should be specpr file name
i1 = i+1
txx=$ARGV[i1]
specpr_file=txx
printf("%s\n","will do output to specpr file " + specpr_file)
i = i1
} else {
printf("%s\n","ERROR: no specpr file name after -specpr")
printf("%s\n","exit 1")
exit_flag = 1
}
}
if($ARGV[i] == "-grid") {
if ( args >= i+1) { # should be grid size
i1 = i+1
txx=$ARGV[i1]
grid=atof(txx)
printf("%s","setting grid size = ")
echo(grid)
i = i1
} else {
printf("%s\n","ERROR: no grid size after -grid")
printf("%s\n","exit 1")
exit_flag = 1
}
}
if($ARGV[i] == "-topo") {
if ( args >= i+1) { # should be grid size
i1 = i+1
txx=$ARGV[i1]
topo=atof(txx)
printf("%s","Mean Topographic elevation = ")
echo(topo)
i = i1
} else {
printf("%s\n","ERROR: no elevation after -topo")
printf("%s\n","exit 1")
exit_flag = 1
}
}
}
}
if ( exit_flag > 0) {
printf("%s\n","NEED TO EXIT---RESULTS will not be correct")
#exit(exit_flag) # no exit because davinci exits always when it sees it
}
printf("%s\n","Mean Elevation = " + topo)
printf("%s\n","Output grid size = " + grid)
ttmp="plot, specpr title= " + title1
printf("%s\n", ttmp)
if (specprout == 1) {
echo ("Output to specpr file enabled: " + specpr_file)
}
#
# READ NAV DATA FROM NAV FILES
#
echo ("reading navigation and enginerring data files")
a=read(filename=file11)
b=read(filename=file12)
#
# SET AVIRIS PARAMETERS
#
pcenter = 307 # center of line
pixels = 614 # number of cross-track pixels
pedge = 307 # center to edge distance in pixels
linestime = 12.0 # scan lines per second
ptop = 0.8745 # pixel to pixel distance in mrad (get from file 5)
# file 5 = .87; aviris engineers say .8745 microsec = int time
ifov = 0.99 # IFOV in mrad (get from file 5)
#
# SET CONSTANTS
#
pi = 3.14159265
rad = 360./(2.*pi) # degrees per radian
# the following earth radii from the book: Astrophysical Quantities, by Allen, 1976
earthpolarradius= 6356779 # meters # NOTE: get exact value
eartheqradius = 6378164 # meters # NOTE: get exact value
# note: mean earth radius = 6371030 meters
#
# EXTRACT AIRCRAFT ATTITUDE DATA FROM NAV DATA
#
roll = b[19,,] # roll in degrees (+ = right)
pitch = b[18,,] # pitch
airvelocity = b[17,,] # speed in m/s # this is zero in some aviris data
altitude = b[16,,]
vspeed = b[15,,]
gvelocity = b[14,,] # ground speed
t = b[13,,]
hours = b[6,,]
minutes = b[7,,]
sec = b[8,,]
latitude = double(b[9,,]) # degrees
longitude = double(b[10,,]) # degrees
#
dhours = hours + minutes/60.0 + sec/3600.0
td = dim(t)
nlines = td[2,,] # number of lines in data set
#
#
x=moment(latitude)
meanlatitude = x[3]
echo ("mean latitude = " + meanlatitude)
# find equivalent earth radius at mean latitude and mean topographic elevation
sinmlat = sin(meanlatitude/rad)
cosmlat = cos(meanlatitude/rad)
e1 = earthpolarradius + (eartheqradius - earthpolarradius) * sinmlat
eradius = e1 + topo
meterdeg = 2.0 * pi * eradius / 360.0 # meters per arc-degree
#
# ADJUST TRUE HEADING FOR 0 = 360 DEGREES
# if heading varies between the two
#
a=moment(t)
if (a[1] < 10.0 && a[2] > 350.0) { # min and max heading bounce around 0, 360 degrees
# so make about 0 (going negative: 359 = -1)
echo ("Heading is north, values range " + a[1] + " to " + a[2] + ", ADJUSTING")
k = nlines+1
i=1
while (i 350.0) {
t[1,i,1] =t[1,i,1] - 360.0
#echo (i+ " " + t[1,i,1])
}
i = i +1
}
}
aa = moment(t)
tmean = aa[3]
# because north-south flights, longitude is a large number with
# round-off in single precision numbers and gplot doen't see those last
# digits well, remove the whole degree and look at relative longitude
r=moment(longitude)
rlong1=double(short(r[2]))
rlong=longitude-rlong1
#
# write data to specpr files
#
if (specprout == 1) {
sptitle=title1 + " Pitch (degrees)"
write(translate(pitch,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Roll (degrees)"
write(translate(roll,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " True Air Speed m/s" # this is zero in some aviris data
write(translate(airvelocity,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Pressure Altitude (meters)"
write(translate(altitude,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Verticle Velocity (m/s)"
write(translate(vspeed,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Ground Speed (m/s)"
write(translate(gvelocity,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " True Heading (degrees)"
write(translate(t,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Latitude (degrees)"
write(translate(latitude,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Longitude (degrees)"
write(translate(longitude,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Rel Long (deg+" + short(rlong1) + ")"
write(translate(rlong,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Time (hours)"
write(translate(dhours,y,z),filename=specpr_file,title=sptitle,type=specpr)
}
#
# SET UP PLOT
#
if(do_plot == 1) {
plot("set xlabel 'Line Number'") # st horixontal axis label
plot("set format y '%10.5f'") # set vertical axis data format
plot("set time") # put date and time on plot
plot("set data style lines") # plot data as lines
plot("set nokey") # turn key off
plot("set terminal postscript landscape") # output in postscript format
plot("set terminal postscript landscape color \"Times-Roman\" 12")
# plot("replot")
# plot("set output \"|lp work_plot.gnu\"") # local tmp file lp or lpr
# plot("replot")
# plot("pause 10") could be needed if you use a loop
# plot("set terminal x11")
# plot("replot")
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'Pitch'")
plot("set format y '%10.3f'")
plot("set output 'pitch.ps'")
plot (pitch)
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'Roll'")
plot("set format y '%10.3f'")
plot("set output 'roll.ps'")
plot (roll)
plot("set ylabel 'Meters / sec.'") # set vertical axis label
plot("set title 'true Air Speed'") # this is zero in some aviris data
plot("set format y '%10.0f'")
plot("set output 'air.speed.ps'")
plot (airvelocity)
plot("set ylabel 'Meters'") # set vertical axis label
plot("set title 'Pressure Altitude'")
plot("set format y '%10.0f'")
plot("set output 'pressure.altitude.ps'")
plot (altitude)
plot("set ylabel 'Meters / sec.'") # set vertical axis label
plot("set title 'Vertical Velocity'")
plot("set format y '%10.3f'")
plot("set output 'vertical.velocity.ps'")
plot (vspeed)
plot("set ylabel 'Meters / sec.'") # set vertical axis label
plot("set title 'Ground Speed'")
plot("set format y '%10.0f'")
plot("set output 'ground.speed.ps'")
plot (gvelocity)
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'True Heading (degrees)'")
plot("set format y '%10.3f'")
plot("set output 'true.heading.ps'")
#plot (b[13,,]) # uncorrected heading
plot (t)
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'Longitude (degrees)'")
plot("set format y '%10.5f'") # set vertical axis data format
plot("set output 'longitude.ps'")
plot (longitude)
plot("set ylabel 'Degrees'") # set vertical axis label
ptitl = "'Relative Longitude (degrees plus " + short(rlong1) + ")'"
plot("set title" + ptitl)
plot("set format y '%7.5f'") # set vertical axis data format
plot("set output 'longitude-rel.ps'")
plot (rlong)
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'Latitude (degrees)'")
plot("set format y '%10.5f'") # set vertical axis data format
plot("set output 'latitude.ps'")
plot (latitude)
plot("set title 'Universal Time (hours)'")
plot("set format y '%10.5f'") # set vertical axis data format
plot("set output 'time.ps'")
plot(dhours)
}
################
# now do some computations.
#############################################################################
####### relative downtrack position from ground speed
echo ("computing relative downtrack position from ground speed")
# do it double precision
ddistance = double(gvelocity) * double(0.0) # create downtrack distance array with zeros
sdistance = gvelocity * 0.0 # create downtrack sample distance array with zeros
# this is the uniform sampling distance.
i=2
while ( i <= nlines ) {
ddistance[1,i,1] = ddistance [1,i-1,1] + double(gvelocity[1,i,1]*gvcal/linestime)
sdistance[1,i,1] = grid * (i - 1) # uniform grid
i = i + 1
}
# offsetdistance is the error in downtrack position relative to uniform grid
offsetdistance = (ddistance - sdistance) / grid
if (specprout == 1) {
sptitle=title1 + " Line Offst frm Grid (pixels)"
write(translate(offsetdistance,y,z),filename=specpr_file,title=sptitle,type=specpr)
}
if(do_plot == 1) {
plot("set ylabel 'Lines'") # set vertical axis label
plot("set format y '%10.4f'")
plot("set title 'Line Offset from Grid due to Ground Velocity (pixels)'")
plot("set output 'grid.velocity.offset.ps'")
plot(offsetdistance)
}
#############################################################################
####### predicted longitude and latitude ####
echo ("predict longitude and latitude based on a straight line between ends of run")
along1=b[10,1,1] # begining longitude
along2=b[10,nlines,1] # ending longitude
alat1=b[9,1,1] # begining latitude
alat2=b[9,nlines,1] # ending latitude
dslong = along2 - along1
dslat = alat2 - alat1
slong = dslong / nlines # slope of longitude vs line
slat = dslat / nlines # slope of latitude vs line
# mean heading from end points of lat,long
xmtp = slat
if (abs(dslat) < 0.1e-5) {
vmeanhead = 90.
} else {
vmeanhead = rad * atan(dslong*cosmlat/dslat)
}
# if latitude is decreasing, we are flying south
# so adjust because atan fcn only workd from -90 to + 90 degrees
if (dslat < 0) { # flying torward south, not north
vmeanhead = vmeanhead + 180.
}
echo ("mean velocity heading from GPS = " + vmeanhead + " degrees")
vmeanheada = clone( vmeanhead, y=nlines)
offlong = along2 - slong * nlines
offlat = alat2 - slat * nlines
predlong = longitude # make a predicted longitude array
predlat = latitude # make a predicted latitude array
i = 1
while (i <= nlines) {
predlong[1,i,1] = offlong + slong * i
predlat[1,i,1] = offlat + slat * i
i = i + 1
}
diflong = longitude - predlong
diflat = latitude - predlat
echo ("smoothing difference longitude and difference latitude")
diflongs = ysmooth(diflong,150) # smooth over 150 channels
diflats = ysmooth(diflat, 150) # smooth over 150 channels
if (specprout == 1) {
sptitle=title1 + " predict Longitude (degrees)"
write(translate(predlong,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " predlat (degrees)"
write(translate(predlat,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Lat -predict (degrees)"
write(translate(diflat,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Long -predict (degrees)"
write(translate(diflong,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " smoothed Lat -predict (degrees)"
write(translate(diflats,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " smoothed Long -predict (degrees)"
write(translate(diflongs,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Mean velocity GPS heading (degrees)"
write(translate(vmeanheada,y,z),filename=specpr_file,title=sptitle,type=specpr)
}
if(do_plot == 1) {
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'Longitude (degrees)'")
plot("set format y '%10.5f'") # set vertical axis data format
plot("set output 'longitude.ps'")
plot (longitude,predlong)
#
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'Latitude (degrees)'")
plot("set format y '%10.5f'") # set vertical axis data format
plot("set output 'latitude.ps'")
plot (latitude,predlat)
plot("set title 'Latitude difference (latitude - predicted, degrees)'")
plot("set output 'latitude-dif.ps'")
plot (diflat)
plot("set title 'Longitude difference (longitude - predicted, degrees)'")
plot("set output 'longitude-dif.ps'")
plot (diflong)
plot("set title 'smoothed Latitude difference (latitude - predicted, degrees)'")
plot("set output 'latitude-difsmooth.ps'")
plot (diflats)
plot("set title 'smoothed Longitude difference (longitude - predicted, degrees)'")
plot("set output 'longitude-difsmooth.ps'")
plot (diflongs)
plot("set title 'Mean velocity heading for line from GPS (degrees)'")
plot("set output 'meanheadgps.ps'")
plot (vmeanheada)
}
#############################################################################
### now interpolate binned values
echo ("interpolating longitude to every scan line")
### Longitude and Longitude are read out once every second, and
# the values repeated vor 12 scan lines until the next value is
# acquired (it should be interpolated, so we must do it here)
### Longitude
along = predlong # get array defined
last = along[1,1,1]
aline = along # create another array
twonum = along[1,1:2,1] # array 2 elements along for two y values
twoxvl = along[1,1:2,1] # array 2 elements along for 2 x values
i = 2
aline[1,1,1]=1
ifirst = 1
while ( i <= nlines ) {
xtmp = along[1,i,1]-last
if ( xtmp < 0.0 ) xtmp = -1.0 * xtmp
if (xtmp < 0.15e-3) {
aline[1,i,1]=i
#echo(i + " " + xtmp)
ilast=i
i = i + 1
} else {
#echo("interpolating at " + i)
aline[1,i,1]=i
twoxvl[1,1,1] = ifirst
twoxvl[1,2,1] = i
twonum[1,1,1] = b[10,ifirst,1]
twonum[1,2,1] = b[10,i,1]
xnew=interp(twonum,twoxvl,aline[1,ifirst:ilast])
xnew
#now set the new interpolated values into the arry to be interpolated
m = ifirst
k = 1
while ( m <= ilast ) {
along[1,m,1] = xnew[1,k,1]
k = k + 1
m = m + 1
}
# now reset to look next position
ifirst = i
last = along[1,i,1]
aline[1,1,1]=i
}
}
if (specprout == 1) {
sptitle=title1 + " interpolated Longitude (degrees)"
write(translate(along,y,z),filename=specpr_file,title=sptitle,type=specpr)
}
if(do_plot == 1) {
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'interpolated Longitude (degrees)'")
plot("set format y '%10.5f'") # set vertical axis data format
plot("set output 'longitude-interp.ps'")
plot (longitude,along)
}
#############################################################################
### Latitude: interpolate to every scan line
echo ("interpolating latitude to every scan line")
alat = b[9,,] # get array defined
last = alat[1,1,1]
aline = alat # create another array
twonum = alat[1,1:2,1] # array 2 elements alat for two y values
twoxvl = alat[1,1:2,1] # array 2 elements alat for 2 x values
i = 2
aline[1,1,1]=1
ifirst = 1
while ( i <= nlines ) {
xtmp = alat[1,i,1]-last
if ( xtmp < 0.0 ) xtmp = -1.0 * xtmp
if (xtmp < 0.5e-4) {
aline[1,i,1]=i
#echo(i + " " + xtmp)
ilast=i
i = i + 1
} else {
if(i - ilast == 1) {
#echo ("ilast = i, skipping at" + i)
ilast = i
alat[1,i,1] = b[9,i,1]
i = i + 1
} else {
echo("interpolating at " + i)
aline[1,i,1]=i
twoxvl[1,1,1] = ifirst
twoxvl[1,2,1] = i
twonum[1,1,1] = b[9,ifirst,1]
twonum[1,2,1] = b[9,i,1]
xnew=interp(twonum,twoxvl,aline[1,ifirst:ilast])
#xnew
#now set the new interpolated values into the arry to be interpolated
m = ifirst
k = 1
while ( m <= ilast ) {
alat[1,m,1] = xnew[1,k,1]
k = k + 1
m = m + 1
}
# now reset to look next position
ifirst = i
last = alat[1,i,1]
aline[1,1,1]=i
}
}
}
if (specprout == 1) {
sptitle=title1 + " interpolated Latitude (degrees)"
write(translate(alat,y,z),filename=specpr_file,title=sptitle,type=specpr)
}
if(do_plot == 1) {
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'interpolated Latitude (degrees)'")
plot("set format y '%10.5f'") # set vertical axis data format
plot("set output 'latitude-interp.ps'")
plot (latitude,alat)
}
#############################################################################
# now compute yaw
echo ("Computing yaw")
m=dim(alat)
nlines=m[2,,] # number of lines
ndel = 50 # channel delta
ndel1 = ndel - 1
ndel2 = ndel / 2
nchandel = nlines - ndel2
navg = 50 # number of channels to average (make even number)
navg2= navg / 2
vhead = double(alat) # create velocity heading array same size as lot, or long
dellat= double(alat)
dellon= double(along)
i = ndel2 + 1
while (i <= nchandel) {
#echo ("Computing yaw, line " + i)
####### average blocks navg in size
m = i - ndel2 - navg2 ### first average
m2= m + navg
xx = moment(alat[1,m:m2,1])
avlat1 = xx[3] # average for interval m to m2
xx = moment(along[1,m:m2,1])
avlon1 = xx[3] # average for interval m to m2
m2 = i + ndel2 + navg2 ### second average
m = m2 - navg
xx = moment(alat[1,m:m2,1])
avlat2 = xx[3] # average for interval m to m2
xx = moment(along[1,m:m2,1])
avlon2 = xx[3] # average for interval m to m2
dellat[1,i,1] = avlat2 - avlat1
dellon[1,i,1] = (avlon2 - avlon1) * cosmlat # adjust diff longitude for mean latitude
# lines getting closer near poles
xtmp = dellat[1,i,1]
if (xtmp < 0) xtmp = -1.0 * xtmp
if ( xtmp < 0.1e-5 ) {
vhead[1,i,1] = 90. # do not divide by zero
} else {
vhead[1,i,1] = atan(dellon[1,i,1]*cosmlat/dellat[1,i,1])*rad
}
#echo ("Value of vhead = " + vhead[1,i,1])
i = i + 1
}
echo ("Computing velicity heading: filling upper array")
i = i - 1
while (i <= nlines) {
vhead[1,i,1] = vhead[1,i-1,1]
dellat[1,i,1] = dellat[1,i-1,1]
dellon[1,i,1] = dellon[1,i-1,1]
i = i + 1
}
m = 1
echo ("Computing velicity heading: filling lower array")
while (m <= ndel2) {
vhead[1,m,1] = vhead[1,ndel2 + 1,1]
dellat[1,m,1] = dellat[1,ndel2 + 1,1]
dellon[1,m,1] = dellon[1,ndel2 + 1,1]
m = m + 1
}
###was: yaw = t - float(vhead) # or is it vhead - t ?
yaw = t - float(vmeanhead) # or is it vhead - t ?
if (specprout == 1) {
sptitle=title1 + " latitude derivative avg=50"
write(translate(dellat,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " longitude derivative avg=50"
write(translate(dellon,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Velocity Heading (degrees)"
write(translate(vhead,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Yaw (degrees)"
write(translate(yaw,y,z),filename=specpr_file,title=sptitle,type=specpr)
}
if(do_plot == 1) {
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'Latitude derivative (degrees, avg=50)'")
plot("set format y '%10.3f'")
plot("set output 'dellat.ps'")
plot(dellat)
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'Longitude derivative (degrees, avg=50)'")
plot("set format y '%10.3f'")
plot("set output 'dellon.ps'")
plot(dellon)
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set title 'Velocity Heading (degrees)'")
plot("set format y '%10.3f'")
plot("set output 'velocity.heading.ps'")
plot(vhead)
plot("set format y '%10.3f'")
plot("set ylabel 'Degrees'") # set vertical axis label
plot("set format y '%10.3f'")
plot("set title 'Yaw (degrees)'")
plot("set output 'yaw.ps'")
plot(yaw)
}
#############################################################################
## OK now compute relative pixel offsets
# because yaw of 1 degree represents only 5 pixels (line direction)
# offset, there is negligable sample offset, so ignore.
# sample offset is simply angle off nadir and distance to ground
echo ("compute relative pixel offsets")
sampoff = create(pixels,nlines,1) # sample offset
lineoff = create(pixels,nlines,1) # line offset
height = altitude-topo # height above surface
mpixel = height*ptop/1000.0
xtmp = moment(mpixel)
meanpixel = xtmp[3,,] # mean pixel spacing
pitchoffpxl = height * tan(pitch/rad)/meanpixel
# compute yaw offset at edge of line
linyawoff = pedge * sin (yaw/rad)
sampyawoff = linyawoff * sin (yaw/rad) # sample offset is generally
# negligable
linoffvel = gvelocity / 12.0 / meanpixel
if (specprout == 1) {
sptitle=title1 + " Pixel Spacing (meters)"
write(translate(mpixel,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Sample Offset (pixels)"
write(translate(pitchoffpxl,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Line Offset (yaw, pixels)"
write(translate(linyawoff,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Sample Offset (yaw, pixels)"
write(translate(sampyawoff,y,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Line Offset (velocity, pixels)"
write(translate(linoffvel,y,z),filename=specpr_file,title=sptitle,type=specpr)
}
if(do_plot == 1) {
plot("set ylabel 'Meters'") # set vertical axis label
plot("set format y '%10.2f'")
plot("set title 'Pixel Spacing (meters)'")
plot("set output 'pixel.spacing.ps'")
plot(mpixel)
plot("set ylabel 'Pixels'") # set vertical axis label
plot("set format y '%10.2f'")
plot("set title 'Relative Pixel Offset due to Pitch'")
plot("set output 'pixel.pitch,offset.ps'")
plot(pitchoffpxl)
plot("set ylabel 'Pixels'") # set vertical axis label
plot("set format y '%10.2f'")
plot("set title 'Relative Line Offset at Edge of Scan due to Yaw (pixels)'")
plot("set output 'pixel.offset.edgeyaw.ps'")
plot(linyawoff)
plot("set ylabel 'Pixels'") # set vertical axis label
plot("set format y '%10.2f'")
plot("set title 'Relative Sample Offset at Edge of Scan due to Yaw (pixels)'")
plot("set output 'sample.offset.edgeyaw.ps'")
plot(sampyawoff)
plot("set ylabel 'Pixels'") # set vertical axis label
plot("set format y '%10.4f'")
plot("set title 'Relative Line Offset due to Ground Velocity (pixels)'")
plot("set output 'grnd.vel.rel.line.offset.ps'")
plot(linoffvel)
}
#############################################################################
# create cross track offest due to scan distortion onto flat surface
echo ("compute cross track offest due to scan angle")
xsample=float(create(pixels,1,1)) # sample offset array
crangle = xsample # array of cross tract angles
# compute now for use later
nadrad = roll[1,1,1] / rad # nadir pointing (aviris starts roll
# compensation here)
i = 1
while (i <= pixels) {
angle = float((i - pcenter)) * ptop/1000.0 - nadrad
crangle[i,1,1] = angle
itmp = i - pcenter
xsample[i,1,1] = (float(itmp)/cos(angle)) - float(itmp)
i = i + 1
}
if (specprout == 1) {
sptitle=title1 + " Sample Offset (scan angle, pxls)"
write(translate(xsample,x,z),filename=specpr_file,title=sptitle,type=specpr)
}
if(do_plot == 1) {
plot("set ylabel 'Pixels'") # set vertical axis label
plot("set format y '%10.2f'")
plot("set title 'Relative Sample Offset due to Cross-Track Scan Angls (pixels)'")
plot("set output 'sample.offset.scan.angle.ps'")
plot(translate(xsample,y,x))
}
######################################################################################
# compute cross track sampling from scan velocity and airplane height to uniform grid
echo ("compute cross track sampling from scan velocity and airplane height to uniform grid")
# notes: ptop is in milliradians
vcross = (altitude - topo) * tan(ptop/1000.0) # nadir pixel spacine along flight
#### later: xvsample = vcross * 0.0 # cross track uniform grid in pixel offsets
xpixel= create(pixels,1,1) * 0.0 # cross track array of pixels
for (i=1; i <= pixels; i = i+1) { # compute array of cross track pixels,
# negative to left of center line
xpixel[i,1,1] = i- pcenter
}
# because we computed the cosine effect above, we only need to compute
# idealized flat surface amount here
xpixeli = clone(float(xpixel),y=nlines) # xpixel image
pvsample = (xpixeli* vcross / grid) - xpixeli
####write (short(pvsample*10), filename="pvsample.v", type=vicar)
if (specprout == 1) {
sptitle=title1 + " Sample Offset (height, grid, pxls)"
write(translate(pvsample,x,z),filename=specpr_file,title=sptitle,type=specpr)
}
if(do_plot == 1) {
plot("set ylabel 'Pixels'") # set vertical axis label
plot("set format y '%10.2f'")
plot("set title 'Relative Sample Offset due to aircraft height, new grid size (pixels)'")
plot("set output 'sample.offset.scan.angle.ps'")
plot(translate(xsample,y,x))
}
#############################################################################
### compute bow in image because aircraft didn't travel in a straight line
echo ("compute bow in image because aircraft didn't travel in a straight line")
x1 = diflongs * sin((90.0 - vmeanhead)/rad)*cos(latitude/rad)
x2 = diflats * cos((90.0 - vmeanhead)/rad)
# The direction of bow can be found from the sign the difference in long, lat
# a + in longitude is to the east; a + in latitude is north.
# Which side of image is positive depends on direction of flight (meanhead).
xx = moment(vhead)
meanhead = xx[3] # mean heading in degrees, 0 = north
xx= vmeanhead - meanhead
echo ("Note: mean heading from GPS - mean velocity heading =" + xx + " degrees")
xdirectlong = -1.0 * cos(vmeanhead/rad)
xdirectlat = sin(meanhead/rad)
if (xdirectlong < 0.0) { # find general direction ignoring magnitude
xdirectlong = -1.0
} else {
xdirectlong = 1.0
}
if (xdirectlat < 0.0) { # find general direction ignoring magnitude
xdirectlat = -1.0
} else {
xdirectlat = 1.0
}
bow = x1 * xdirectlong + x2 * xdirectlat # in degrees
mbow = bow * meterdeg # in meters
pbow = mbow / grid # in pixels in the output grid
if (specprout == 1) {
sptitle=title1 + " Sample Bow Offset (m)"
write(translate(mbow,x,z),filename=specpr_file,title=sptitle,type=specpr)
sptitle=title1 + " Sample Bow Offset (pxl)"
write(translate(pbow,x,z),filename=specpr_file,title=sptitle,type=specpr)
}
if(do_plot == 1) {
plot("set ylabel 'Meters'") # set vertical axis label
plot("set format y '%10.1f'")
plot("set title 'Relative Sample Bow Offset from curved flight (meters)'")
plot("set output 'sample.offset.bow.m.ps'")
plot(mbow)
plot("set ylabel 'Pixels'") # set vertical axis label
plot("set format y '%10.1f'")
plot("set title 'Relative Sample Bow Offset from curved flight (pixels)'")
plot("set output 'sample.offset.bow.pxl.ps'")
plot(pbow)
}
### create image of sample offset
echo ("computing sample offset image")
ixsample = clone(short(xsample*10.0 + 0.5),y=nlines) # output offsets are 0.1 pixels/DN
ixsample = ixsample + short(pbow*10.0 + 0.5) + short(pvsample*10 + 0.5)
write (ixsample,filename="sample.offset.10.v",type=vicar)
#write (ixsample,filename="sample.offset.10.isis",type=isis)
#############################################################################
### now compute line offsets
ixline = ixsample * 0 # create ixline image same size as ixsample
il = 1
echo ("computing line offset image")
### compute cross track offset temp variable
xtmp = float(create(pixels,1,1))
is = 1
while ( is <= pixels ) {
xtmp[is,1,1] = float((is - pcenter))/float(pcenter) # normalized angle from edge of line
is = is + 1
}
xtmpi = float(clone(object=xtmp,y=nlines)) # make an image of normalized angles nlines long
# yaw and pitch offset for image + offset due to downtrack velicity variations
xline = linyawoff * xtmpi + pitchoffpxl + offsetdistance
ixline = short(xline * 10.0 + 0.5) # output offsets are 0.1 pixels/DN
write (ixline,filename="line.offset.10.v",type=vicar)
xline0 = linyawoff * xtmpi + pitchoffpxl # same as above but no resample to
# uniform grid, so we can see
# relative offsets ignoring velocity
# changes and gris resampling
ixline0 = short(xline0 * 10.0 + 0.5) # output offsets are 0.1 pixels/DN
write (ixline0,filename="line.offset.10.nodis.v",type=vicar)
echo ("done")
U.S. Geological Survey,
a bureau of the U.S. Department of the Interior
This page URL= http://speclab.cr.usgs.gov/aviris.geom.software/cmd.davinci.avirisgeom.8.html
This page is maintained by: Dr. Roger N. Clark rclark@speclab.cr.usgs.gov
Last modified December 15, 1998.