[AE] AE处理Landsat8

【内容】

• 获取影像头文件中的数据
• 打开指定波段的landsat8数据
• 计算指定波段的真实大气顶部反射率（landsat8）

using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Text;
using System.Text.RegularExpressions;
using ESRI.ArcGIS.DataSourcesRaster;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.SpatialAnalyst;

namespace PM25Application {
    class Landsat8Tools {

        /// <summary>
        ///     获取影像头文件中的数据
        /// </summary>
        /// <param name="imgDir">Landsat8的路径</param>
        /// <param name="regex_str">正则规则</param>
        ///     【举例】"DATE_ACQUIRED = (.*)"
        /// <returns></returns>
        public static string GetStringFromMetadataFile_landsat8 (string imgDir, string regex_str) {
            // 获得头文件内容
            var files = Directory.GetFiles (imgDir, "*MTL.txt");
            if (files.Length > 1) { //出错
                return null;
            }
            string file_str = File.ReadAllText (files[0]);

            // 获得日期
            Match match = null;
            match = Regex.Match (file_str, regex_str);
            if (match.Success == true) {
                return match.Groups[1].Value;
            }
            return null;
        }

        /// <summary>
        ///     打开landsat8的指定波段
        /// </summary>
        /// <param name="landsat_dir"></param>
        /// <param name="band"></param>
        /// <returns></returns>
        public static IGeoDataset OpenRasterByBand_landsat8 (string landsat_dir, int band) {
            var files = Directory.GetFiles (landsat_dir, "*B" + band + ".TIF");
            if (files.Length > 1) {
                return null;
            }
            string fp = files[0]; //文件路径
            return Tools.OpenRaster (fp);
        }

        /// <summary>
        ///     计算指定波段的真实大气顶部反射率（landsat8）
        /// </summary>
        ///     Debug OK
        /// <param name="landsat_dir">landsat8数据所在文件夹</param>
        /// <param name="band">指定波段</param>
        /// <returns>计算结果IGeoDataset</returns>
        public static IGeoDataset CalTOAReality_landsat8 (string landsat_dir, int band) {
            string REFLECTANCE_MULT_BAND_rule = string.Concat ("REFLECTANCE_MULT_BAND_", band, " = (.*)");
            string SUN_ELEVATION_rule = "SUN_ELEVATION = (.*)";

            // 从头文件中读取REFLECTANCE_MULT_BAND
            string REFLECTANCE_MULT_BAND_str = Landsat8Tools.GetStringFromMetadataFile_landsat8 (landsat_dir, REFLECTANCE_MULT_BAND_rule);
            Decimal REFLECTANCE_MULT_BAND;
            Decimal.TryParse (REFLECTANCE_MULT_BAND_str, System.Globalization.NumberStyles.Any, null, out REFLECTANCE_MULT_BAND);
            // 从头文件中读取SUM_ELEVATION
            string SUN_ELEVATION_str = Landsat8Tools.GetStringFromMetadataFile_landsat8 (landsat_dir, SUN_ELEVATION_rule);
            double SUN_ELEVATION = System.Convert.ToDouble (SUN_ELEVATION_str);
            SUN_ELEVATION = SUN_ELEVATION * Math.PI / 180; //转换为弧度
            // 读取指定波段
            IGeoDataset pGD_band = Landsat8Tools.OpenRasterByBand_landsat8 (landsat_dir, band);
            // 栅格计算器：计算真实TOAz
            IMapAlgebraOp pMapAlgebraOp = new RasterMapAlgebraOpClass ();
            pMapAlgebraOp.BindRaster (pGD_band, "band");
            string formula = string.Concat ("( [band] * ", REFLECTANCE_MULT_BAND, " - 0.1 ) / Sin( ", SUN_ELEVATION, " )");
            formula = string.Concat ("Con( [band] == 0, 0, ", formula, " )"); //去掉背景值
            IGeoDataset tmp = pMapAlgebraOp.Execute (formula);
            // 小于等于0 --> 0 ; 超出1的为1
            IMapAlgebraOp pMapAlgebraOp2 = new RasterMapAlgebraOpClass ();
            pMapAlgebraOp2.BindRaster (tmp, "band");
            string formula2 = "Con( [band] <= 0 , 0 , Con( [band] >= 1 , 1, [band] ) )";
            IGeoDataset pGD_TOAr = pMapAlgebraOp2.Execute (formula2);
            return pGD_TOAr;
        }
    }
}