cmd.davinci.avirisgeom.8

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
# $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.