#!/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
#
# Revision 1.1 97/10/21 11:24:20 11:24:20 raymond (Ray Kokaly)
# Initial revision
#
# 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 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
#
# version 9: Additions to handle different AVIRIS 1998 navigation file.
# Changes to comments.
# - Raymond F. Kokaly 12/18/2000
#
###############################################################################
verbose=0 # BE SILENT
NUM_ARGS_REQD = 2
#
# READ COMMAND LINE OPTIONS
#
args=$ARGC # number of command line arguments
echo("command line arguments = " + args)
if (args < NUM_ARGS_REQD) {
printf("%s\n","ERROR: insufficient input")
Usage = $0 + " Flight_Year Navigation_File"
Usage = Usage + " [-eng_file Engineering_File]"
Usage = Usage + " [-title Plot_Title]"
Usage = Usage + " [-plot] [-specpr Specpr_Filename]"
Usage = Usage + " [-grid Grid_Size_In_Meters]"
Usage = Usage + " [-topo Height_In_Meters]"
printf("%s\n","Usage: " + Usage)
printf("%s\n","NOTE: for 1995 AVIRIS data an engineering data file must be specified.")
printf("%s\n","exit 1")
}
if (args < NUM_ARGS_REQD) { # BUG. THIS MUST BE SEPARATE OR
# LAST IF STATEMENT EXISTS BEFORE PRINTING.
#exit (1) # NO EXIT BECAUSE DAVINCI EXITS ALWAYS WHEN IT SEES IT
}
flt_year = atoi($1) # Flight Year (required)
file12 = $2 # Navigation File (required)
#
# ECHO COMMAND LINE ARGUMENTS TO USER
#
printf("%s\n","User Specified Command Line Arguments:")
printf("%s %i\n","Flight Year:",flt_year)
printf("%s %s\n","Navigation File:",file12)
#
# SET DEFAULTS
#
file11 = " "
title1 = "Navigation Data" # TITLE FOR PLOTS AND SPECPR FILES
do_plot = 0 # DEFAULT: NO PLOTS
do_eng = 0 # DEFAULT - DO NOT READ ENGINEERING DATA
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
#
# READ OPTIONAL COMMAND LINE ARGUMENTS
#
if (args > NUM_ARGS_REQD) {
for (i = NUM_ARGS_REQD+1; 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] == "-eng_file") {
do_eng = 1
if ( args >= i+1) { # SHOULD BE NAME OF ENGINEERING DATA FILE
i1 = i+1
txx=$ARGV[i1]
file11=txx
printf("%s\n","Engineering Data File is " + file11)
i = i1
} else {
printf("%s\n","ERROR: no filename after -eng_file")
printf("%s\n","exit 1")
exit_flag = 1
}
}
if($ARGV[i] == "-title") {
if ( args >= i+1) { # SHOULD BE PLOT TITLE
i1 = i+1
txx=$ARGV[i1]
title1=txx
printf("%s\n","Plot title is " + title1)
i = i1
} else {
printf("%s\n","ERROR: no plot title after -title")
printf("%s\n","exit 1")
exit_flag = 1
}
}
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 MEAN ELEVATION OF THE LINE
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) }
if (flt_year == 1995) {
if (do_eng == 0) {
printf("%s\n","ERROR: for 1995 AVIRIS data an engineering data file must be specified.\n")
printf("%s\n","Usage: " + Usage)
printf("%s\n","NOTE: for 1995 AVIRIS data an engineering data file must be specified.")
printf("%s\n","exit 1")
}
}
#
# FOR 1995 AVIRIS READ NAV DATA FROM NAV FILES
#
if (flt_year == 1995) {
echo ("reading navigation and engineering data files")
a=read(filename=file11)
b=read(filename=file12)
#
# 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
} else {
#
# FOR 1998 AVIRIS DATA SET UP THE NAV FILE PARAMETER LOCATIONS
#
echo ("reading navigation data file")
navdata = read_lines(filename=file12)
gpsstatus_pos = 2
gpsstatus_nchars = 1
day_pos = 4
day_nchars = 3
hr_pos = 8
hr_nchars = 2
min_pos = 11
min_nchars = 2
sec_pos = 14
sec_nchars = 2
lat_pos = 17
lat_nchars = 9
lon_pos = 27
lon_nchars = 10
thead_pos = 38
thead_nchars = 6
pitch_pos = 45
pitch_nchars = 8
roll_pos = 54
roll_nchars = 8
gspeed_pos = 63
gspeed_nchars = 6
trkangl_pos = 70
trkangl_nchars = 6
windspd_pos = 77
windspd_nchars = 4
winddir_pos = 82
winddir_nchars = 5
blongacc_pos = 88
blongacc_nchars = 6
blatacc_pos = 95
blatacc_nchars = 6
bnormacc_pos = 102
bnormacc_nchars = 6
trkangrate_pos = 109
trkangrate_nchars = 5
pitchrate_pos = 115
pitchrate_nchars = 5
rollrate_pos = 121
rollrate_nchars = 5
vertspd_pos = 127
vertspd_nchars = 6
gpsalt_pos = 134
gpsalt_nchars = 7
gpslat_pos = 142
gpslat_nchars = 9
gpslon_pos = 152
gpslon_nchars = 10
statpress_pos = 163
statpress_nchars = 8
totalpress_pos = 172
totalpress_nchars = 8
diffpress_pos = 181
diffpress_nchars = 6
totaltemp_pos = 188
totaltemp_nchars = 6
stattemp_pos = 195
stattemp_nchars = 6
baroalt_pos = 202
baroalt_nchars = 7
machno_pos = 210
machno_nchars = 5
airspd_pos = 216
airspd_nchars = 4
#
# EXTRACT AIRCRAFT ATTITUDE DATA FROM AVIRIS 1998 NAV DATA
#
navdata = read_lines(filename=file12)
roll = atof(navdata[roll_pos :roll_pos +roll_nchars -1,])
pitch = atof(navdata[pitch_pos :pitch_pos +pitch_nchars -1,])
airvelocity = atof(navdata[airspd_pos :airspd_pos +airspd_nchars -1,])
altitude = atof(navdata[gpsalt_pos :gpsalt_pos +gpsalt_nchars -1,])
vspeed = atof(navdata[vertspd_pos:vertspd_pos +vertspd_nchars-1,])
gvelocity = atof(navdata[gspeed_pos :gspeed_pos +gspeed_nchars -1,])
t = atof(navdata[thead_pos :thead_pos +thead_nchars -1,])
hours = atof(navdata[hr_pos :hr_pos +hr_nchars -1,])
minutes = atof(navdata[min_pos :min_pos +min_nchars -1,])
sec = atof(navdata[sec_pos :sec_pos +sec_nchars -1,])
latitude = double(atof(navdata[gpslat_pos:gpslat_pos+gpslat_nchars-1,]))
longitude = double(atof(navdata[gpslon_pos:gpslon_pos+gpslon_nchars-1,]))
dhours = hours + minutes/60.0 + sec/3600.0
td = dim(t)
nlines = td[2,,] # number of lines in data set
}
#
# SET AVIRIS PARAMETERS
#
ifov = 0.99 # IFOV in mrad (get from file 5)
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
#
# SET CONSTANTS
# the following earth radii from the book: Astrophysical Quantities, by Allen, 1976
# note: mean earth radius = 6371030 meters
#
pi = 3.14159265
rad = 360./(2.*pi) # degrees per radian
earthpolarradius= 6356779 # meters # NOTE: get exact value
eartheqradius = 6378164 # meters # NOTE: get exact value
#
# CALCULATE THE MEAN LATITUDE
#
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 DOESN'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)
}
#
# COMPUTE 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)
}
#
# CALCULATE PREDICTED LONGITUDE AND LATITUDE
#
echo ("predict longitude and latitude based on a straight line between ends of run")
alon_start = longitude[,1,1] # BEGINING LONGITUDE
alon_stop = longitude[,nlines,1] # ENDING LONGITUDE
alat_start = latitude[,1,1] # BEGINING LATITUDE
alat_stop = latitude[,nlines,1] # ENDING LATITUDE
dslong = alon_stop - alon_start
dslat = alat_stop - alat_start
slong = dslong / nlines # SLOPE OF LONGITUDE VS LINE
slat = dslat / nlines # SLOPE OF LATITUDE VS LINE
#
# CALCULATE MEAN HEADING FROM END POINTS OF LAT,LON
#
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) { vmeanhead = vmeanhead + 180. }
echo ("mean velocity heading from GPS = " + vmeanhead + " degrees")
vmeanheada = clone( vmeanhead, y=nlines)
offlong = alon_stop - slong * nlines
offlat = alat_stop - 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
}
diflon = longitude - predlong
diflat = latitude - predlat
echo ("smoothing difference longitude and difference latitude")
diflons = ysmooth(diflon,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(diflon,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(diflons,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 (diflon)
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 (diflons)
plot("set title 'Mean velocity heading for line from GPS (degrees)'")
plot("set output 'meanheadgps.ps'")
plot (vmeanheada)
}
#
# LATITUDE AND LONGITUDE ARE READ OUT ONCE EVERY SECOND, AND
# THE VALUES REPEATED FOR 12 SCAN LINES UNTIL THE NEXT VALUE IS
# ACQUIRED (IT SHOULD BE INTERPOLATED, SO WE MUST DO IT HERE)
#
### LONGITUDE
echo ("interpolating longitude to every scan line")
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
ifirst = 1
aline[1,1,1] = 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] = longitude[10,ifirst,1]
twonum[1,2,1] = longitude[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
echo ("interpolating latitude to every scan line")
alat = latitude[,,] # 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] = latitude[,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] = latitude[,ifirst,1]
twonum[1,2,1] = latitude[,i,1]
xnew=interp(twonum,twoxvl,aline[1,ifirst:ilast])
# 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)
}
#
# 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 LAT OR LON
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 velocity 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 velocity 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 TRACK 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
# NOTES: ptop is in milliradians
#
echo ("compute cross track sampling from scan velocity and airplane height to uniform grid")
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
# 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).
#
echo ("compute bow in image because aircraft didn't travel in a straight line")
x1 = diflons * sin((90.0 - vmeanhead)/rad)*cos(latitude/rad)
x2 = diflats * cos((90.0 - vmeanhead)/rad)
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 VELOCITY 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 GRID 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.9.html
This page is maintained by: Dr. Roger N. Clark rclark@speclab.cr.usgs.gov
Last modified Jan. 10, 2001.