cmd.davinci.avirisgeom.9

http://speclab.cr.usgs.gov


#!/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.