using System;
using System.Collections.Generic;
using System.Text;

namespace GeoFly
{
    public class PrismInter : SpatialInterpolate
    {
        public GridLayer DEMLayer;
        /// <summary>
        /// 降雨量随高程变化的梯度
        /// </summary>
        public double Grads;

        public GridLayer IDWPrismOut(DateTime date)
        {
            if (this.pStationinfo == null || this.pMeteoData == null)
                return null;

            //栅格数据结构初始化
            m_gridLayer = HydroSimulate.g_GridLayerPara.g_DemLayer.AttributesCopy();

            ProgressBar bar = new ProgressBar();
            bar.Show();
            bar.Text = "正在内插栅格数据";
            if (this.m_CoordType == CoordType.UTM_Coord)
            {
                this.PrepareFromUTM(date);
            }
            else
            {
                this.PrepareFromLongLa(date);
            }
            for (int row = 0; row < m_gridLayer.rowCount; row++)		//开始新一日的空间插值计算
            {
                for (int col = 0; col < m_gridLayer.colCount; col++)
                {
                    double value = -9999;
                    //获取当前计算栅格中心坐标
                    //获取当前计算栅格的坐标
                    LPoint p = m_gridLayer.CellPosition(row, col);
                    double CurrentX = p.X;
                    double CurrentY = p.Y;
                    double sum1 = 0;
                    double sum2 = 0;
                    bool IsOnStation = false;
                    for (int i = 0; i < this.pStationinfo.Count; i++)
                    {
                        double distX = X[i] - p.X;
                        double distY = Y[i] - p.Y;
                        double distSqare = distX * distX + distY * distY;
                        //在站点上
                        if (distSqare < 0.00001)
                        {
                            value = Z[i];
                            IsOnStation = true;
                            break;
                        }
                        int colStation = (int)((X[i] - m_gridLayer.DownLeft_X) / m_gridLayer.resolution) + 1;
                        int rowStation = (int)((Y[i] - m_gridLayer.DownLeft_Y) / m_gridLayer.resolution) + 1;
                        sum1 += (Grads * (DEMLayer[row, col] - DEMLayer[rowStation, colStation]) + Z[i]) / distSqare;
                        sum2 += 1 / distSqare;
                    }
                    //不在站点上
                    if (IsOnStation == false)
                        value = sum1 / sum2;
                    this.m_gridLayer[row, col] = value;
                }
                bar.progressBar1.Value = (int)(row * 100.0 / m_gridLayer.rowCount);
            }
            //图层中间过程输出
            bar.Close();
            return this.m_gridLayer;
        }
        
        /// <summary>
        /// 通用Prism法：这个插值法需要提供插值区域的DEM栅格数据，以及待插值属性随高程变化的梯度值(即回归方程系数，
        ///可以由样本数据拟合)。另外还要提供综合考虑高程、距离、坡度、坡向、植被、下垫面等因子对气象要素的影响的
        ///各样本点综合的权重系数。可见，该插值法需要具体分析才能获得合适的权重系数Pi，具有一定的应用难度。
        /// </summary>
        /// <param name="date"></param>
        /// <returns></returns>
        public GridLayer TruePrismOut(DateTime date,double[] P)
        {
            if (this.pStationinfo == null || this.pMeteoData == null)
                return null;

            //栅格数据结构初始化
            m_gridLayer = HydroSimulate.g_GridLayerPara.g_DemLayer.AttributesCopy();

            ProgressBar bar = new ProgressBar();
            bar.Show();
            bar.Text = "正在内插栅格数据";
            if (this.m_CoordType == CoordType.UTM_Coord)
            {
                this.PrepareFromUTM(date);
            }
            else
            {
                this.PrepareFromLongLa(date);
            }
            for (int row = 0; row < m_gridLayer.rowCount; row++)		//开始新一日的空间插值计算
            {
                for (int col = 0; col < m_gridLayer.colCount; col++)
                {
                    double sum1 = 0;
                    for (int i = 0; i <this.pStationinfo.Count; i++)
                    {
                        int colStation = (int)((X[i] - m_gridLayer.DownLeft_X) / m_gridLayer.resolution) + 1;
                        int rowStation = (int)((Y[i] - m_gridLayer.DownLeft_Y) / m_gridLayer.resolution) + 1;
                        sum1 += P[i] * (Grads * (DEMLayer[row, col] - DEMLayer[rowStation, colStation]) + Z[i]);
                    }
                    this.m_gridLayer[row, col] = sum1;
                }
                bar.progressBar1.Value = (int)(row * 100.0 / m_gridLayer.rowCount);
            }
            //图层中间过程输出
            bar.Close();
            return this.m_gridLayer;
        }
    }
}