# ------------------------------------------------------------
# SRFI.sml
# ------------------------------------------------------------
# SET WARNING LEVEL: Refer to B5 & B6.
$warnings 3
# ------------------------------------------------------------
# DEFINE PROCEDURE writeTitle: Refer to B12 & B13.
# PURPOSE: WRITES TITLE & AUTHOR INFO TO CONSOLE WINDOW.
proc writeTitle() begin
printf("SRFI.sml:\n");
printf(" SML VERSION: November 15, 2005\n");
printf(" TNTmips VERSION: Released 7.1\n");
printf(" PURPOSE: CONVERTS MS DNs ");
printf("TO SRFI VALUES.\n");
printf(" MAY PRODUCE: PVI & PBI RASTERS\n");
printf(" DETAILS: FAQs_by_Jack A & B\n");
printf(" AUTHOR: Dr. Jack F. Paris\n");
printf(" CONTACT INFO: [email protected] ");
printf(" 303-775-1195\n");
printf(" ALLOWED USE: ONLY NON-COMMERCIAL\n\n");
end
# ------------------------------------------------------------
# DEFINE VARIABLES FOR GENERAL CHARACTER STRINGS: Refer to B7.
string t$,p$,p1$,p2$,p3$,p4$,p5$,p6$,p7$,p8$,p9$;
string p10$,p11$,p12$,p13$,p14$,p15$,p16$,p17$,p18$,p19$;
string p20$,p21$;
# ------------------------------------------------------------
# CLEAR CONSOLE WINDOW & ASK USER TO RE-POSITION IT:
# Refer to B8 & B9.
clear();
p1$ = "CONSOLE-WINDOW ADJUSTMENT:\n";
p2$ = "* REPOSITION the CONSOLE WINDOW.\n";
p3$ = "* Then, CLICK OK.";
p$ = p1$ + p2$ + p3$;
PopupMessage(p$);
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO DATETIME & CONTRAST CLASSES:
# Refer to B10 - B11.
class DATETIME cdate$;
class CONTRAST smlContrast;
string clut1$,clut2$,clut3$;
# ------------------------------------------------------------
# ASSIGN VALUES TO VARIABLES RELATED TO smlContrast:
clut1$ = "exp"; clut2$ = "Exponential";
clut3$ = "Exponential contrast table";
smlContrast.Power = 0.8;
smlContrast.InputLowerLimit = 1;
smlContrast.OutputLowerLimit = 1;
smlContrast.OutputUpperLimit = 255;
# ------------------------------------------------------------
# WRITE TITLE & AUTHOR INFORMATION: Refer to B13.
writeTitle();
# ------------------------------------------------------------
# DECLARE USER-INPUT VARIABLES: Refer to B14 & B15.
string site$,gc$;
numeric cdate,pdate,dd,d1,d2;
numeric sunelevang,imager,vnir;
numeric atcor,delcf,msfac,icRL;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO DATA TYPES OF RASTERS:
# Refer to B14 & B15.
string dntype$,srfitype$,ptype$;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO USER-INPUT VARIABLES:
# Refer to B14 & B15.
numeric he,doy,sc,dnlow,dnhigh;
string imager$,atcor$,gi$;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO SRFapc-to-SRFsfc CONVERSION:
# Refer to B14 & B15.
numeric cCB,cBL,cGL,cYL,cRL,cRE,cNA,cNB;
numeric cMA,cMB,cMC,cMD,cME,cMF,cMG,cRLm1,pc;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO BOOLEAN ENABLERS:
# Refer to B14 & B15.
numeric pF,pCB,pBL,pYL,pRE,pNA,pNB;
numeric pMA,pMB,pMC,pMD,pME,pMF,pMG,nc;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO WAVELENGTH-BASED MODELS:
# Refer to B14 & B15.
numeric p,p2,wLenCB,wLenBL,wLenGL,wLenYL,wLenRL,wLenRE;
numeric wLenNA,wLenNB,wLenMA,wLenMB,wLenMC,wLenMD,wLenME;
numeric wLenMF,wLenMG;
numeric wLenV,logwLenRLdwLenV,wtBL;
numeric srfpathV,logsrfpathV;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO DN-to-SR CONVERSION:
# Refer to B14 & B15.
numeric skCB,skBL,skGL,skYL,skRL,skRE,skNA,skNB;
numeric skMA,skMB,skMC,skMD,skME,skMF,skMG;
numeric dnbCB,dnbBL,dnbGL,dnbYL,dnbRL,dnbRE,dnbNA,dnbNB;
numeric dnbMA,dnbMB,dnbMC,dnbMD,dnbME,dnbMF,dnbMG;
# ------------------------------------------------------------
# DECLARE DIRECT-SOLAR-SPECTRAL-IRRADIANCE VARIABLES:
# Refer to B14 & B15.
numeric dssiCB,dssiBL,dssiGL,dssiYL,dssiRL,dssiRE;
numeric dssiNA,dssiNB,dssiMA,dssiMB,dssiMC,dssiMD,dssiME;
numeric dssiMF,dssiMG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO PVI & PBI: Refer to B14 & B15.
numeric pfac,angrad,sinang,cosang;
numeric pbioff,pvioff,srfiNRint,bslineslope,maxpvipbi;
numeric srfiNR,tsrfiNR,pbif,pbi,pvif,pvi;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO ARRAY PARAMETERS:
# Refer to B14 & B15.
numeric maxi,maxip1;
# ------------------------------------------------------------
# DECLARE LIST OF POSSIBLE INPUT RASTERS: Refer to B14 & B15.
raster CB,BL,GL,YL,RL,RE,NA,NB,MA,MB,MC,MD,ME,MF,MG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO LOOPS: Refer to B14 & B15.
numeric i,nlins,ncols;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO INPUT RASTERS:
# Refer to B14 & B15.
numeric dnCB,dnBL,dnGL,dnYL,dnRL,dnRE,dnNA,dnNB;
numeric dnMA,dnMB,dnMC,dnMD,dnME,dnMF,dnMG;
numeric dnminCB,dnminBL,dnminGL,dnminYL,dnminRL,dnminRE;
numeric dnminNA,dnminNB,dnminMA,dnminMB,dnminMC,dnminMD;
numeric dnminME,dnminMF,dnminMG;
numeric dnmaxCB,dnmaxBL,dnmaxGL,dnmaxYL,dnmaxRL,dnmaxRE;
numeric dnmaxNA,dnmaxNB,dnmaxMA,dnmaxMB,dnmaxMC,dnmaxMD;
numeric dnmaxME,dnmaxMF,dnmaxMG;
# ------------------------------------------------------------
# DECLARE VARIABLES FOR HISTOGRAM-EDGE FINDING ALGORITHM:
# Refer to B14 & B15.
numeric dnheCB,dnheBL,dnheGL,dnheYL,dnheRL,dnheRE;
numeric dnheNA,dnheNB,dnheMA,dnheMB,dnheMC,dnheMD;
numeric dnheME,dnheMF,dnheMG;
numeric logsrfpathRL,difoflogs;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO SItoa & SRFtoa:
# Refer to B14 & B15.
numeric pi100,esd,sitoafac;
numeric sitoaCB,sitoaBL,sitoaGL,sitoaYL,sitoaRL,sitoaRE;
numeric sitoaNA,sitoaNB,sitoaMA,sitoaMB,sitoaMC,sitoaMD;
numeric sitoaME,sitoaMF,sitoaMG;
numeric aCB,aBL,aGL,aYL,aRL,aRE,aNA,aNB;
numeric aMA,aMB,aMC,aMD,aME,aMF,aMG;
numeric mCB,mBL,mGL,mYL,mRL,mRE,mNA,mNB;
numeric mMA,mMB,mMC,mMD,mME,mMF,mMG;
numeric srftoaCB,srftoaBL,srftoaGL,srftoaYL,srftoaRL;
numeric srftoaRE,srftoaNA,srftoaNB,srftoaMA,srftoaMB;
numeric srftoaMC,srftoaMD,srftoaME,srftoaMF,srftoaMG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO DNpath: Refer to B14 & B15.
numeric dnpathCB1,dnpathBL1,dnpathGL1,dnpathYL1,dnpathRE1;
numeric dnpathNA1,dnpathNB1,dnpathMA1,dnpathMB1,dnpathMC1;
numeric dnpathMD1,dnpathME1,dnpathMF1,dnpathMG1;
numeric dnpathCB,dnpathBL,dnpathGL,dnpathYL,dnpathRL;
numeric dnpathRE,dnpathNA,dnpathNB,dnpathMA,dnpathMB;
numeric dnpathMC,dnpathMD,dnpathME,dnpathMF,dnpathMG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO SRFpath: Refer to B14 & B15.
numeric srfpathCB1,srfpathBL1,srfpathGL1,srfpathYL1;
numeric srfpathRE1,srfpathNA1,srfpathNB1,srfpathMA1;
numeric srfpathMB1,srfpathMC1,srfpathMD1,srfpathME1;
numeric srfpathMF1,srfpathMG1;
numeric srfpathCB2,srfpathBL2,srfpathGL2,srfpathYL2;
numeric srfpathRE2,srfpathNA2,srfpathNB2,srfpathMA2;
numeric srfpathMB2,srfpathMC2,srfpathMD2,srfpathME2;
numeric srfpathMF2,srfpathMG2;
numeric srfpathCB,srfpathBL,srfpathGL,srfpathYL,srfpathRL;
numeric srfpathRE,srfpathNA,srfpathNB,srfpathMA,srfpathMB;
numeric srfpathMC,srfpathMD,srfpathME,srfpathMF,srfpathMG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO SRFapc: Refer to B14 & B15.
numeric srfapcCB,srfapcBL,srfapcGL,srfapcYL,srfapcRL;
numeric srfapcRE,srfapcNA,srfapcNB,srfapcMA,srfapcMB;
numeric srfapcMC,srfapcMD,srfapcME,srfapcMF,srfapcMG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO SRFsfc & SRFI:
# Refer to B14 & B15.
numeric srfsfcCB,srfsfcBL,srfsfcGL,srfsfcYL,srfsfcRL;
numeric srfsfcRE,srfsfcNA,srfsfcNB,srfsfcMA,srfsfcMB;
numeric srfsfcMC,srfsfcMD,srfsfcME,srfsfcMF,srfsfcMG;
numeric srfiCB,srfiBL,srfiGL,srfiYL,srfiRL,srfiRE,srfiNA;
numeric srfiNB,srfiMA,srfiMB,srfiMC,srfiMD,srfiME,srfiMF;
numeric srfiMG;
# ------------------------------------------------------------
# DECLARE LIST OF POSSIBLE OUTPUT RASTERS: Refer to B14 & B15.
raster SRFICB,SRFIBL,SRFIGL,SRFIYL,SRFIRL,SRFIRE,SRFINA;
raster SRFINB,SRFIMA,SRFIMB,SRFIMC,SRFIMD,SRFIME,SRFIMF;
raster SRFIMG,PBI,PVI;
# ------------------------------------------------------------
# GET NAME OF SITE FROM THE USER: Refer to B16.
t$ = "SITE NAME";
p1$ = "SITE-NAME:\n";
p2$ = "* ENTER SITE NAME.\n";
p3$ = "* Then, CLICK OK.\n";
p4$ = "NAME ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$;
site$ = PopupString(p$,"NAME",t$);
printf(" SITE NAME: %s\n",site$);
# ------------------------------------------------------------
# GET COLLECTION DATE FROM THE USER: Refer to B17.
p1$ = "COLLECTION-DATE:\n";
p2$ = " FORMAT: YYYYMMDD\n";
p3$ = "* ACCEPT the Default DATE,\n";
p4$ = "* Or, ENTER the Correct DATE.\n";
p5$ = "* Then, Click OK.\n\n";
p6$ = "DATE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$;
dd = 20010509; d1 = 19710723; d2 = 29991231;
cdate = PopupNum(p$,dd,d1,d2,0);
dd = cdate; d1 = cdate;
cdate$.SetDateYYYYMMDD(cdate);
doy = cdate$.GetDayOfYear();
printf(" COLLECTION DATE: %8d\n",cdate);
printf(" DAY OF THE YEAR: %d\n",doy);
# ------------------------------------------------------------
# GET SUN ELEVATION ANGLE FROM THE USER: Refer to B18.
p1$ = "SUN-ELEVATION-ANGLE:\n";
p2$ = " UNITS: deg above the Horizon\n";
p3$ = " FORMAT: NN.NN\n";
p4$ = " RANGE: 1.00 to 90.00 deg\n";
p5$ = "* ACCEPT the Default ANGLE,\n";
p6$ = "* Or, ENTER the Correct ANGLE.\n";
p7$ = "* Then, CLICK OK.\n\n";
p8$ = "ANGLE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
sunelevang = PopupNum(p$,66.99,1,90,2);
printf(" SUN ELEV. ANGLE: %5.2f deg\n",sunelevang);
# ------------------------------------------------------------
# GET IMAGER NUMBER FROM THE USER: Refer to B19.
p1$ = "IMAGER-NUMBER:\n";
p2$ = " IMAGER\n";
p3$ = " NUMBER: NAME OF IMAGER STATUS\n";
p4$ = " 1: QuickBird MS Available\n";
p5$ = " 2: Ikonos MS Available\n";
p6$ = " 3: OrbView-3 MS Not\n";
p7$ = " 4: Landsat-7 ETM+ Available\n";
p8$ = " 5: Landsat-5 TM Available\n";
p9$ = " 6: Landsat-5 MSS Not\n";
p10$ = " 7: Landsat-4 TM Available\n";
p11$ = " 8: Landsat-4 MSS Not\n";
p12$ = " 9: Landsat-3 MSS Not\n";
p13$ = " 10: Landsat-2 MSS Not\n";
p14$ = " 11: Landsat-1 MSS Not\n";
p15$ = " 12: Terra ASTER Available\n";
p16$ = " 13: Terra MODIS Not\n";
p17$ = " 14: Aqua MODIS Not\n";
p18$ = "* ACCEPT the Default IMAGER NUMBER,\n";
p19$ = "* Or, SELECT the Correct IMAGER NUMBER.\n";
p20$ = "* Then, CLICK OK.\n\n";
p21$ = "NUMBER ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$ + p9$;
p$ = p$ + p10$ + p11$ + p12$ + p13$ + p14$ + p15$ + p16$;
p$ = p$ + p17$ + p18$ + p19$ + p20$ + p21$;
imager = PopupNum(p$,12,14,1,0);
# ------------------------------------------------------------
# GET ATMOSPHERIC CORRECTION (atcor) LEVEL FROM THE USER:
# Refer to B20.
p1$ = "ATMOSPHERIC-CORRECTION LEVEL:\n";
p2$ = " LEVEL: TYPE OF ATMOSPHERIC CORRECTION\n";
p3$ = " 1: No Atmospheric Correction\n";
p4$ = " 2: Atmospheric Reflectance Only\n";
p5$ = " 3: All Atmospheric Effects\n";
p6$ = "* ACCEPT the Default LEVEL,\n";
p7$ = "* Or, SELECT a Differernt LEVEL.\n";
p8$ = "* Then, CLICK OK.\n\n";
p9$ = "LEVEL ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$ + p9$;
atcor = PopupNum(p$,3,1,3,0);
# ------------------------------------------------------------
# WRITE CORRECTION LEVEL TO THE CONSOLE WINDOW:
# Refer to B20.
printf("CORRECTION LEVEL: ");
p4$ = "No Atmospheric Corrections";
p5$ = "Corrected for Atmospheric Reflectance";
p6$ = "Corrected for All Atmospheric Effects";
if (atcor==1) then printf("%s\n",p4$);
if (atcor==2) then printf("%s\n",p5$);
if (atcor==3) then printf("%s\n",p6$);
# ------------------------------------------------------------
# IF atcor > 1, GET delcf PARAMETER FROM USER:
# Refer to B22.
if (atcor > 1) then begin
p1$ = "HISTOGRAM-EDGE FINDING PARAMETER:\n";
p2$ = " FORMAT: N.NN (Units: Percentage Points)\n";
p3$ = " RANGE: 0.01 to 1.00 Percentage Points\n";
p4$ = "* Either, ACCEPT the Default VALUE,\n";
p5$ = "* Or, ENTER a Better VALUE.\n";
p6$ = "* Then, CLICK OK.\n\n";
p7$ = "VALUE ENTERED:";
p$ = p1$ + p2$ +p3$ + p4$ + p5$ + p6$ + p7$;
delcf = PopupNum(p$,0.05,0.01,1.00,2);
printf(" delcf: %4.2f Percentage",delcf);
printf(" Points\n");
end
he = 0.01 * delcf;
# ------------------------------------------------------------
# IF atcor = 3, GET the msfac VALUE: Refer B23.
if (atcor == 3) then begin
p1$ = "msfac-VALUE ENTRY:\n";
p2$ = " msfac: Scaling Factor for All MS Bands.\n";
p3$ = " FORMAT: N.NNNN\n";
p4$ = " RANGE: 0.1000 to 2.0000\n";
p5$ = "* ACCEPT the Default VALUE,\n";
p6$ = "* Or, ENTER a Better VALUE.\n";
p7$ = " NOTE: Don't Change msfac Unless Justified.\n";
p8$ = "* Then, CLICK OK.\n\n";
p9$ = "VALUE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
p$ = p$ + p9$;
msfac = PopupNum(p$,1,0.1,2,4);
printf(" msfac: %6.4f\n",msfac);
# ------------------------------------------------------------
# IF atcor = 3, GET icRL VALUE FROM THE USER:
# Refer to B24.
p1$ = "icRL-VALUE ENTRY:\n";
p2$ = " icRL: Input c-Factor for the RL Band.\n";
p3$ = " FORMAT: N.NNNN\n";
p4$ = " RANGE: 1.0000 to 3.0000\n";
p5$ = "* ACCEPT the Default Value,\n";
p6$ = "* Or, ENTER a Better Value.\n";
p7$ = "* Then, CLICK OK.\n\n";
p8$ = "VALUE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
icRL = PopupNum(p$,1.34,1,2,4);
printf(" icRL: %6.4f\n",icRL);
cRL = icRL;
end
# ------------------------------------------------------------
# CALCULATE cRLm1 & ASSIGN THE pc VALUE:
# Refer to B26.
cRLm1 = cRL - 1; pc = 2.27;
# ------------------------------------------------------------
# GENERATE IMAGER-SPECIFIC PARAMETERS: Refer to B25.
if (imager == 1) then begin
imager$ = "QuickBird MS";
pBL=1; pNA=1; maxi=2100; dnlow=1; dnhigh=2047;
wLenBL=0.482; wLenGL=0.548; wLenRL=0.654; wLenNA=0.809;
skBL=0.2359; skGL=0.1453; skRL=0.1785; skNA=0.1353;
dnbBL=0; dnbGL=0; dnbRL=0; dnbNA=0;
dssiBL=1925; dssiGL=1843; dssiRL=1575; dssiNA=1250;
end
if (imager == 2) then begin
imager$ = "Ikonos MS";
pBL=1; pNA=1; maxi=2100; dnlow=1; dnhigh=2047;
wLenBL=0.480; wLenGL=0.551; wLenRL=0.665; wLenNA=0.805;
skBL=0.2216; skGL=0.1739; skRL=0.1809; skNA=0.1405;
if (cdate > 20010222) then begin
skBL=0.1935; skGL=0.1847; skRL=0.1536; skNA=0.1148;
end
dnbBL=0; dnbGL=0; dnbRL=0; dnbNA=0;
dssiBL=1939; dssiGL=1843; dssiRL=1575; dssiNA=1250;
end
if (imager == 3) then begin
p1$ = "WARNING!\n";
p2$ = "SRFI.sml Cannot Process MS Data from\n";
p3$ = "OrbView 3.\n";
p4$ = "Please Stop this SML by\n";
p5$ = "CLICKing CANCEL.";
p$ = p1$ + p2$ + p3$ + p4$ + p5$;
PopupMessage(p$);
imager$ = "OrbView-3 MS";
pBL=1; pNA=1; maxi=2100; dnlow=1; dnhigh=2047;
wLenBL=0.482; wLenGL=0.548; wLenRL=0.654; wLenNA=0.809;
# skBL=0.2359; skGL=0.1453; skRL=0.1785; skNA=0.1353;
dnbBL=0; dnbGL=0; dnbRL=0; dnbNA=0;
dssiBL=1925; dssiGL=1843; dssiRL=1575; dssiNA=1250;
end
if (imager == 4) then begin
imager$ = "Landsat-7 ETM+";
pBL=1; pNA=1; pMB=1; pMC=1; maxi=255;
dnlow=0; dnhigh=255;
wLenBL=0.482; wLenGL=0.565; wLenRL=0.660; wLenNA=0.825;
wLenMB=1.650; wLenMC=2.220;
# GET PROCESSING DATE FROM THE USER:
p1$ = "PROCESSING-DATE ENTRY\n";
p2$ = " FORMAT: YYYYMMDD\n";
p3$ = "* Either ACCEPT the Default DATE,\n";
p4$ = "* Or, ENTER the Correct DATE.\n";
p5$ = "* Then, Click OK.\n\n";
p6$ = "PROCESSING-DATE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$;
pdate = PopupNum(p$,dd,d1,d2,0);
# GET PRODUCT-SOURCE CODE FROM USER:
p1$ = "LANDSAT 7 ETM+ PRODUCT SOURCE CODE:\n";
p2$ = " SOURCE CODE: SOURCE NAME\n";
p3$ = " 1: LPGS (EOS Data Gateway)\n";
p4$ = " 2: NLAPS (EarthExplorer)\n";
p5$ = " Either, ACCEPT the Default CODE,\n";
p6$ = " Or, ENTER the Correct CODE.\n";
p7$ = " Then, CLICK OK.\n";
p8$ = "CODE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
sc = PopupNum(p$,2,1,2,0);
# GET GAIN-CODE STRING FROM USER:
p1$ = "LANDSAT 7 ETM+ GAIN-CODE:\n";
p2$ = " GAIN-CODE BAND SEQUENCE:\n";
p3$ = " BL GAIN: H or L : Default=H\n";
p4$ = " GL GAIN: H or L : Default=H\n";
p5$ = " RL GAIN: H or L : Default=H\n";
p6$ = " NA GAIN: H or L : Default=L\n";
p7$ = " MB GAIN: H or L : Default=H\n";
p8$ = " MC GAIN: H or L : Default=H\n";
p9$ = " Default GAIN-CODE = HHHLHH\n";
p10$ = " ACCEPT the Default GAIN-CODE,\n";
p11$ = " Or, ENTER the Correct GAIN-CODE.\n";
p12$ = " Then, CLICK OK.\n";
p13$ = "GAIN-CODE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
p$ = p$ + p9$ + p10$ + p11$ + p12$ + p13$;
gc$ = "HHHLHH";
t$ = "GAIN-CODE QUERY";
gc$ = PopupString(p$,gc$,t$);
# Low-Gain Values for Data Processed BEFORE 20000701 by
# NLAPS (EarthExplorer), i.e., qcalmin = 0:
skBL=1.19098; skGL=1.21333; skRL=0.94118; skNA=0.93922;
skMB=0.19098; skMC=0.06647;
dnbBL=5.21; dnbGL=4.95; dnbRL=4.78; dnbNA=4.79;
dnbMB=5.24; dnbMC=5.27;
# High-Gain Values for Data Processed BEFORE 20000701 by
# NLAPS (EarthExplorer), i.e., qcalmin = 0:
gi$ = mid$(gc$,1,1);
if (gi$ == "H") then begin
skBL=0.78627; dnbBL=7.89;
end
gi$ = mid$(gc$,2,1);
if (gi$ == "H") then begin
skGL=0.81725; dnbGL=7.34;
end
gi$ = mid$(gc$,3,1);
if (gi$ == "H") then begin
skRL=0.63961; dnbRL=7.04;
end
gi$ = mid$(gc$,4,1);
if (gi$ == "H") then begin
skNA=0.63529; dnbNA=7.08;
end
gi$ = mid$(gc$,5,1);
if (gi$ == "H") then begin
skMB=0.12847; dnbMB=7.78;
end
gi$ = mid$(gc$,6,1);
if (gi$ == "H") then begin
skMC=0.04424; dnbMC=7.91;
end
if (pdate > 20000701) then begin
# Low-Gain Values for Data Processed AFTER 20000701 by
# NLAPS (EarthExplorer), i.e., qcalmin = 0:
skBL=1.17608; skGL=1.20510; skRL=0.93882; skNA=0.96549;
skMB=0.19047; skMC=0.06624;
dnbBL=5.27; dnbGL=5.31; dnbRL=5.33; dnbNA=5.28;
dnbMB=5.25; dnbMC=5.28;
# High-Gain Values for Data Processed AFTER 20000701 by
# NLAPS (EarthExplorer), i.e., qcalmin = 0:
gi$ = mid$(gc$,1,1);
if (gi$ == "H") then begin
skBL=0.77569; dnbBL=7.99;
end
gi$ = mid$(gc$,2,1);
if (gi$ == "H") then begin
skGL=0.79569; dnbGL=8.04;
end
gi$ = mid$(gc$,3,1);
if (gi$ == "H") then begin
skRL=0.61922; dnbRL=8.07;
end
gi$ = mid$(gc$,4,1);
if (gi$ == "H") then begin
skNA=0.63725; dnbNA=8.00;
end
gi$ = mid$(gc$,5,1);
if (gi$ == "H") then begin
skMB=0.12573; dnbMB=7.95;
end
gi$ = mid$(gc$,6,1);
if (gi$ == "H") then begin
skMC=0.04373; dnbMC=8.00;
end
end
# ---------------------------------------------------------
# Corrections Related to Changing qcalmin from 0 to 1:
# This Change is Associated with Either of Two Cases:
# * sc == 1 (LPGS: All Dates)
# * pdate > 20040701 (LPGS and NLAPS)
if (sc == 1 or pdate > 20040404) then begin
dnlow = 1;
skBL = skBL*255/254; dnbBL = 1 + dnbBL*254/255;
skGL = skGL*255/254; dnbGL = 1 + dnbGL*254/255;
skRL = skRL*255/254; dnbRL = 1 + dnbRL*254/255;
skNA = skNA*255/254; dnbNA = 1 + dnbNA*254/255;
skMB = skMB*255/254; dnbMB = 1 + dnbMB*254/255;
skMC = skMC*255/254; dnbMC = 1 + dnbMC*254/255;
end
dssiBL=1970; dssiGL=1843; dssiRL=1555; dssiNA=1047;
dssiMB=215.0; dssiMC=80.67;
end
if (imager == 5) then begin
imager$ = "Landsat-5 TM";
pBL=1; pNA=1; pMB=1; pMC=1; maxi=255;
dnlow=0; dnhigh=255;
wLenBL=0.482; wLenGL=0.565; wLenRL=0.660; wLenNA=0.825;
wLenMB=1.650; wLenMC=2.220;
skBL=0.60243; skGL=1.17510; skRL=0.80576; skNA=0.81455;
skMB=0.10808; skMC=0.05698;
dnbBL=2.52; dnbGL=2.42; dnbRL=1.45; dnbNA=1.85;
dnbMB=3.42;dnbMC=2.63;
dssiBL=1957; dssiGL=1829; dssiRL=1554; dssiNA=1036;
dssiMB=215.0; dssiMC=80.67;
# GET PROCESSING DATE FROM THE USER:
p1$ = "PROCESSING-DATE:\n";
p2$ = " FORMAT: YYYYMMDD\n";
p3$ = "* ACCEPT the Default DATE,\n";
p4$ = "* Or, ENTER the Correct DATE.\n";
p5$ = "* Then, Click OK.\n\n";
p6$ = "DATE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$;
pdate = PopupNum(p$,dd,d1,d2,0);
if (pdate > 20030504) then begin
skBL=0.76282; skGL=1.44251; skRL=1.03988;
skNA=0.87259; skMB=0.11988; skMC=0.06529;
dnbBL=1.99; dnbGL=1.97; dnbRL=1.13; dnbNA=1.73;
dnbMB=3.09; dnbMC=2.30;
end
end
if (imager == 6) then begin
p1$ = "WARNING!\n";
p2$ = "SRFI.sml Cannot Process MS Data from\n";
p3$ = "Landsat-5 MSS.\n";
p4$ = "Please Stop this SML by\n";
p5$ = "CLICKing CANCEL.";
p$ = p1$ + p2$ + p3$ + p4$ + p5$;
PopupMessage(p$);
imager$ = "Landsat-5 MSS";
pRE=1; pNB=1; maxi=200; dnlow=0; dnhigh=127;
wLenGL=0.550; wLenRL=0.650; wLenRE=0.750;
wLenNB=0.950;
# skGL=1.6000; skRL=1.6000; skRE=1.6000; skNB=1.6000;
# dnbGL=0; dnbRL=0; dnbRE=0; dnbNA=0;
# dssiGL=1829; dssiRL=1557; dssiRE=1300; dssiNB=950.0;
end
if (imager == 7) then begin
p1$ = "WARNING!\n";
p2$ = "SRFI.sml Cannot Process MS Data from\n";
p3$ = "Landsat 4 TM.\n";
p4$ = "Please Stop this SML by\n";
p5$ = "CLICKing CANCEL.";
p$ = p1$ + p2$ + p3$ + p4$ + p5$;
PopupMessage(p$);
imager$ = "Landsat-4 TM";
pBL=1; pNA=1; pMB=1; pMC=1; maxi=255;
dnlow=0; dnhigh=255;
wLenBL=0.482; wLenGL=0.565; wLenRL=0.660; wLenNA=0.825;
wLenMB=1.650; wLenMC=2.220;
# skBL=1.6000; skGL=1.6000; skRL=1.6000; skNA=1.6000;
# skMB=1.6000; skMC=1.6000;
# dnbBL=0; dnbGL=0; dnbRL=0; dnbNA=0;
# dnbMB=0; dnbMC=0;
dssiBL=1957; dssiGL=1825; dssiRL=1557; dssiNA=1033;
dssiMB=215.0; dssiMC=80.67;
end
if (imager == 8) then begin
p1$ = "WARNING!\n";
p2$ = "SRFI.sml Cannot Process MS Data from\n";
p3$ = "Landsat 4 MSS.\n";
p4$ = "Please Stop this SML by\n";
p5$ = "CLICKing CANCEL.";
p$ = p1$ + p2$ + p3$ + p4$ + p5$;
PopupMessage(p$);
imager$ = "Landast-4 MSS"; pRE=1; pNB=1; maxi=200;
dnlow=0; dnhigh=127;
wLenGL=0.550; wLenRL=0.650; wLenRE=0.750;
wLenNB=0.950;
# skGL=1.6000; skRL=1.6000; skRE=1.6000;
# skNB=1.6000;
# dnbGL=0; dnbRL=0; dnbRE=0; dnbNA=0;
# dssiGL=1829; dssiRL=1557; dssiRE=1300;
# dssiNB=950.0;
end
if (imager == 9) then begin
p1$ = "WARNING!\n";
p2$ = "SRFI.sml Cannot Process MS Data from\n";
p3$ = "Landsat 3 MSS.\n";
p4$ = "Please Stop this SML by\n";
p5$ = "CLICKing CANCEL.";
p$ = p1$ + p2$ + p3$ + p4$ + p5$;
PopupMessage(p$);
imager$ = "Landsat-3 MSS"; pRE=1; pNB=1; maxi=255;
dnlow=0; dnhigh=127;
wLenGL=0.550; wLenRL=0.650; wLenRE=0.750;
wLenNB=0.950;
# skGL=1.6000; skRL=1.6000; skRE=1.6000;
# skNB=1.6000;
# dnbGL=0; dnbRL=0; dnbRE=0; dnbNA=0;
# dssiGL=1829; dssiRL=1557; dssiRE=1300;
# dssiNB=950.0;
end
if (imager == 10) then begin
p1$ = "WARNING!\n";
p2$ = "SRFI.sml Cannot Process MS Data from\n";
p3$ = "Landsat 2 MSS.\n";
p4$ = "Please Stop this SML by\n";
p5$ = "CLICKing CANCEL.";
p$ = p1$ + p2$ + p3$ + p4$ + p5$;
PopupMessage(p$);
imager$ = "Landsat-2 MSS"; pRE=1; pNB=1; maxi=100;
dnlow=0; dnhigh=63;
wLenGL=0.550; wLenRL=0.650; wLenRE=0.750;
wLenNB=0.950;
# skGL=1.6000; skRL=1.6000; skRE=1.6000;
# skNB=1.6000;
# dnbGL=0; dnbRL=0; dnbRE=0; dnbNA=0;
# dssiGL=1829; dssiRL=1557; dssiRE=1300;
# dssiNB=950.0;
end
if (imager == 11) then begin
p1$ = "WARNING!\n";
p2$ = "SRFI.sml Cannot Process MS Data from\n";
p3$ = "Landsat 1 MSS.\n";
p4$ = "Please Stop this SML by\n";
p5$ = "CLICKing CANCEL.";
p$ = p1$ + p2$ + p3$ + p4$ + p5$;
PopupMessage(p$);
imager$ = "Landsat-1 MSS"; pRE=1; pNB=1; maxi=100;
dnlow=0; dnhigh=63;
wLenGL=0.550; wLenRL=0.650; wLenRE=0.750;
wLenNB=0.950;
# skGL=1.6000; skRL=1.6000; skRE=1.6000;
# skNB=1.6000;
# dnbGL=0; dnbRL=0; dnbRE=0; dnbNA=0;
# dssiGL=1829; dssiRL=1557; dssiRE=1300;
# dssiNB=950.0;
end
# ---------------------------------------------------------
if (imager == 12) then begin
imager$ = "Terra ASTER"; pNA=1; pMB=1; pMC=1; pMD=1;
pME=1; pMF=1; pMG=1; maxi=255;
dnlow=1; dnhigh=255;
wLenGL=0.556; wLenRL=0.659; wLenNA=0.807;
wLenMB=1.657; wLenMC=2.169; wLenMD=2.209; wLenME=2.263;
wLenMF=2.334; wLenMG=2.400;
p1$ = "VNIR BANDS ONLY versus ALL BANDS OPTION:\n";
p2$ = " OPTION 1: PROCESS VNIR BANDS ONLY\n";
p3$ = " OPTION 2: PROCESS ALL BANDS\n";
p4$ = " NOTE: INPUT BANDS MUST MATCH.\n";
p5$ = "OPTION ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$;
vnir = PopupNum(p$,2,1,2,0);
# DEFAULT (NORMAL-GAIN) sk VALUES:
skGL=1.6880; skRL=1.4150; skNA=0.8620; skMB=0.2174;
skMC=0.0696; skMD=0.0625; skME=0.0597; skMF=0.0417;
skMG=0.0318;
dnbGL=1; dnbRL=1; dnbNA=1; dnbMB=1; dnbMC=1; dnbMD=1;
dnbME=1; dnbMF=1; dnbMG=1;
if (vnir == 1) then begin
pMB=0; pMC=0; pMD=0; pME=0; pMF=0; pMG=0;
# GET SHORT GAIN-CODE STRING FROM USER:
p1$ = "Terra ASTER GAIN-CODE:\n";
p2$ = " GAIN-CODE BAND SEQUENCE:\n";
p3$ = " GL GAIN: 1, N, or H: Default=N\n";
p4$ = " RL GAIN: 1, N, or H: Default=N\n";
p5$ = " NA GAIN: 1, N, or H: Default=N\n";
p6$ = " Default GAIN-CODE = NNN\n";
p7$ = " ACCEPT the Default GAIN-CODE,\n";
p8$ = " Or, ENTER Correct GAIN-CODE.\n";
p9$ = " Then, CLICK OK.\n";
p10$ = "GAIN-CODE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
p$ = p$ + p9$ + p10$;
gc$ = "HHN";
t$ = "GAIN-CODE QUERY";
gc$ = PopupString(p$,gc$,t$);
gi$ = mid$(gc$,1,1);
if (gi$ == "1") then skGL=2.2500;
if (gi$ == "H") then skGL=0.6760;
gi$ = mid$(gc$,2,1);
if (gi$ == "1") then skRL=1.8900;
if (gi$ == "H") then skRL=0.7080;
gi$ = mid$(gc$,3,1);
if (gi$ == "1") then skNA=1.1500;
if (gi$ == "H") then skNA=0.4230;
end
if (vnir <> 1) then begin
# GET LONG GAIN-CODE STRING FROM USER:
p1$ = "Terra ASTER GAIN-CODE:\n";
p2$ = " GAIN-CODE BAND SEQUENCE:\n";
p3$ = " GL GAIN: 1, N, or H: Default=N\n";
p4$ = " RL GAIN: 1, N, or H: Default=N\n";
p5$ = " NA GAIN: 1, N, or H: Default=N\n";
p6$ = " MB GAIN: 1, 2, N, or H: Default=N\n";
p7$ = " MC GAIN: 1, 2, N, or H: Default=N\n";
p8$ = " MD GAIN: 1, 2, N, or H: Default=N\n";
p9$ = " ME GAIN: 1, 2, N, or H: Default=N\n";
p10$ = " MF GAIN: 1, 2, N, or H: Default=N\n";
p11$ = " MG GAIN: 1, 2, N, or H: Default=N\n";
p12$ = " Default GAIN-CODE = NNNNNNNNN\n";
p13$ = " ACCEPT the Default GAIN-CODE,\n";
p14$ = " Or, ENTER the Correct GAIN-CODE.\n";
p15$ = " Then, CLICK OK.\n";
p16$ = "GAIN-CODE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
p$ = p$ + p9$ + p10$ + p11$ + p12$ +p13$ + p14$;
p$ = p$ + p15$ + p16$;
gc$ = "HHNNNNNNN";
t$ = "GAIN-CODE QUERY";
gc$ = PopupString(p$,gc$,t$);
gi$ = mid$(gc$,1,1);
if (gi$ == "1") then skGL=2.2500;
if (gi$ == "H") then skGL=0.6760;
gi$ = mid$(gc$,2,1);
if (gi$ == "1") then skRL=1.8900;
if (gi$ == "H") then skRL=0.7080;
gi$ = mid$(gc$,3,1);
if (gi$ == "1") then skNA=1.1500;
if (gi$ == "H") then skNA=0.4230;
gi$ = mid$(gc$,4,1);
if (gi$ == "1") then skMB=0.2900;
if (gi$ == "2") then skMB=0.2900;
if (gi$ == "H") then skMB=0.1087;
gi$ = mid$(gc$,5,1);
if (gi$ == "1") then skMC=0.0925;
if (gi$ == "2") then skMC=0.4090;
if (gi$ == "H") then skMC=0.0348;
gi$ = mid$(gc$,6,1);
if (gi$ == "1") then skMD=0.0830;
if (gi$ == "2") then skMD=0.3900;
if (gi$ == "H") then skMD=0.0313;
gi$ = mid$(gc$,7,1);
if (gi$ == "1") then skME=0.0795;
if (gi$ == "2") then skME=0.3320;
if (gi$ == "H") then skME=0.0299;
gi$ = mid$(gc$,8,1);
if (gi$ == "1") then skMF=0.0556;
if (gi$ == "2") then skMF=0.2450;
if (gi$ == "H") then skMF=0.0209;
gi$ = mid$(gc$,9,1);
if (gi$ == "1") then skMG=0.0424;
if (gi$ == "2") then skMG=0.2650;
if (gi$ == "H") then skMG=0.0159;
end
dssiGL=1829 ; dssiRL=1557 ; dssiNA=1047 ; dssiMB=215.0;
dssiMC=84.67; dssiMD=82.67; dssiME=80.67; dssiMF=78.67;
dssiMG=76.67;
end
if (imager == 13) then begin
p1$ = "WARNING!\n";
p2$ = "SRFI.sml Cannot Process MS Data from\n";
p3$ = "Terra MODIS.\n";
p4$ = "Please Stop this SML by\n";
p5$ = "CLICKing CANCEL.";
p$ = p1$ + p2$ + p3$ + p4$ + p5$;
PopupMessage(p$);
imager$ = "Terra MODIS";
pBL=1; pNA=1; pMA=1; pMB=1; pMC=1; maxi=4200;
dnlow=1; dnhigh=4095;
# wLenBL=0.469; wLenGL=0.555; wLenRL=0.645; wLenNA=0.8585;
# wLenMA=1.240; wLenMB=1.640; wLenMC=2.013;
# skBL=0.1000; skGL=0.1000; skRL=0.1000; skNA=0.1000;
# skMA=0.1000; skMB=0.1000; skMC=0.1000;
# dnbBL=0; dnbGL=0; dnbRL=0; dnbNA=0;
# dnbMA=0; dnbMB=0; dnbMC=0;
# dssiBL=1957; dssiGL=1829; dssiRL=1557; dssiNA=1047;
# dssiMA= 400; dssiMB=215.0; dssiMC=80.67;
end
if (imager == 14) then begin
p1$ = "WARNING!\n";
p2$ = "SRFI.sml Cannot Process MS Data from\n";
p3$ = "Aqua MODIS.\n";
p4$ = "Please Stop this SML by\n";
p5$ = "CLICKing CANCEL.";
p$ = p1$ + p2$ + p3$ + p4$ + p5$;
PopupMessage(p$);
imager$ = "Aqua MODIS";
pBL=1; pNA=1; pMA=1; pMB=1; pMC=1; maxi=4200;
dnlow=1; dnhigh=4095;
# wLenBL=0.469; wLenGL=0.555; wLenRL=0.645; wLenNA=0.8585;
# wLenMA=1.240; wLenMB=1.640; wLenMC=2.013;
# skBL=0.1000; skGL=0.1000; skRL=0.1000; skNA=0.1000;
# skMA=0.1000; skMB=0.1000; skMC=0.1000;
# dnbBL=0; dnbGL=0; dnbRL=0; dnbNA=0;
# dnbMA=0; dnbMB=0; dnbMC=0;
# dssiBL=1957; dssiGL=1829; dssiRL=1557; dssiNA=1047;
# dssiMA= 400; dssiMB=215.0; dssiMC=80.67;
end
printf(" IMAGER: %s\n",imager$);
printf(" IMAGE DN SCALE: %d to %d\n",dnlow,dnhigh);
if (imager == 4) then begin
printf(" GAIN-CODE: %s\n",gc$);
p$ = "LPGS (EOS Data Gateway)";
if (sc == 2) then p$ = "NLAPS (EarthExplorer)";
printf(" DATA VENDOR: %s\n",p$);
end
# ------------------------------------------------------------
# CALCULATE C-FACTORS & FOR ENABLED MS BANDS: Refer to B26.
if (pCB) then cCB = 1 + cRLm1 * (wLenRL/wLenCB)^pc;
if (pBL) then cBL = 1 + cRLm1 * (wLenRL/wLenBL)^pc;
cGL = 1 + cRLm1 * (wLenRL/wLenGL)^pc;
if (pYL) then cYL = 1 + cRLm1 * (wLenRL/wLenYL)^pc;
pc = 1.1;
if (pRE) then cRE = 1 + cRLm1 * (wLenRL/wLenRE)^pc;
if (pNA) then cNA = 1 + cRLm1 * (wLenRL/wLenNA)^pc;
if (pNB) then cNB = 1 + cRLm1 * (wLenRL/wLenNB)^pc;
if (pMA) then cMA = 1 + cRLm1 * (wLenRL/wLenMA)^pc;
if (pMB) then cMB = 1 + cRLm1 * (wLenRL/wLenMB)^pc;
if (pMC) then cMC = 1 + cRLm1 * (wLenRL/wLenMC)^pc;
if (pMD) then cMD = 1 + cRLm1 * (wLenRL/wLenMD)^pc;
if (pME) then begin
cME = 1 + cRLm1 * (wLenRL/wLenME)^pc;
cME = cME * 1.1;
end
if (pMF) then begin
cMF = 1 + cRLm1 * (wLenRL/wLenMF)^pc;
cMF = cMF * 1.5;
end
if (pMG) then begin
cMG = 1 + cRLm1 * (wLenRL/wLenMG)^pc;
cMG = cMG * 2.4;
end
# ------------------------------------------------------------
# PARAMETERS RELATED TO WVI & WBI CALCULATIONS: Refer to B36.
srfiNRint = 254; bslineslope = 1.086;
pbioff = 100; pvioff = 1000;
pfac = 0.2723659; maxpvipbi = 3000;
# ------------------------------------------------------------
# IF atcor < 3, SET C-FACTORS EQUAL TO 1: Refer to B26.
if (atcor < 3) then begin
cCB = 1; cBL = 1; cGL = 1; cYL = 1; cRL = 1; cRE = 1;
cNA = 1; cNB = 1; cMA = 1; cMB = 1; cMC = 1; cMD = 1;
cME = 1; cMF = 1; cMG = 1;
end
# ------------------------------------------------------------
# SET UP ARRAYS: Refer to B22, B28 & B29.
maxip1 = maxi + 1;
array numeric vCB[maxip1]; array numeric vBL[maxip1];
array numeric vGL[maxip1]; array numeric vYL[maxip1];
array numeric vRL[maxip1]; array numeric vRE[maxip1];
array numeric vNA[maxip1]; array numeric vNB[maxip1];
array numeric vMA[maxip1]; array numeric vMB[maxip1];
array numeric vMC[maxip1]; array numeric vMD[maxip1];
array numeric vME[maxip1]; array numeric vMF[maxip1];
array numeric vMG[maxip1]; array numeric v[maxip1];
# ------------------------------------------------------------
# DEFINE proc checkHisto(Band)
# PURPOSE: Check that Input Raster has a Full (Unsampled)
# Histogram. If Not, Then Recompute Histogram.
proc checkHisto (raster Band) begin
local numeric sampleInt;
sampleInt = HistogramGetSampleInterval(Band);
if (sampleInt > 1) begin
DeleteHistogram(Band); CreateHistogram(Band,0);
end
else if (sampleInt == -1) then begin
CreateHistogram(Band,0);
end
end
# ------------------------------------------------------------
# DEFFINE func dnHE(heloc): Refer to B22 * B28.
# PURPOSE: Find Dark Edge of Raster Histogram Based on the
# DN Value at Which the Cumulative Frequency
# Distributoin Rises Rapidly (Faster than HEP).
func dnHE(heloc) begin
local numeric j, jm1, jp1, sumfq, sumcfq, del, dnhe;
sumfq = 0;
for j=1 to maxi begin
sumfq = sumfq + v[j];
end
v[1] = v[1] * 100 / sumfq;
# v[1] equals percent frequency for DN = 1
for j=2 to maxi begin
jm1 = j - 1; v[j] = v[jm1] + v[j] * 100 / sumfq;
# v[j] equals percent frequency for DN = j
end
dnhe = 0;
for j=1 to maxi begin
jp1 = j + 1; del = v[jp1] - v[j];
if (dnhe == 0) then begin
if (del > heloc) then dnhe = jp1;
end
end
return(dnhe);
end
# ------------------------------------------------------------
# IF atcor = 3, RESCALE C-FACTORS BASED ON msfac:
# Refer to B23.
if (atcor == 3) then begin
if (pCB) then cCB = cCB * msfac;
if (pBL) then cBL = cBL * msfac;
cGL = cGL * msfac; cRL = cRL * msfac;
if (pYL) then cYL = cYL * msfac;
if (pRE) then cRE = cRE * msfac;
if (pNA) then cNA = cNA * msfac;
if (pNB) then cNB = cNB * msfac;
if (pMA) then cMA = cMA * msfac;
if (pMB) then cMB = cMB * msfac;
if (pMC) then cMC = cMC * msfac;
if (pMD) then cMD = cMD * msfac;
if (pME) then cME = cME * msfac;
if (pMF) then cMF = cMF * msfac;
if (pMG) then cMG = cMG * msfac;
end
# ------------------------------------------------------------
# OPEN APPROPRIATE INPUT RASTERS & GET MIN & MAX VALUES:
# Refer to B27.
printf("\nOPEN INPUT RASTERS:\n");
if (pCB) then begin
printf(" CB");
GetInputRaster(CB); nlins = NumLins(CB);
ncols = NumCols(CB); dntype$ = RastType(CB);
checkHisto(CB);
dnminCB = GlobalMin(CB); dnmaxCB = GlobalMax(CB);
end
if (pBL) then begin
if (pCB) then begin
if (pBL) then begin
printf(" BL");
GetInputRaster(BL,nlins,ncols,dntype$);
end
end
else begin
printf(" BL");
GetInputRaster(BL); nlins = NumLins(BL);
ncols = NumCols(BL); dntype$ = RastType(BL);
end
checkHisto(BL);
dnminBL = GlobalMin(BL); dnmaxBL = GlobalMax(BL);
end
printf(" GL");
if (pCB or pBL) then begin
GetInputRaster(GL,nlins,ncols,dntype$);
end
else begin
GetInputRaster(GL);
nlins = NumLins(GL); ncols = NumCols(GL);
dntype$ = RastType(GL);
end
checkHisto(GL);
dnminGL = GlobalMin(GL); dnmaxGL = GlobalMax(GL);
if (pYL) then begin
printf(" YL");
GetInputRaster(YL,nlins,ncols,dntype$);
checkHisto(YL);
dnminYL = GlobalMin(YL); dnmaxYL = GlobalMax(YL);
end
printf(" RL");
GetInputRaster(RL,nlins,ncols,dntype$);
checkHisto(RL);
dnminRL = GlobalMin(RL); dnmaxRL = GlobalMax(RL);
if (pRE) then begin
printf(" RE");
GetInputRaster(RE,nlins,ncols,dntype$);
checkHisto(RE);
dnminRE = GlobalMin(RE); dnmaxRE = GlobalMax(RE);
end
if (pNA) then begin
printf(" NA");
GetInputRaster(NA,nlins,ncols,dntype$);
checkHisto(NA);
dnminNA = GlobalMin(NA); dnmaxNA = GlobalMax(NA);
end
if (pNB) then begin
printf(" NB");
GetInputRaster(NB,nlins,ncols,dntype$);
checkHisto(NB);
dnminNB = GlobalMin(NB); dnmaxNB = GlobalMax(NB);
end
if (pMA) then begin
printf(" MA");
GetInputRaster(MA,nlins,ncols,dntype$);
checkHisto(MA);
dnminMA = GlobalMin(MA); dnmaxMA = GlobalMax(MA);
end
if (pMB) then begin
printf(" MB");
GetInputRaster(MB,nlins,ncols,dntype$);
checkHisto(MB);
dnminMB = GlobalMin(MB); dnmaxMB = GlobalMax(MB);
end
if (pMC) then begin
printf(" MC");
GetInputRaster(MC,nlins,ncols,dntype$);
checkHisto(MC);
dnminMC = GlobalMin(MC); dnmaxMC = GlobalMax(MC);
end
if (pMD) then begin
printf(" MD");
GetInputRaster(MD,nlins,ncols,dntype$);
checkHisto(MD);
dnminMD = GlobalMin(MD); dnmaxMD = GlobalMax(MD);
end
if (pME) then begin
printf(" ME");
GetInputRaster(ME,nlins,ncols,dntype$);
checkHisto(ME);
dnminME = GlobalMin(ME); dnmaxME = GlobalMax(ME);
end
if (pMF) then begin
printf(" MF");
GetInputRaster(MF,nlins,ncols,dntype$);
checkHisto(MF);
dnminMF = GlobalMin(MF); dnmaxMF = GlobalMax(MF);
end
if (pMG) then begin
printf(" MG");
GetInputRaster(MG,nlins,ncols,dntype$);
checkHisto(MG);
dnminMG = GlobalMin(MG); dnmaxMG = GlobalMax(MG);
end
printf("\n");
# ------------------------------------------------------------
# IF atcor > 1, FIND DARK-EDGES OF DN HISTOGRAMS:
# Refer to B28.
if (atcor > 1) then begin
if (pCB) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(CB,v);
dnheCB = dnHE(he);
end
if (pBL) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(BL,v);
dnheBL = dnHE(he);
if (dnheBL < dnbBL) then dnheBL = ceil(dnbBL);
end
ResizeArrayClear(v,maxip1); ReadHistogram(GL,v);
dnheGL = dnHE(he);
if (dnheGL < dnbGL) then dnheGL = ceil(dnbGL);
if (pYL) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(YL,v);
dnheYL = dnHE(he);
if (dnheYL < dnbYL) then dnheYL = ceil(dnbYL);
end
ResizeArrayClear(v,maxip1); ReadHistogram(RL,v);
dnheRL = dnHE(he);
if (dnheRL < dnbRL) then dnheRL = ceil(dnbRL);
if (pRE) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(RE,v);
dnheRE = dnHE(he);
if (dnheRE < dnbRE) then dnheRE = ceil(dnbRE);
end
if (pNA) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(NA,v);
dnheNA = dnHE(he);
if (dnheNA < dnbNA) then dnheNA = ceil(dnbNA);
end
if (pNB) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(NB,v);
dnheNB = dnHE(he);
if (dnheNB < dnbNB) then dnheNB = ceil(dnbNB);
end
if (pMA) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(MA,v);
dnheMA = dnHE(he);
if (dnheMA < dnbMA) then dnheMA = ceil(dnbMA);
end
if (pMB) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(MB,v);
dnheMB = dnHE(he);
if (dnheMB < dnbMB) then dnheMB = ceil(dnbMB);
end
if (pMC) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(MC,v);
dnheMC = dnHE(he);
if (dnheMC < dnbMC) then dnheMC = ceil(dnbMC);
end
if (pMD) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(MD,v);
dnheMD = dnHE(he);
if (dnheMD < dnbMD) then dnheMD = ceil(dnbMD);
end
if (pME) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(ME,v);
dnheME = dnHE(he);
if (dnheME < dnbME) then dnheME = ceil(dnbME);
end
if (pMF) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(MF,v);
dnheMF = dnHE(he);
if (dnheMF < dnbMF) then dnheMF = ceil(dnbMF);
end
if (pMG) then begin
ResizeArrayClear(v,maxip1); ReadHistogram(MG,v);
dnheMG = dnHE(he);
if (dnheMG < dnbMG) then dnheMG = ceil(dnbMG);
end
end
# ------------------------------------------------------------
# CALCULATE SItoa & DN-to-SRF CONVERSION-COEFFICIENTS: aBAND
# Refer to B29.
esd = 1 - 0.01672 * cosd(0.9856 * (doy - 4));
sitoafac = sind(sunelevang) / esd^2;
if (pCB) then sitoaCB = dssiCB * sitoafac;
if (pBL) then sitoaBL = dssiBL * sitoafac;
sitoaGL = dssiGL * sitoafac; sitoaRL = dssiRL * sitoafac;
if (pYL) then sitoaYL = dssiYL * sitoafac;
if (pRE) then sitoaRE = dssiRE * sitoafac;
if (pNA) then sitoaNA = dssiNA * sitoafac;
if (pNB) then sitoaNB = dssiNB * sitoafac;
if (pMA) then sitoaMA = dssiMA * sitoafac;
if (pMB) then sitoaMB = dssiMB * sitoafac;
if (pMC) then sitoaMC = dssiMC * sitoafac;
if (pMD) then sitoaMD = dssiMD * sitoafac;
if (pME) then sitoaME = dssiME * sitoafac;
if (pMF) then sitoaMF = dssiMF * sitoafac;
if (pMG) then sitoaMG = dssiMG * sitoafac;
pi100 = pi * 100;
if (pCB) then aCB = pi100 * skCB / sitoaCB;
if (pBL) then aBL = pi100 * skBL / sitoaBL;
aGL = pi100 * skGL / sitoaGL; aRL = pi100 * skRL / sitoaRL;
if (pYL) then aYL = pi100 * skYL / sitoaYL;
if (pRE) then aRE = pi100 * skRE / sitoaRE;
if (pNA) then aNA = pi100 * skNA / sitoaNA;
if (pNB) then aNB = pi100 * skNB / sitoaNB;
if (pMA) then aMA = pi100 * skMA / sitoaMA;
if (pMB) then aMB = pi100 * skMB / sitoaMB;
if (pMC) then aMC = pi100 * skMC / sitoaMC;
if (pMD) then aMD = pi100 * skMD / sitoaMD;
if (pME) then aME = pi100 * skME / sitoaME;
if (pMF) then aMF = pi100 * skMF / sitoaMF;
if (pMG) then aMG = pi100 * skMG / sitoaMG;
# ------------------------------------------------------------
# IF atcor > 1, ESTIMATE SRFpath & DNpath: Refer to B29.
if (atcor > 1) then begin
dnpathGL1 = dnheGL;
srfpathGL1 = (dnpathGL1 - dnbGL) * aGL;
dnpathRL = dnheRL;
srfpathRL = (dnpathRL - dnbRL) * aRL;
logsrfpathRL = log(srfpathRL);
wLenV = wLenGL; srfpathV = srfpathGL1;
if (pBL) then begin
wtBL = 0.55114;
dnpathBL1 = dnheBL;
srfpathBL1 = (dnpathBL1 - dnbBL) * aBL;
wLenV = wtBL * wLenBL + (1 - wtBL) * wLenGL;
srfpathV = (srfpathBL1 + srfpathGL1) / 2;
end
logwLenRLdwLenV = log(wLenRL / wLenV);
logsrfpathV = log(srfpathV);
difoflogs = logsrfpathV - logsrfpathRL;
p = difoflogs / logwLenRLdwLenV;
if (pCB) then begin
srfpathCB2 = srfpathRL * (wLenRL / wLenCB)^p;
srfpathCB = SetMin(srfpathCB1, srfpathCB2);
dnpathCB = round(srfpathCB / aCB) + dnbCB;
end
if (pBL) then begin
srfpathBL2 = srfpathRL * (wLenRL / wLenBL)^p;
srfpathBL = SetMin(srfpathBL1, srfpathBL2);
dnpathBL = round(srfpathBL / aBL) + dnbBL;
end
srfpathGL2 = srfpathRL * (wLenRL / wLenGL)^p;
srfpathGL = SetMin(srfpathGL1, srfpathGL2);
dnpathGL = round(srfpathGL / aGL) + dnbGL;
if (pYL) then begin
srfpathYL2 = srfpathRL * (wLenRL / wLenYL)^p;
dnpathYL1 = dnheYL;
srfpathYL1 = (dnpathYL1 - dnbYL) * aYL;
srfpathYL = SetMin(srfpathYL1, srfpathYL2);
dnpathYL = round(srfpathYL / aYL) + dnbYL;
end
p2 = 0.2 * p;
if (pRE) then begin
srfpathRE2 = srfpathRL * (wLenRL / wLenRE)^p2;
dnpathRE1 = dnheRE;
srfpathRE1 = (dnpathRE1 - dnbRE) * aRE;
srfpathRE = SetMin(srfpathRE1, srfpathRE2);
dnpathRE = round(srfpathRE / aRE) + dnbRE;
end
if (pNA) then begin
srfpathNA2 = srfpathRL * (wLenRL / wLenNA)^p2;
dnpathNA1 = dnheNA;
srfpathNA1 = (dnpathNA1 - dnbNA) * aNA;
srfpathNA = SetMin(srfpathNA1, srfpathNA2);
dnpathNA = round(srfpathNA / aNA) + dnbNA;
end
if (pNB) then begin
srfpathNB2 = srfpathRL * (wLenRL / wLenNB)^p2;
dnpathNB1 = dnheNB;
srfpathNB1 = (dnpathNB1 - dnbNB) * aNB;
srfpathNB = SetMin(srfpathNB1, srfpathNB2);
dnpathNB = round(srfpathNB / aNB) + dnbNB;
end
if (pMA) then begin
srfpathMA2 = srfpathRL * (wLenRL / wLenMA)^p2;
dnpathMA1 = dnheMA;
srfpathMA1 = (dnpathMA1 - dnbMA) * aMA;
srfpathMA = SetMin(srfpathMA1, srfpathMA2);
dnpathMA = round(srfpathMA / aMA) + dnbMA;
end
if (pMB) then begin
srfpathMB2 = srfpathRL * (wLenRL / wLenMB)^p2;
dnpathMB1 = dnheMB;
srfpathMB1 = (dnpathMB1 - dnbMB) * aMB;
srfpathMB = SetMin(srfpathMB1, srfpathMB2);
dnpathMB = round(srfpathMB / aMB) + dnbMB;
end
if (pMC) then begin
srfpathMC2 = srfpathRL * (wLenRL / wLenMC)^p2;
dnpathMC1 = dnheMC;
srfpathMC1 = (dnpathMC1 - dnbMC) * aMC;
srfpathMC = SetMin(srfpathMC1, srfpathMC2);
dnpathMC = round(srfpathMC / aMC) + dnbMC;
end
if (pMD) then begin
srfpathMD2 = srfpathRL * (wLenRL / wLenMD)^p2;
dnpathMD1 = dnheMD;
srfpathMD1 = (dnpathMD1 - dnbMD) * aMD;
srfpathMD = SetMin(srfpathMD1, srfpathMD2);
dnpathMD = round(srfpathMD / aMD) + dnbMD;
end
if (pME) then begin
srfpathME2 = srfpathRL * (wLenRL / wLenME)^p2;
dnpathME1 = dnheME;
srfpathME1 = (dnpathME1 - dnbME) * aME;
srfpathME = SetMin(srfpathME1, srfpathME2);
dnpathME = round(srfpathME / aME) + dnbME;
end
if (pMF) then begin
srfpathMF2 = srfpathRL * (wLenRL / wLenMF)^p2;
dnpathMF1 = dnheMF;
srfpathMF1 = (dnpathMF1 - dnbMF) * aMF;
srfpathMF = SetMin(srfpathMF1, srfpathMF2);
dnpathMF = round(srfpathMF / aMF) + dnbMF;
end
if (pMG) then begin
srfpathMG2 = srfpathRL * (wLenRL / wLenMG)^p2;
dnpathMG1 = dnheMG;
srfpathMG1 = (dnpathMG1 - dnbMG) * aMG;
srfpathMG = SetMin(srfpathMG1, srfpathMG2);
dnpathMG = round(srfpathMG / aMG) + dnbMG;
end
end
# ------------------------------------------------------------
# SET UP LOOK-UP TABLES FOR SRFI VALUES (FUNCTION OF DNs):
# Refer to B30.
if (atcor == 1) then begin
srfpathCB = 0; srfpathBL = 0; srfpathGL = 0;
srfpathYL = 0; srfpathRL = 0; srfpathRE = 0;
srfpathNA = 0; srfpathNB = 0; srfpathMA = 0;
srfpathMB = 0; srfpathMC = 0; srfpathMD = 0;
srfpathME = 0; srfpathMF = 0; srfpathMG = 0;
end
if (cCB) then ResizeArrayClear(vCB,maxip1);
if (pBL) then ResizeArrayClear(vBL,maxip1);
ResizeArrayClear(vGL,maxip1); ResizeArrayClear(vRL,maxip1);
if (pYL) then ResizeArrayClear(vYL,maxip1);
if (pRE) then ResizeArrayClear(vRE,maxip1);
if (pNA) then ResizeArrayClear(vNA,maxip1);
if (pNB) then ResizeArrayClear(vNB,maxip1);
if (pMA) then ResizeArrayClear(vMA,maxip1);
if (pMB) then ResizeArrayClear(vMB,maxip1);
if (pMC) then ResizeArrayClear(vMC,maxip1);
if (pMD) then ResizeArrayClear(vMD,maxip1);
if (pME) then ResizeArrayClear(vME,maxip1);
if (pMF) then ResizeArrayClear(vMF,maxip1);
if (pMG) then ResizeArrayClear(vMG,maxip1);
for i=1 to maxi begin
if (pCB) then begin
srftoaCB = (i - dnbCB) * aCB;
srfapcCB = srftoaCB - srfpathCB;
srfsfcCB = srfapcCB * cCB;
srfiCB = round(srfsfcCB * 100);
if (srfiCB < 1) then srfiCB = 1;
vCB[i] = srfiCB;
end
if (pBL) then begin
srftoaBL = (i - dnbBL) * aBL;
srfapcBL = srftoaBL - srfpathBL;
srfsfcBL = srfapcBL * cBL;
srfiBL = round(srfsfcBL * 100);
if (srfiBL < 1) then srfiBL = 1;
vBL[i] = srfiBL;
end
srftoaGL = (i - dnbGL) * aGL;
srfapcGL = srftoaGL - srfpathGL;
srfsfcGL = srfapcGL * cGL;
srfiGL = round(srfsfcGL * 100);
if (srfiGL < 1) then srfiGL = 1;
vGL[i] = srfiGL;
if (pYL) then begin
srftoaYL = (i - dnbYL) * aYL;
srfapcYL = srftoaYL - srfpathYL;
srfsfcYL = srfapcYL * cYL;
srfiYL = round(srfsfcYL * 100);
if (srfiYL < 1) then srfiYL = 1;
vYL[i] = srfiYL;
end
srftoaRL = (i - dnbRL) * aRL;
srfapcRL = srftoaRL - srfpathRL;
srfsfcRL = srfapcRL * cRL;
srfiRL = round(srfsfcRL * 100);
if (srfiRL < 1) then srfiRL = 1;
vRL[i] = srfiRL;
if (pRE) then begin
srftoaRE = (i - dnbRE) * aRE;
srfapcRE = srftoaRE - srfpathRE;
srfsfcRE = srfapcRE * cRE;
srfiRE = round(srfsfcRE * 100);
if (srfiRE < 1) then srfiRE = 1;
vRE[i] = srfiRE;
end
if (pNA) then begin
srftoaNA = (i - dnbNA) * aNA;
srfapcNA = srftoaNA - srfpathNA;
srfsfcNA = srfapcNA * cNA;
srfiNA = round(srfsfcNA * 100);
if (srfiNA < 1) then srfiNA = 1;
vNA[i] = srfiNA;
end
if (pNB) then begin
srftoaNB = (i - dnbNB) * aNB;
srfapcNB = srftoaNB - srfpathNB;
srfsfcNB = srfapcNB * cNB;
srfiNB = round(srfsfcNB * 100);
if (srfiNB < 1) then srfiNB = 1;
vNB[i] = srfiNB;
end
if (pMA) then begin
srftoaMA = (i - dnbMA) * aMA;
srfapcMA = srftoaMA - srfpathMA;
srfsfcMA = srfapcMA * cMA;
srfiMA = round(srfsfcMA * 100);
if (srfiMA < 1) then srfiMA = 1;
vMA[i] = srfiMA;
end
if (pMB) then begin
srftoaMB = (i - dnbMB) * aMB;
srfapcMB = srftoaMB - srfpathMB;
srfsfcMB = srfapcMB * cMB;
srfiMB = round(srfsfcMB * 100);
if (srfiMB < 1) then srfiMB = 1;
vMB[i] = srfiMB;
end
if (pMC) then begin
srftoaMC = (i - dnbMC) * aMC;
srfapcMC = srftoaMC - srfpathMC;
srfsfcMC = srfapcMC * cMC;
srfiMC = round(srfsfcMC * 100);
if (srfiMC < 1) then srfiMC = 1;
vMC[i] = srfiMC;
end
if (pMD) then begin
srftoaMD = (i - dnbMD) * aMD;
srfapcMD = srftoaMD - srfpathMD;
srfsfcMD = srfapcMD * cMD;
srfiMD = round(srfsfcMD * 100);
if (srfiMD < 1) then srfiMD = 1;
vMD[i] = srfiMD;
end
if (pME) then begin
srftoaME = (i - dnbME) * aME;
srfapcME = srftoaME - srfpathME;
srfsfcME = srfapcME * cME;
srfiME = round(srfsfcME * 100);
if (srfiME < 1) then srfiME = 1;
vME[i] = srfiME;
end
if (pMF) then begin
srftoaMF = (i - dnbMF) * aMF;
srfapcMF = srftoaMF - srfpathMF;
srfsfcMF = srfapcMF * cMF;
srfiMF = round(srfsfcMF * 100);
if (srfiMF < 1) then srfiMF = 1;
vMF[i] = srfiMF;
end
if (pMG) then begin
srftoaMG = (i - dnbMG) * aMG;
srfapcMG = srftoaMG - srfpathMG;
srfsfcMG = srfapcMG * cMG;
srfiMG = round(srfsfcMG * 100);
if (srfiMG < 1) then srfiMG = 1;
vMG[i] = srfiMG;
end
end
# ------------------------------------------------------------
# IF atcor > 1, THEN GENERATE A HISTOGRAM ANALYSES REPORT:
# Refer to B31.
if (atcor > 1) then begin
printf("\nHISTOGRAM ANALYSES REPORT:\n");
printf("BAND DNmin DNhe DNmax");
printf(" DNpath SRFpath1 SRFpath\n");
if (pCB) then begin
printf(" CB %3d %3d",dnminCB,dnheCB);
printf(" %4d %3d",dnmaxCB, dnpathCB);
printf(" %8.3f %7.3f\n",srfpathCB1,srfpathCB);
end
if (pBL) then begin
printf(" BL %3d %3d",dnminBL,dnheBL);
printf(" %4d %3d",dnmaxBL, dnpathBL);
printf(" %8.3f %7.3f\n",srfpathBL1,srfpathBL);
end
printf(" GL %3d %3d",dnminGL,dnheGL);
printf(" %4d %3d",dnmaxGL, dnpathGL);
printf(" %8.3f %7.3f\n",srfpathGL1,srfpathGL);
if (pYL) then begin
printf(" YL %3d %3d",dnminYL,dnheYL);
printf(" %4d %3d",dnmaxYL, dnpathYL);
printf(" %8.3f %7.3f\n",srfpathYL1,srfpathYL);
end
printf(" RL %3d %3d",dnminRL,dnheRL);
printf(" %4d %3d",dnmaxRL, dnpathRL);
printf(" %8.3f %7.3f\n",srfpathRL,srfpathRL);
if (pRE) then begin
printf(" RE %3d %3d",dnminRE,dnheRE);
printf(" %4d %3d",dnmaxRE, dnpathRE);
printf(" %8.3f %7.3f\n",srfpathRE1,srfpathRE);
end
if (pNA) then begin
printf(" NA %3d %3d",dnminNA,dnheNA);
printf(" %4d %3d",dnmaxNA, dnpathNA);
printf(" %8.3f %7.3f\n",srfpathNA1,srfpathNA);
end
if (pNB) then begin
printf(" NB %3d %3d",dnminNB,dnheNB);
printf(" %4d %3d",dnmaxNB, dnpathNB);
printf(" %8.3f %7.3f\n",srfpathNB1,srfpathNB);
end
if (pMA) then begin
printf(" MA %3d %3d",dnminMA,dnheMA);
printf(" %4d %3d",dnmaxMA, dnpathMA);
printf(" %8.3f %7.3f\n",srfpathMA1,srfpathMA);
end
if (pMB) then begin
printf(" MB %3d %3d",dnminMB,dnheMB);
printf(" %4d %3d",dnmaxMB, dnpathMB);
printf(" %8.3f %7.3f\n",srfpathMB1,srfpathMB);
end
if (pMC) then begin
printf(" MC %3d %3d",dnminMC,dnheMC);
printf(" %4d %3d",dnmaxMC, dnpathMC);
printf(" %8.3f %7.3f\n",srfpathMC1,srfpathMC);
end
if (pMD) then begin
printf(" MD %3d %3d",dnminMD,dnheMD);
printf(" %4d %3d",dnmaxMD, dnpathMD);
printf(" %8.3f %7.3f\n",srfpathMD1,srfpathMD);
end
if (pME) then begin
printf(" ME %3d %3d",dnminME,dnheME);
printf(" %4d %3d",dnmaxME, dnpathME);
printf(" %8.3f %7.3f\n",srfpathME1,srfpathME);
end
if (pMF) then begin
printf(" MF %3d %3d",dnminMF,dnheMF);
printf(" %4d %3d",dnmaxMF, dnpathMF);
printf(" %8.3f %7.3f\n",srfpathMF1,srfpathMF);
end
if (pMG) then begin
printf(" MG %3d %3d",dnminMG,dnheMG);
printf(" %4d %3d",dnmaxMG, dnpathMG);
printf(" %8.3f %7.3f\n",srfpathMG1,srfpathMG);
end
printf("SRFIpath Model Parameter: %6.4f\n",p);
# ---------------------------------------------------------
# PRINT HISTOGRAM ANALYSIS REPORT TO CONSOLE WINDOW:
printf("\n");
p1$ = "HISTOGRAM ANALYSES REPORT ASSESSMENT\n";
p2$ = "* EXAMINE the HISTOGRAM ANALYSES REPORT.\n";
p3$ = " - DNheXX Should be Somewhat Higher";
p4$ = " than the DNminXX.\n";
p5$ = " - SRFpath1XX Should be Like SRFpathXX\n";
p6$ = " EXCEPT for NA and/or NB Band.\n";
p7$ = " - SRFpath Should Decrease Band by Band.\n";
p8$ = " - SRFIpath Model Parameter Value Should\n";
p9$ = " be in the RANGE: 1.5000 to 5.0000.\n\n";
p10$ = "* If the REPORT Looks Reasonable, Continue.\n";
p11$ = "* If the REPORT Looks Unreasonable, STOP.\n\n";
p12$ = " - TO STOP (in the NEXT Popup WINDOW),\n";
p13$ = " * CLICK CANCEL.\n";
p14$ = " * Then, Re-Run SRFI.sml with a\n";
p15$ = " Better delcf Value.\n\n";
p16$ = " - To CLOSE This Window, CLICK OK.";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
p$ = p$ + p9$ + p10$ + p11$ + p12$ + p13$ + p14$;
p$ = p$ + p15$ + p16$;
PopupMessage(p$);
end
# ------------------------------------------------------------
# ASSIGN NAMES & LOCATIONS FOR OUTPUT RASTERS:
# Refer to B33.
printf("ASSIGN NAMES & LOCATIONS FOR OUTPUT RASTERS:\n");
srfitype$ = "16-bit unsigned"; ptype$ = "16-bit unsigned";
if (pCB) then begin
printf(" SRFICB");
GetOutputRaster(SRFICB,nlins,ncols,srfitype$);
SetNull(SRFICB,0);
end
if (pBL) then begin
printf(" SRFIBL");
GetOutputRaster(SRFIBL,nlins,ncols,srfitype$);
SetNull(SRFIBL,0);
end
printf(" SRFIGL");
GetOutputRaster(SRFIGL,nlins,ncols,srfitype$);
SetNull(SRFIGL,0);
if (pYL) then begin
printf(" SRFIYL");
GetOutputRaster(SRFIYL,nlins,ncols,srfitype$);
SetNull(SRFIYL,0);
end
printf(" SRFIRL");
GetOutputRaster(SRFIRL,nlins,ncols,srfitype$);
SetNull(SRFIRL,0);
if (pRE) then begin
printf(" SRFIRE");
GetOutputRaster(SRFIRE,nlins,ncols,srfitype$);
SetNull(SRFIRE,0);
end
if (pNA) then begin
printf(" SRFINA");
GetOutputRaster(SRFINA,nlins,ncols,srfitype$);
SetNull(SRFINA,0);
end
if (pNB) then begin
printf(" SRFINB");
GetOutputRaster(SRFINB,nlins,ncols,srfitype$);
SetNull(SRFINB,0);
end
if (pMA) then begin
printf(" SRFIMA");
GetOutputRaster(SRFIMA,nlins,ncols,srfitype$);
SetNull(SRFIMA,0);
end
if (pMB) then begin
printf(" SRFIMB");
GetOutputRaster(SRFIMB,nlins,ncols,srfitype$);
SetNull(SRFIMB,0);
end
if (pMC) then begin
printf(" SRFIMC");
GetOutputRaster(SRFIMC,nlins,ncols,srfitype$);
SetNull(SRFIMC,0);
end
if (pMD) then begin
printf(" SRFIMD");
GetOutputRaster(SRFIMD,nlins,ncols,srfitype$);
SetNull(SRFIMD,0);
end
if (pME) then begin
printf(" SRFIME");
GetOutputRaster(SRFIME,nlins,ncols,srfitype$);
SetNull(SRFIME,0);
end
if (pMF) then begin
printf(" SRFIMF");
GetOutputRaster(SRFIMF,nlins,ncols,srfitype$);
SetNull(SRFIMF,0);
end
if (pMG) then begin
printf(" SRFIMG");
GetOutputRaster(SRFIMG,nlins,ncols,srfitype$);
SetNull(SRFIMG,0);
end
if (atcor == 3) then begin
printf(" PVI");
GetOutputRaster(PVI,nlins,ncols,ptype$);
printf(" PBI");
GetOutputRaster(PBI,nlins,ncols,ptype$);
SetNull(PVI,0); SetNull(PBI,0);
end
printf("\n\n");
# ------------------------------------------------------------
# FILL SRFI RASTERS WITH VALUES: Refer to B35.
if (atcor == 1) then begin
printf("PRODUCING SRFItoa RASTERS:\n");
end
if (atcor == 2) then begin
printf("PRODUCING SRFIapc RASTERS:\n");
end
if (atcor == 3) then begin
printf("PRODUCING SRFIsfc RASTERS:\n");
end
if (pCB) then begin
printf(" CB");
for each CB begin
dnCB = CB; nc = IsNull(dnCB);
if (nc==0) then SRFICB = vCB[dnCB];
end
CreateHistogram(SRFICB,0); CreatePyramid(SRFICB,0);
CopySubobjects(RL,SRFICB,"GEOREF");
smlContrast.InputUpperLimit = 2820;
smlContrast.Compute(SRFICB,clut1$,clut2$,clut3$);
CloseRaster(CB); CloseRaster(SRFICB);
end
if (pBL) then begin
printf(" BL");
for each BL begin
dnBL = BL; nc = IsNull(dnBL);
if (nc==0) then SRFIBL = vBL[dnBL];
end
CreateHistogram(SRFIBL,0); CreatePyramid(SRFIBL,0);
CopySubobjects(RL,SRFIBL,"GEOREF");
smlContrast.InputUpperLimit = 2820;
smlContrast.Compute(SRFIBL,clut1$,clut2$,clut3$);
CloseRaster(BL); CloseRaster(SRFIBL);
end
printf(" GL");
for each GL begin
dnGL = GL; nc = IsNull(dnGL);
if (nc==0) then SRFIGL = vGL[dnGL];
end
CreateHistogram(SRFIGL,0); CreatePyramid(SRFIGL,0);
CopySubobjects(RL,SRFIGL,"GEOREF");
smlContrast.InputUpperLimit = 3230;
smlContrast.Compute(SRFIGL,clut1$,clut2$,clut3$);
CloseRaster(GL); CloseRaster(SRFIGL);
if (pYL) then begin
printf(" YL");
for each YL begin
dnYL = YL; nc = IsNull(dnYL);
if (nc==0) then SRFIYL = vYL[dnYL];
end
CreateHistogram(SRFIYL,0); CreatePyramid(SRFIYL,0);
CopySubobjects(RL,SRFIYL,"GEOREF");
smlContrast.InputUpperLimit = 3230;
smlContrast.Compute(SRFIYL,clut1$,clut2$,clut3$);
CloseRaster(YL); CloseRaster(SRFIYL);
end
printf(" RL");
for each RL begin
dnRL = RL; nc = IsNull(dnRL);
if (nc==0) then SRFIRL = vRL[dnRL];
end
CreateHistogram(SRFIRL,0); CreatePyramid(SRFIRL,0);
CopySubobjects(RL,SRFIRL,"GEOREF");
smlContrast.InputUpperLimit = 3500;
smlContrast.Compute(SRFIRL,clut1$,clut2$,clut3$);
CloseRaster(RL);
if (pRE) then begin
printf(" RE");
for each RE begin
dnRE = RE; nc = IsNull(dnRE);
if (nc==0) then SRFIRE = vRE[dnRE];
end
CreateHistogram(SRFIRE,0); CreatePyramid(SRFIRE,0);
CopySubobjects(RL,SRFIRE,"GEOREF");
smlContrast.InputUpperLimit = 4000;
smlContrast.Compute(SRFIRE,clut1$,clut2$,clut3$);
CloseRaster(RE); CloseRaster(SRFIRE);
end
if (pNA) then begin
printf(" NA");
for each NA begin
dnNA = NA; nc = IsNull(dnNA);
if (nc==0) then SRFINA = vNA[dnNA];
end
CreateHistogram(SRFINA,0); CreatePyramid(SRFINA,0);
CopySubobjects(RL,SRFINA,"GEOREF");
smlContrast.InputUpperLimit = 5000;
smlContrast.Compute(SRFINA,clut1$,clut2$,clut3$);
CloseRaster(NA);
end
if (pNB) then begin
printf(" NB");
for each NB begin
dnNA = NB; nc = IsNull(dnNB);
if (nc==0) then SRFINB = vNB[dnNB];
end
CreateHistogram(SRFINB,0); CreatePyramid(SRFINB,0);
CopySubobjects(RL,SRFINB,"GEOREF");
smlContrast.InputUpperLimit = 5000;
smlContrast.Compute(SRFINB,clut1$,clut2$,clut3$);
CloseRaster(NB);
end
if (pMA) then begin
printf(" MA");
for each MA begin
dnMA = MA; nc = IsNull(dnMA);
if (nc==0) then SRFIMA = vMA[dnMA];
end
CreateHistogram(SRFIMA,0); CreatePyramid(SRFIMA,0);
CopySubobjects(RL,SRFIMA,"GEOREF");
smlContrast.InputUpperLimit = 4200;
smlContrast.Compute(SRFIMA,clut1$,clut2$,clut3$);
CloseRaster(MA); CloseRaster(SRFIMA);
end
if (pMB) then begin
printf(" MB");
for each MB begin
dnMB = MB; nc = IsNull(dnMB);
if (nc==0) then SRFIMB = vMB[dnMB];
end
CreateHistogram(SRFIMB,0); CreatePyramid(SRFIMB,0);
CopySubobjects(RL,SRFIMB,"GEOREF");
smlContrast.InputUpperLimit = 4200;
smlContrast.Compute(SRFIMB,clut1$,clut2$,clut3$);
CloseRaster(MB); CloseRaster(SRFIMB);
end
if (pMC) then begin
printf(" MC");
for each MC begin
dnMC = MC; nc = IsNull(dnMC);
if (nc==0) then SRFIMC = vMC[dnMC];
end
CreateHistogram(SRFIMC,0); CreatePyramid(SRFIMC,0);
CopySubobjects(RL,SRFIMC,"GEOREF");
smlContrast.InputUpperLimit = 3500;
smlContrast.Compute(SRFIMC,clut1$,clut2$,clut3$);
CloseRaster(MC); CloseRaster(SRFIMC);
end
if (pMD) then begin
printf(" MD");
for each MD begin
dnMD = MD; nc = IsNull(dnMD);
if (nc==0) then SRFIMD = vMD[dnMD];
end
CreateHistogram(SRFIMD,0); CreatePyramid(SRFIMD,0);
CopySubobjects(RL,SRFIMD,"GEOREF");
smlContrast.InputUpperLimit = 3500;
smlContrast.Compute(SRFIMD,clut1$,clut2$,clut3$);
CloseRaster(MD); CloseRaster(SRFIMD);
end
if (pME) then begin
printf(" ME");
for each ME begin
dnME = ME; nc = IsNull(dnME);
if (nc==0) then SRFIME = vME[dnME];
end
CreateHistogram(SRFIME,0); CreatePyramid(SRFIME,0);
CopySubobjects(RL,SRFIME,"GEOREF");
smlContrast.InputUpperLimit = 3500;
smlContrast.Compute(SRFIME,clut1$,clut2$,clut3$);
CloseRaster(ME); CloseRaster(SRFIME);
end
if (pMF) then begin
printf(" MF");
for each MF begin
dnMF = MF; nc = IsNull(dnMF);
if (nc==0) then SRFIMF = vMF[dnMF];
end
CreateHistogram(SRFIMF,0); CreatePyramid(SRFIMF,0);
CopySubobjects(RL,SRFIMF,"GEOREF");
smlContrast.InputUpperLimit = 3500;
smlContrast.Compute(SRFIMF,clut1$,clut2$,clut3$);
CloseRaster(MF); CloseRaster(SRFIMF);
end
if (pMG) then begin
printf(" MG");
for each MG begin
dnMG = MG; nc = IsNull(dnMG);
if (nc==0) then SRFIMG = vMG[dnMG];
end
CreateHistogram(SRFIMG,0); CreatePyramid(SRFIMG,0);
CopySubobjects(RL,SRFIMG,"GEOREF");
smlContrast.InputUpperLimit = 3500;
smlContrast.Compute(SRFIMG,clut1$,clut2$,clut3$);
CloseRaster(MG); CloseRaster(SRFIMG);
end
# ------------------------------------------------------------
# IF atcor = 3, PRODUCE VALUES FOR PVI & PBI RASTERS:
# Refer to B36
if (atcor == 3) then begin
printf("\n\nPRODUCING PVI & PBI RASTERS:");
angrad = -atan(bslineslope);
sinang = sin(angrad); cosang = cos(angrad);
for each SRFIRL begin
srfiRL = SRFIRL;
if (srfiRL > 0) then begin
if (pNA) then srfiNR = SRFINA;
if (pNB and pNA==0) then srfiNR = SRFINB;
tsrfiNR = srfiNR - srfiNRint;
pbif = srfiRL * cosang - tsrfiNR * sinang;
pbi = round(pfac * pbif)+ pbioff;
pvif = srfiRL * sinang + tsrfiNR * cosang;
pvi = round(pfac * pvif) + pvioff;
if (pbi < 1) then pbi = 1;
if (pvi < 1) then pvi = 1;
if (pbi > maxpvipbi) then pbi = maxpvipbi;
if (pvi > maxpvipbi) then pvi = maxpvipbi;
PBI = pbi; PVI = pvi;
end
end
CreateHistogram(PBI,0); CreateHistogram(PVI,0);
CreatePyramid(PBI,0); CreatePyramid(PVI,0);
CopySubobjects(RL,PVI,"GEOREF");
CopySubobjects(RL,PBI,"GEOREF");
smlContrast.Power = 1;
smlContrast.InputLowerLimit = 800;
smlContrast.InputUpperLimit = 2000;
smlContrast.Compute(PVI,clut1$,clut2$,clut3$);
smlContrast.InputLowerLimit = 1;
smlContrast.InputUpperLimit = 1500;
smlContrast.Compute(PBI,clut1$,clut2$,clut3$);
CloseRaster(PBI); CloseRaster(PVI);
end
if (pNA) then CloseRaster(SRFINA);
if (pNB) then CloseRaster(SRFINB);
CloseRaster(SRFIRL);
# ------------------------------------------------------------
# CLEAR CONSOLE WINDOW, WRITE TITLE & WRITE PROCESSING-REPORT
# TO CONSOLE WINDOW: Refer to B37.
clear(); writeTitle();
printf("SRFI.sml PROCESSING-REPORT:\n\n");
printf(" SITE NAME: %s\n",site$);
printf(" COLLECTION DATE: %8d\n",cdate);
printf(" DAY OF YEAR: %d\n",doy);
if (pdate>0) then begin
printf(" PROCESSING DATE: %8d\n",pdate);
end
printf("SUN ELEVATION ANGLE: %5.2f deg\n",sunelevang);
printf(" IMAGER: %s\n",imager$);
printf(" IMAGE DN SCALE: %d to %d\n",dnlow,dnhigh);
if (imager == 4) then begin
printf(" GAIN-CODE: %s\n",gc$);
p$ = "LPGS (EOS Data Gateway)";
if (sc == 2) then p$ = "NLAPS (EarthExplorer)";
printf(" DATA VENDOR: %s\n",p$);
end
printf(" CORRECTION LEVEL: ");
p4$ = "Make No Atmospheric Corrections";
p5$ = "Correct for Atmospheric Reflectance Only";
p6$ = "Correct for All Atmospheric Effects";
if (atcor == 1) then printf("%s\n",p4$);
if (atcor == 2) then printf("%s\n",p5$);
if (atcor == 3) then printf("%s\n",p6$);
printf(" delcf: %4.2f Percentage",delcf);
printf(" Points\n");
printf(" msfac: %6.4f\n",msfac);
printf(" icRL: %6.4f\n\n",icRL);
if (imager == 4) then begin
printf(" GAIN-CODE: %s\n",gc$);
p$ = "LPGS (EOS Data Gateway)";
if (sc == 2) then p$ = "NLAPS (EarthExplorer)";
printf(" DATA VENDOR: %s\n",p$);
end
printf(" BAND dnb dnpath srfpath a");
printf(" c m\n");
if (pCB) then begin
mCB = aCB * cCB;
printf(" CB %4d %6d %7.3f",dnbCB,dnpathCB,srfpathCB);
printf(" %8.6f %5.3f %9.7f\n",aCB,cCB,mCB);
end
if (pBL) then begin
mBL = aBL * cBL;
printf(" BL %4d %6d %7.3f",dnbBL,dnpathBL,srfpathBL);
printf(" %8.6f %5.3f %9.7f\n",aBL,cBL,mBL);
end
mGL = aGL * cGL;
if (pYL) then begin
mYL = aYL * cYL;
printf(" YL %4d %6d %7.3f",dnbYL,dnpathYL,srfpathYL);
printf(" %8.6f %5.3f %9.7f\n",aYL,cYL,mYL);
end
mRL = aRL * cRL;
printf(" GL %4d %6d %7.3f",dnbGL,dnpathGL,srfpathGL);
printf(" %8.6f %5.3f %9.7f\n",aGL,cGL,mGL);
printf(" RL %4d %6d %7.3f",dnbRL,dnpathRL,srfpathRL);
printf(" %8.6f %5.3f %9.7f\n",aRL,cRL,mRL);
if (pRE) then begin
mRE = aRE * cRE;
printf(" RE %4d %6d %7.3f",dnbRE,dnpathRE,srfpathRE);
printf(" %8.6f %5.3f %9.7f\n",aRE,cRE,mRE);
end
if (pNA) then begin
mNA = aNA * cNA;
printf(" NA %4d %6d %7.3f",dnbNA,dnpathNA,srfpathNA);
printf(" %8.6f %5.3f %9.7f\n",aNA,cNA,mNA);
end
if (pNB) then begin
mNB = aNB * cNB;
printf(" NB %4d %6d %7.3f",dnbNB,dnpathNB,srfpathNB);
printf(" %8.6f %5.3f %9.7f\n",aNB,cNB,mNB);
end
if (pMA) then begin
mMA = aMA * cMA;
printf(" MA %4d %6d %7.3f",dnbMA,dnpathMA,srfpathMA);
printf(" %8.6f %5.3f %9.7f\n",aMA,cMA,mMA);
end
if (pMB) then begin
mMB = aMB * cMB;
printf(" MB %4d %6d %7.3f",dnbMB,dnpathMB,srfpathMB);
printf(" %8.6f %5.3f %9.7f\n",aMB,cMB,mMB);
end
if (pMC) then begin
mMC = aMC * cMC;
printf(" MC %4d %6d %7.3f",dnbMC,dnpathMC,srfpathMC);
printf(" %8.6f %5.3f %9.7f\n",aMC,cMC,mMC);
end
if (pMD) then begin
mMD = aMD * cMD;
printf(" MD %4d %6d %7.3f",dnbMD,dnpathMD,srfpathMD);
printf(" %8.6f %5.3f %9.7f\n",aMD,cMD,mMD);
end
if (pME) then begin
mME = aME * cME;
printf(" ME %4d %6d %7.3f",dnbME,dnpathME,srfpathME);
printf(" %8.6f %5.3f %9.7f\n",aME,cME,mME);
end
if (pMF) then begin
mMF = aMF * cMF;
printf(" MF %4d %6d %7.3f",dnbMF,dnpathMF,srfpathMF);
printf(" %8.6f %5.3f %9.7f\n",aMF,cMF,mMF);
end
if (pMG) then begin
mMG = aMG * cMG;
printf(" MG %4d %6d %7.3f",dnbMG,dnpathMG,srfpathMG);
printf(" %8.6f %5.3f %9.7f\n",aMG,cMG,mMG);
end
printf(" SRFpath Model Parameter: %6.4f\n\n",p);
printf(" CONVERSION EQUATIONS for Band XX:\n");
printf(" mXX = aXX * cXX\n");
printf(" srftoaXX = (dnXX - dnbXX) * aXX\n");
printf(" srfpathXX = (dnpathXX - dnbXX) * aXX\n");
printf(" srfapcXX = srftoaXX - srfpathXX\n");
printf(" srfsfcXX = srfapcXX * cXX\n");
printf(" srfsfcXX = (dnXX - dnbXX - dnpathXX)");
printf(" * mXX\n");
printf(" srfitoaXX = round(100 * srftoaXX)\n");
printf(" srfiapcXX = round(100 * srfapcXX)\n");
printf(" srfisfcXX = round(100 * srfsfcXX)\n");
printf(" dnXX = dnpathXX + round");
printf("(0.01 * srfiXX / mXX)\n");
# ------------------------------------------------------------
# GIVE REFERENCE FOR ANALYSIS HINTS: Refer to B37.
printf("\nFOR ANALYSIS HINTS: Refer to B37.\n\n");
printf("TO SAVE THE CONSOLE WINDOW TEXT AS A REPORT:\n");
printf(" 1. RIGHT-CLICK IN THE CONSOLE WINDOW.\n");
printf(" 2. SELECT THE Save As... OPTION.\n");
printf(" 3. NAVIGATE TO THE DESIRED LOCATION.\n");
printf(" 4. PROVIDE A REPORT NAME (or OVERWRITE).\n");
printf(" 5. CLICK OK.");
Download