[AE] ArcGIS Engine处理Landsat8数据

秒速五厘米 2021-12-11 19:49 416阅读 0赞

【内容】

*  获取影像头文件中的数据
 *  打开指定波段的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;
            }
        }
    }