AE大气扩散模型算法.ppt

基于Visual C#和ArcEngine的二次开发,GIS空间分析基于GIS的大气环境影响预测与模拟,大气污染模型高斯大气扩散模式 单源大气污染扩散模拟,1 大气扩散模式 大气扩散模式是对大气污染物散布过程的模拟。按污染源的几何特征可分为电源、线源、面源扩散模式；按污染源排放特性分有正常工况模式和非正常排放模式；按假设条件不同分有烟流扩散模式和烟团扩散模式；按模式适用范围不同分有短距离扩散模式和中长距离输送模式等。具体操作时应根据当地地形、气候等环境要素以及污染源亲何况正确选取并适当修正，不能盲目照搬。 主要的大气扩散模式有： （1）高架电源扩散模式（静风模式、有效源高等） （2）面源扩散模式（虚电源模式、窄烟云模式：在多面源情形下计算效果较好） （3）线源扩散模式（电源求和法、风向与线源平行时 扩散模式）,一、大气污染模型,2 高斯大气扩散模式 高斯大气扩散模式属于高架电源的烟羽大气扩散模式，它是计算释入大气中的气载污染物下风向浓度的应用最广的方法。高斯扩散模式所描述的扩散过程主要有以下几种：下垫面平坦、 开阔、 性质均匀 ,平均流场平直、 稳定 ，不考虑风场的切变；扩散过程中 ,污染物本身是被动、 保守的 ,即污染物和空气无相对运动 ,且扩散过程中污染物无损失、 无转化 ,污染物在地面被反射；扩散在同一温度层集中发生 ,平均风速 1. 0 m/ s ;适用范围一般 1020 km。在均匀、定常的湍流大气中污染物浓度满足正态分布，但实际大气一般不满足均匀、定常条件，因此高斯模式应用于下垫面均匀平坦、气流稳定的小尺度扩散问题更为有效。,3高斯计算公式 对于在恒定气象条件（指风向、风速、大气稳定度不随时间而变）下的高架点源的连续排放，在考虑了烟羽在地面的全反射后，下风向任一点的污染物浓度C(x,y,z)可由下式计算： 式中：C(x,y,z)：下风向某点(x,y,z)处的空气中污染物的浓度（mg/m3）； x：下风向距离（m）； y：横截风向距离（m）； z：距地面的高度（m）； Q：气载污染物源强，即释放率（mg/s）； ：排气筒出口处的平均风速（m/s）； 、 ：分别为水平方向和垂直方向扩散参数（m），它是下风向距离x及大气稳定度的函数； He：有效排放高度（m），其中He=H+ ； H：排气筒距地面的几何高度（m）； ：烟气抬升高度（m）。,表1 太阳辐射等级数,表2 大气稳定度等级,2高架源回归系数（Py、qy、Pz、qz）确定 国际原子能机构（IAEA）根据相关实验数据，推荐了一组适用于大粗糙度（Z01m）和高架源的弥散系数公式，如表3，据此求得相应地扩散参数。,表3 与排放高度有关的扩散参数（ 、 ）,、 是下风距离X及大气稳定度的函数，而下风方向的距离X是一变量，它是风速和时间的函数。高架电源烟气扩散是三维扩散过程，x方向上风速较大，以对流扩散为主，y方向和z方向则以弥散扩散为主。本实例编程只是实现了二维扩散模型。即只是实现了xy方向上的扩散模拟。,二、单源大气污染扩散模拟,1 相关菜单的添加 在空间分析主菜单中添加如下图所示菜单项。,2 大气污染扩散分析窗口的设计 设计一个大气污染扩散分析窗口。,1,2,3,4,5,6,7,8,9,10,11,12,这些Label的AutoSize属性设置为：False BorderStyle属性设置为：Fixed3D,一些主要控件的属性设置如下：,其中大气污染扩散分析窗口的AutoSizeMode属性设置为：GrowAndShrink；MaximizeBox属性设置为：False,3菜单与窗口连接的实现 双击“单源大气污染”菜单项，在响应的Click事件中添加如下代码： private void 单源大气污染ToolStripMenuItem_Click(object sender, EventArgs e) 大气污染扩散分析 form = new 大气污染扩散分析(); form.Show(); ,4 高斯大气扩散功能的模拟实现 步骤1：重写“大气污染扩散分析”类的构造函数， 鉴于在该模型实现过程中多次需要调用主窗口中的函数或控件，所以为构造函数添加参数：主窗口 MainForm，同时在单源大气污染ToolStripMenuItem_Click事件中修改代码。 大气污染扩散分析 form = new 大气污染扩散分析(); 修改为： 大气污染扩散分析 form = new 大气污染扩散分析(this);,步骤2：在“大气污染扩散分析”窗口类中对相关参数进行初始化。 using ESRI.ArcGIS.Geometry; /相关接口所在的命名空间 using ESRI.ArcGIS.Carto; /窗口类中定义的相关全局变量 public DiffusedParameter diffusedParameter50; /定义大气扩散参数结构体，该结构体是自己编写的，这里放在初始化函数后 public DiffusedParameter diffusedParameter100; public DiffusedParameter diffusedParameter180; const double PI = 3.1415926; private 主窗口 parentForm; private int radius = -1; /辐射等级 private int atmosphereStability = -1; /初始化大气稳定度等级 private double Py, Pz, Qy, Qz; /计算中使用的扩散参数 private double Q = 500; /气载污染物源强，即释放率（mg/s） IEnvelope envelope;,public 大气污染扩散分析(主窗口 MainForm) InitializeComponent(); /初始化“风向”、“风速”、“释放高度”、“计算高度” windspeed.Text = “5“; height.Text = “100“; interpolateHeight.Text = “100“; ILayer lyr; IFeatureLayer featureLyr; parentForm = MainForm; envelope = new EnvelopeClass();,/初始化图层组合框 for (int i = 0; i MainForm.GetLayerCount(); i+) /GetLayerCount()是在主窗口中编写的获取图层数函数。 lyr = MainForm.GetLayer(i); /GetLayer()是在主窗口编写的根据索引获取图层函数 try featureLyr = (IFeatureLayer)lyr; catch (Exception e) featureLyr = null; /当图层为点图层时，即将其对应的图层名添加到 comboBox_layer组合框中 if (featureLyr != null ,/初始化扩散参数表，此时的0、1、2、3、4、5对应大气稳定度等级A、B、C、D、E、F diffusedParameter50 = new DiffusedParameter6; diffusedParameter500.Py = 1.503; diffusedParameter500.Qy = 0.833; diffusedParameter500.Pz = 0.151; diffusedParameter500.Qz = 1.219; diffusedParameter501.Py = 0.876; diffusedParameter501.Qy = 0.823; diffusedParameter501.Pz = 0.127; diffusedParameter501.Qz = 1.108; diffusedParameter502.Py = 0.659; diffusedParameter502.Qy = 0.807; diffusedParameter502.Pz = 0.165; diffusedParameter502.Qz = 0.996;,diffusedParameter503.Py = 0.640; diffusedParameter503.Qy = 0.784; diffusedParameter503.Pz = 0.215; diffusedParameter503.Qz = 0.885; diffusedParameter504.Py = 0.801; diffusedParameter504.Qy = 0.754; diffusedParameter504.Pz = 0.264; diffusedParameter504.Qz = 0.774; diffusedParameter505.Py = 1.294; diffusedParameter505.Qy = 0.718; diffusedParameter505.Pz = 0.241; diffusedParameter505.Qz = 0.662;,diffusedParameter100 = new DiffusedParameter6; diffusedParameter1000.Py = 0.170; diffusedParameter1000.Qy = 1.296; diffusedParameter1000.Pz = 0.051; diffusedParameter1000.Qz = 1.317; diffusedParameter1001.Py = 0.324; diffusedParameter1001.Qy = 1.025; diffusedParameter1001.Pz = 0.070; diffusedParameter1001.Qz = 1.151; diffusedParameter1002.Py = 0.466; diffusedParameter1002.Qy = 0.866; diffusedParameter1002.Pz = 0.137; diffusedParameter1002.Qz = 0.985;,diffusedParameter1003.Py = 0.504; diffusedParameter1003.Qy = 0.818; diffusedParameter1003.Pz = 0.265; diffusedParameter1003.Qz = 0.818; diffusedParameter1004.Py = 0.411; diffusedParameter1004.Qy = 0.882; diffusedParameter1004.Pz = 0.487; diffusedParameter1004.Qz = 0.652; diffusedParameter1005.Py = 0.253; diffusedParameter1005.Qy = 1.057; diffusedParameter1005.Pz = 0.717; diffusedParameter1005.Qz = 0.486;,diffusedParameter180 = new DiffusedParameter6; diffusedParameter1800.Py = 0.671; diffusedParameter1800.Qy = 0.903; diffusedParameter1800.Pz = 0.0245; diffusedParameter1800.Qz = 0.500; diffusedParameter1801.Py = 0.415; diffusedParameter1801.Qy = 0.903; diffusedParameter1801.Pz = 0.0330; diffusedParameter1801.Qz = 1.320; diffusedParameter1802.Py = 0.232; diffusedParameter1802.Qy = 0.903; diffusedParameter1802.Pz = 0.104; diffusedParameter1802.Qz = 0.997;,diffusedParameter1803.Py = 0.208; diffusedParameter1803.Qy = 0.903; diffusedParameter1803.Pz = 0.307; diffusedParameter1803.Qz = 0.734; diffusedParameter1804.Py = 0.345; diffusedParameter1804.Qy = 0.903; diffusedParameter1804.Pz = 0.546; diffusedParameter1804.Qz = 0.557; diffusedParameter1805.Py = 0.671; diffusedParameter1805.Qy = 0.903; diffusedParameter1805.Pz = 0.484; diffusedParameter1805.Qz = 0.500; ,public struct DiffusedParameter /自定义一个大气扩散参数结构体 public double Py; public double Pz; public double Qy; public double Qz; /在主窗口编写的函数 public int GetLayerCount() return this.axMapControl1.LayerCount; public ILayer GetLayer(int index) return this.axMapControl1.get_Layer(index); ,运行程序，显示如下：,comboBox_layer中添加的点图层,步骤3：响应comboBox_layer的SelectIndexChanged函数，根据comboBox中的图层文本获取相应的特征图层，并将图层中的“名称”字段值添加到comboBox_Polluter组合框中。 private void comboBox_layer_SelectedIndexChanged(object sender, EventArgs e) /重置大气污染源集合 comboBox_polluter.Items.Clear(); comboBox_polluter.Text = “; string lyrName = boBox_layer.Text; ILayer lyr = parentForm.GetLayerByname(lyrName); if (null = lyr) MessageBox.Show(“没有找到“ + lyrName + “图层！“); return; ITable table = (ITable)lyr; IQueryFilter queryFilter = new QueryFilterClass();,ICursor cursor = table.Search(queryFilter, true); int fieldIndex = table.FindField(“名称“); if (fieldIndex = 0) /当图层属性表中存在“名称”属性列时，即将该列的所有属性字段值添加到combBox_polluter组合框中 IRow row; while (row = cursor.NextRow() != null) comboBox_polluter.Items.Add(row.get_Value(fieldIndex); else MessageBox.Show(“所选择的数据层不能提供污染源信息“); /初始化label_address、label_name、label_person、label_type label_address.Text = “; label_name.Text = “; label_person.Text = “; label_type.Text = “; ,运行程序，显示如下：,步骤4：响应comboBox_polluter的SelectedIndexChanged函数，当选择comboBox_polluter中任意一数据项时，在label_name、 label_type、 label_address、 label_person中添加从数据库中获取的污染源的名称、所属行业名称、地址和负责人信息。 using System.Data.OleDb; private OleDbConnection conn; private void comboBox_polluter_SelectedIndexChanged(object sender, EventArgs e) string name = boBox_polluter.Text; string sql = “SELECT 代码, 名称, 8#行业类别_名称, 县(区、市、旗), 乡(镇), 街(村)、门牌号, 单位负责人 FROM (SELECT 代码, 名称, 8#行业类别_名称, 县(区、市、旗), 乡(镇), 街(村)、门牌号, 单位负责人 FROM G101 UNION ALL SELECT 代码, 名称, 8#行业类别_名称, 县(区、市、旗), 乡(镇), 街(村)、门牌号, 单位负责人 FROM G201) AS derivedtbl_1 WHERE (名称 = “ + name + “)“;,conn = new OleDbConnection(“Provider=Microsoft.Jet.OLEDB.4.0;Data Source=“ + Application.StartupPath + “Database.mdb“); OleDbCommand cmd = new OleDbCommand(sql, conn); /申明的OleDbCommand对象，表示要对数据源执行的SQL语句 conn.Open(); OleDbDataReader reader = cmd.ExecuteReader(); int cols = reader.FieldCount; object values = new objectcols; if (reader.Read() reader.GetValues(values); label_name.Text = values1.ToString(); label_type.Text = values2.ToString(); label_address.Text = values3.ToString() + values4.ToString() + values5.ToString(); label_person.Text = values6.ToString(); reader.Close(); conn.Close(); ,运行程序，显示如下：,步骤5：选中大气污染扩散分析窗口中所有的RadioButton按钮，在radioButton1中响应其Click事件，即可以实现对所有RadioButton控件的Click事件的响应。 private void radioButton1_Click(object sender, EventArgs e) RadioButton radioButton = (RadioButton)sender; radius = Convert.ToInt32(radioButton.Text); /将radioButton的Text值赋给radius变量 ,步骤6：响应btnOK的Click事件，实现单源大气污染源的高斯扩散模拟。 using ESRI.ArcGIS.SpatialAnalyst; private void btnOK_Click(object sender, EventArgs e) double speed; try speed = Convert.ToDouble(this.windspeed.Text); catch (Exception e1) MessageBox.Show(“风速填写错误，请查看！“); return; if (speed 0) MessageBox.Show(“风速填写错误，请查看！“); return; ,switch (radius) /根据地面风速和太阳辐射等级计算大气稳定度 case -2: if (speed 2.9 ,case 0: atmosphereStability = 4; break; case 1: if (speed 1.9 ,case 2: if (speed 1.9 ,/设定扩散参数 if (height.Text = “) MessageBox.Show(“请设置释放高度参数！“); return; if (height.Text.Trim() = “50“) Py = diffusedParameter50atmosphereStability - 1.Py; Pz = diffusedParameter50atmosphereStability - 1.Pz; Qy = diffusedParameter50atmosphereStability - 1.Qy; Qz = diffusedParameter50atmosphereStability - 1.Qz; if (height.Text.Trim() = “100“) Py = diffusedParameter100atmosphereStability - 1.Py; Pz = diffusedParameter100atmosphereStability - 1.Pz; Qy = diffusedParameter100atmosphereStability - 1.Qy; Qz = diffusedParameter100atmosphereStability - 1.Qz; ,if (height.Text.Trim() = “180“) Py = diffusedParameter180atmosphereStability - 1.Py; Pz = diffusedParameter180atmosphereStability - 1.Pz; Qy = diffusedParameter180atmosphereStability - 1.Qy; Qz = diffusedParameter180atmosphereStability - 1.Qz; /获取污染源 ILayer pLayer = parentForm.GetLayerByname(comboBox_layer.Text); if (null = pLayer) MessageBox.Show(“没有找到“ + comboBox_layer.Text + “图层！“); return; if (label_name.Text = “) MessageBox.Show(“请选择污染源！“); return; ,ITable ptable = (ITable)pLayer; IQueryFilter queryFilter = new QueryFilterClass(); queryFilter.WhereClause = “名称=“ + label_name.Text + “; ICursor cursor = ptable.Search(queryFilter, true); IRow prow = cursor.NextRow(); int index = ptable.FindField(“shape“); IPoint wuranyuan = (IPoint)prow.get_Value(index); /设置边界范围 envelope.XMax = wuranyuan.X + 5000; envelope.XMin = wuranyuan.X - 5000; envelope.YMax = wuranyuan.Y + 5000; envelope.YMin = wuranyuan.Y - 5000; /采样点图层 IEnvelope pEnvelope = new EnvelopeClass(); IPoint pPoint = new ESRI.ArcGIS.Geometry.PointClass(); IPolygon pPolygon = GetInterpolateArea(); /根据风向确定污染影响的可能范围，该函数的原理实现在后面有说明,IPointCollection pPointCollection = (IPointCollection)pPolygon; /根据确定的多边形生成采样点集，IPointCollection是一个点集接口，IPolygon、IPointCollection接口都位于Geometry命名空间中，通过两者之间的强制转换，可以根据Polygon生成相应的点集对象。 PointF pPointF = new PointFpPointCollection.PointCount; for (int num = 0; num pPointF.Length; num+) pPointFnum.X = (float)(IPoint)pPointCollection.get_Point(num).X; pPointFnum.Y = (float)(IPoint)pPointCollection.get_Point(num).Y; Polygon ClipPolygon = new Polygon(pPointF);,pPolygon.QueryEnvelope(pEnvelope); /将pPolygon的QueryEnelope赋给pEnvelop /为即将生成的采样点矢量图层，添加属性字段。包括一个“FID“字段、一个“c“字段 IFields PropertyFields = new FieldsClass(); IFieldsEdit FieldsEdit = (IFieldsEdit)PropertyFields; IField Field = new FieldClass(); IFieldEdit FieldEdit = (IFieldEdit)Field; FieldEdit.Name_2 = “FID“; FieldEdit.AliasName_2 = “FID“; FieldEdit.Type_2 = esriFieldType.esriFieldTypeOID; FieldsEdit.AddField(FieldEdit); Field = new FieldClass(); FieldEdit = (IFieldEdit)Field; FieldEdit.Name_2 = “c“; FieldEdit.AliasName_2 = “c“; FieldEdit.Type_2 = esriFieldType.esriFieldTypeDouble; FieldsEdit.AddField(FieldEdit);,/生成cayangdian矢量图层 IFeatureLayer lyr = parentForm.CreateFeatureLayerInmemeory(“caiyangdian“, “caiyangdian“, parentForm.getMapControl().Map.SpatialReference, esriGeometryType.esriGeometryPoint, PropertyFields); / CreateFeatureLayerInmemeory函数在主窗口中 ITable CYDTable = (ITable)lyr; IRow CYDRow; int shapeIndex = CYDTable.FindField(“SHAPE“); /“SHAPE“字段是在CreateFeatureLayerInmemory函数中生成的 int nongduIndex = CYDTable.FindField(“c“); /为采样点附浓度，添加采样点 int direction = Convert.ToInt32(winddirection.Value); /风向 Random rnd = new Random(); PointF tempPt = new PointF(); double interpolatehigh; try interpolatehigh = Convert.ToDouble(interpolateHeight.Text); catch (Exception e2) MessageBox.Show(“计算高度填写有误，请检查！“); return; ,for (double x = pEnvelope.XMin; x = pEnvelope.XMax; x = x + 200) for (double y = pEnvelope.YMin; y = pEnvelope.YMax; y = y + 200) tempPt.X = (float)x; tempPt.Y = (float)y; if (ClipPolygon.PtInPolygon(tempPt) != -1) / PtInPolygon函数在Polygon自定义类中用于判断某一点是在Polygon内Polygon边界上还是在Polygon外 pPoint.X = x; pPoint.Y = y; IPoint ptemPt = new PointClass(); ptemPt.X = pPoint.X; ptemPt.Y = pPoint.Y;,/IPoint接口中的ConstrainAngle(double constrainAngle,IPoint anchor,bool allowOpposite)这个方法的运算结果如下所述: 首先用anchor和constrainAngle在一个以anchor为极点的极坐标系中构造一个无限的射线(矢量)然后从IPoint当前的坐标出发作一个垂直于这条射线的直线，交点就是IPoint的新的当前值.如果allowOpposite为true，那么还可以向射线的反向延长线作垂线。这个方法可以用来约束新点必须在某条直线上。如图所示,pPoint.ConstrainAngle(PI * direction) / (float)180, wuranyuan, true); double c_x = System.Math.Sqrt(Math.Pow(Math.Abs(wuranyuan.X - pPoint.X), 2) + Math.Pow(Math.Abs(wuranyuan.Y - pPoint.Y), 2); double c_y = System.Math.Sqrt(Math.Pow(Math.Abs(pPoint.X - ptemPt.X), 2) + Math.Pow(Math.Abs(pPoint.Y - ptemPt.Y), 2); double c = Calculate(c_x, c_y, interpolatehigh); /Calculate为高斯模型计算 if (c 0.0000000001) /过滤掉c浓度值小于0.0000000001的采样点，再将c 0.0000000001的采样点的值存储下来 CYDRow = CYDTable.CreateRow(); CYDRow.set_Value(shapeIndex, ptemPt); CYDRow.set_Value(nongduIndex, c); CYDRow.Store(); ,IRasterLayer IdwRasterLyr = parentForm.IDW(lyr, “c“, 6, 20, 30, envelope); /利用反距离加权插值法根据浓度”c”值进行插值，IDW函数在主窗口中定义 IExtractionOp2 extract = new RasterExtractionOpClass(); /切割出插值后的栅格图，因为原本插值后的图形覆盖整个envelope区域，由于风向原因，这里用pPolygon进行切割，即只显示主导风向的多边形区域 IGeoDataset geoDataset = extract.Polygon(IGeoDataset)IdwRasterLyr.Raster, pPolygon, true); IRaster pOutRast = geoDataset as IRaster; IRasterLayer pOutLayer = new RasterLayerClass(); pOutLayer = parentForm.SetRsLayerClassifiedColor(pOutRast, 10); / SetRsLayerClassifiedColor函数在主窗口中定义，用于实现栅格图形的分级设色，分级设色后有一大部分区域显示为白色 ILayerEffects lyrEffects = (ILayerEffects)pOutLayer; lyrEffects.Transparency = 70; /设置透明度 parentForm.getMapControl().AddLayer(lyr); parentForm.getMapControl().AddLayer(pOutLayer); parentForm.getMapControl().Extent = envelope; ,/高斯模型计算 public double Calculate(double x, double y, double z) double y, z; y = this.Py * Math.Pow(x, Qy); z = this.Pz * Math.Pow(x, Qz); double speed = Convert.ToDouble(windspeed.Text); /风速 double He = Convert.ToDouble(this.height.Text) + 10; /计算有效排放高度 double F = Math.Exp(-(z - He) * (z - He) / (2 * (z * z) + Math.Exp(-(z + He) * (z + He) / (2 * (z * z); return (this.Q * Math.Exp(-y * y / (2 * (y * y) * F) / (2 * PI * y * z * speed); ,/根据污染源和风向计算可能被污染的区域。这里自定义了一个算法函数。 其中记pt为污染源所在点；pt1和pt2为与风向垂直方向垂线与矩形边框（evelope）的交点；分别用1、2、3、4来记录上、右、下、左边；用pt1Flag、pt2Flag来记录第一各交点和第二个交点在哪条边上；用pts0、pts1、pts2、pts3来记录右上、右下、左下、左上四个顶点。 （1）计算夹角angle的斜率，用tan表示。,（2）确定交点所在边 记 tempY = envelope.YMax; / 即矩形边界的Y的最大值 tempX = (tempY - pt.Y) / tan + pt.X; 若 tempX = envelope.XMin 若 tempY = envelope.YMin & tempY = envelope.YMax且此时pt1Flag = 0 则 此时该交点为第一个交点且标识该交点在右边上，否则该交点为第二个交点，且标识第二个交点在第二条边上 依此，分别确定两个交点所在的边。,（3）用顺时针、逆时针来区别污染影响的可能区域及其相反区域。为后续存点做好准备。 我们定义与风向相反向的，两个交点中从第2个点到第1个点的方向为污染影响区域所在方位。即图中影响的多边形区域即为顺时针方向。 （4）根据确定的顺时针、逆时针方向来确定多边形区域的顶点，并将其存储下来。 如上图中，方向为顺时针，则除了存储两个交点pt1、pt2外，还要存储pts0、pts1两个顶点。因为可能存在凸凹多边形，所以可能需要存储的顶点有多种可能，如下图，此时即需要存储三个顶点。 从pt2到pt1则为顺时针方向，此时除存储pt1和pt2外，还要存储pts3、pts0、pts1。具体的实现见后面代码。,public IPolygon GetInterpolateArea() ILayer lyr = this.parentForm.GetLayerByname(comboBox_layer.Text); if (null = lyr) MessageBox.Show(“没有找到“ + comboBox_layer.Text + “图层“); return null; ITable table = (ITable)lyr; IQueryFilter queryFilter = new QueryFilterClass(); queryFilter.WhereClause = “名称=“ + label_name.Text + “; ICursor cursor = table.Search(queryFilter, true); IRow row = cursor.NextRow(); int index = table.FindField(“shape“);,IPoint pt = (IPoint)row.get_Value(index); int direction = Convert.ToInt32(winddirection.Value); /风向 double angle = (direction + 90) % 180; /与风向垂直方向的夹角 double tan = Math.Tan(PI * angle / 180); /与风向垂直的斜率 IPoint pt1 = new PointClass(); int pt1Flag = 0; /记录第一个交点在那条边上 1234分别表示上右下左边 IPoint pt2 = new PointClass(); int pt2Flag = 0; /记录第二个交点在那条边上 double tempX, tempY;,/依次判断 ,上边 tempY = envelope.YMax; tempX = (tempY - pt.Y) / tan + pt.X; if (tempX = envelope.XMin ,/右边 tempX = envelope.XMax; tempY = tan * (tempX - pt.X) + pt.Y; if (tempY = envelope.YMin /若该交点为第二个交点，则标识第二个交点在第2条边上 ,/下边 tempY = envelope.YMin; tempX = (tempY - pt.Y) / tan + pt.X; if (tempX = envelope.XMin ,/左边 tempX = envelope.XMin; tempY = tan * (tempX - pt.X) + pt.Y; if (tempY = envelope.YMin ,bool IsClockWise = true; /多边形方向 switch (pt1Flag) case 1: if (Math.Abs(direction - 180) 90) /判断顺时针、逆时针都是判断的与风向相反向的两个交点中从第2个点到第1个点的方向 IsClockWise = false; else IsClockWise = true; break; case 2: if (Math.Abs(direction - 90) 90) IsClockWise = false; else IsClockWise = true; break;,case 3: if (direction 90) /此时第一个点在第3条边上，则第二个点在第4条边上 IsClockWise = false; else IsClockWise = true; break; ,/从右上角的顺时针顺序开始存储 PointClass pts = new PointClass4; pts0 = new PointClass(); pts1 = new PointClass(); pts2 = new PointClass(); pts3 = new PointClass(); pts0.X = envelope.XMax; pts0.Y = envelope.YMax; pts1.X = envelope.XMax; pts1.Y = envelope.YMin; pts2.X = envelope.XMin; pts2.Y = envel