C#实现纵横断面求解
测绘技能大赛纵横断面,发现自身还有很多不足,以至于最后不得不添加多余的代码来进行实现,不过好在完成了求解。
point.cs
public class point
{
public string name;
public double x;
public double y;
public double z;
public double d;
}
Calculate.cs
public class Calculate
{
public point K0 = new point();
public point K1 = new point();
public point K2 = new point();
public point H0 = new point();
public point A = new point();
public point B = new point();
public List<point> data = new List<point>();
//存储内插点最临近点
public List<point> pointss;
public Calculate()
{
}
public Calculate(point testA,point testB,point H,List<point> readdata)
{
H0 = H;
A = testA;
B = testB;
for (int i = 0; i < readdata.Count; i++)
{
if (readdata[i].name == "K0")
{
K0 = readdata[i];
}
else if (readdata[i].name == "K1")
{
K1 = readdata[i];
}
else if (readdata[i].name == "K2")
{
K2 = readdata[i];
}
if(readdata[i].name!="K0"&& readdata[i].name != "K1"&& readdata[i].name != "K2")
{
data.Add(readdata[i]);
}
}
}
/// <summary>
/// 计算方位角,返回弧度
/// </summary>
/// <param name="p1"></param>
/// <param name="p2"></param>
/// <returns></returns>
public double azimuthangle(point p1,point p2)
{
double y = p1.y - p2.y;
double x = p1.x - p2.x;
double angle = Math.Atan2(y,x);
if(x>=0&&y>=0)
{
angle = angle;
}
if(x<=0&&y>=0)
{
angle = angle;
}
if(x <= 0 && y <= 0)
{
angle = Math.PI + angle + Math.PI;
}
if(x>=0&&y<=0)
{
angle = Math.PI + angle + Math.PI;
}
return angle;
}
public double azimuthangle(double y1,double x1)
{
double y = y1;
double x = x1;
double angle = Math.Atan2(y, x);
if (x >= 0 && y >= 0)
{
angle = angle;
}
if (x <= 0 && y >= 0)
{
angle = angle;
}
if (x <= 0 && y <= 0)
{
angle = Math.PI + angle + Math.PI;
}
if (x >= 0 && y <= 0)
{
angle = Math.PI + angle + Math.PI;
}
return angle;
}
/// <summary>
/// 由弧度换算为角度
/// </summary>
/// <param name="hudu"></param>
/// <returns></returns>
public double huzhuandu(double hudu)
{
double angle = 180 * (hudu / Math.PI);
return angle;
}
/// <summary>
/// 由角度返回度分秒
/// </summary>
/// <param name="angle"></param>
/// <returns></returns>
public string duzhuan(double angle)
{
string du;
string fen;
string miao;
double ldu = angle * 3600;
int i = 0;
int j = 0;
int k = 0;
while(ldu>=3600)
{
ldu = ldu - 3600;
i++;
}
while(ldu>=60)
{
ldu = ldu - 60;
j++;
}
if(j<10)
{
fen = "0" + j.ToString();
}
else
{
fen = j.ToString();
}
if(ldu<10)
{
miao = "0" + ldu.ToString("00");
}
else
{
miao = ldu.ToString("00");
}
du = i.ToString();
string line = du + "°" + fen + "'" + miao + "\"";
return line;
}
/// <summary>
/// 计算内插点高程
/// </summary>
/// <param name="points"></param>
/// <param name="number"></param>
/// <returns></returns>
public double neica(point points,int number)
{
pointss = new List<point>();
double h = 0;
int num = 0;
List<point> p = new List<point>();
List<double> d = new List<double>();
for(int i=0;i<data.Count;i++)
{
if(points==data[i])
{
num = 1;
return data[i].z;
}
p.Add(data[i]);
d.Add(Math.Sqrt((points.x - data[i].x) * (points.x - data[i].x) + (points.y - data[i].y) * (points.y - data[i].y)));
p[i].d = d[i];
}
double k = 0;
point pp = new point();
for(int i=0;i<d.Count-1;i++)
{
for(int j=0;j<d.Count-1-i;j++)
{
if(d[j]>d[j+1])
{
k = d[j];
d[j] = d[j + 1];
d[j + 1] = k;
pp = p[j];
p[j] = p[j + 1];
p[j + 1] = pp;
}
}
}
double he = 0;
double ld = 0;
if (num == 0)
{
for (int i = 0; i < number; i++)
{
pointss.Add(p[i]);
he = he + (pointss[i].z / pointss[i].d);
ld = ld + 1 / pointss[i].d;
}
}
h = he / ld;
return h;
}
/// <summary>
/// 梯形面积
/// </summary>
/// <param name="h1"></param>
/// <param name="h2"></param>
/// <param name="p1"></param>
/// <param name="p2"></param>
/// <param name="h0"></param>
/// <returns></returns>
public double trapezoidsize(double h1,double h2,point p1,point p2,double h0)
{
double s = 0;
double l = Math.Sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
s = (h1 + h2 - 2 * h0) * l / 2;
return s;
}
/// <summary>
/// 长度
/// </summary>
/// <param name="p1"></param>
/// <param name="p2"></param>
/// <returns></returns>
public double longe(point p1, point p2)
{
double l = Math.Sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
return l;
}
/// <summary>
/// 内插点计算
/// </summary>
/// <param name="p2"></param>
/// <param name="p1"></param>
/// <param name="l"></param>
/// <returns></returns>
public List<point> neichapoint(point p2,point p1,double l)
{
double xmin = 0;
double xmax = 0;
if(p1.x>=p2.x)
{
xmin = p2.x;
xmax = p1.x;
}
if(p1.x <= p2.x)
{
xmax = p2.x;
xmin = p1.x;
}
List<point> neichapoints = new List<point>();
//point point = new point();
double angle = azimuthangle( p2, p1);
double x = p1.x;
double y = p1.y;
do
{
point point = new point();
x = x + l * Math.Cos(angle);
y = y + l * Math.Sin(angle);
point.x = x;
point.y = y;
point.z = neica(point, 5);
neichapoints.Add(point);
} while (x >= xmin && x <= xmax);
neichapoints.Remove(neichapoints[neichapoints.Count - 1]);
return neichapoints;
}
public List<point> neichapoint(point p2, point p1, double xmin,double xmax,double l)
{
List<point> neichapoints = new List<point>();
//point point = new point();
double angle = azimuthangle(p2, p1);
double x = p1.x;
double y = p1.y;
do
{
point point = new point();
x = x + l * Math.Cos(angle);
y = y + l * Math.Sin(angle);
point.x = x;
point.y = y;
point.z = neica(point, 5);
neichapoints.Add(point);
} while (x >= xmin && x <= xmax);
neichapoints.Remove(neichapoints[neichapoints.Count - 1]);
return neichapoints;
}
}
一些公共变量
List<string> showdata = new List<string>();
List<string> savedata = new List<string>();
List<point> points = new List<point>();
point H0 = new point();
point testA = new point();
point testB = new point();
读取数据
string FileName = "";
OpenFileDialog openFileDialog = new OpenFileDialog();
openFileDialog.InitialDirectory = Environment.GetFolderPath(Environment.SpecialFolder.Personal);
openFileDialog.Filter = "文本文件(*.txt)|*.txt|所有文件(*.*)|*.*";
if (openFileDialog.ShowDialog(this) == DialogResult.OK)
{
FileName = openFileDialog.FileName;
}
if (FileName != "")
{
try
{
using (StreamReader re = new StreamReader(FileName))
{
string linshidata =null;
while((linshidata=re.ReadLine())!=null)
{
showdata.Add(linshidata);
}
}
}
catch
{
}
}
else
{
showdata.Add("请读取数据!");
}
写入数据
H0.z= double.Parse(showdata[0].Split(',')[1]);
testA.x= double.Parse(showdata[2].Split(',')[1]);
testA.y = double.Parse(showdata[2].Split(',')[2]);
testB.x = double.Parse(showdata[3].Split(',')[1]);
testB.y = double.Parse(showdata[3].Split(',')[2]);
for (int i=5;i<showdata.Count;i++)
{
point p = new point();
p.name = showdata[i].Split(',')[0];
p.x = double.Parse(showdata[i].Split(',')[1]);
p.y = double.Parse(showdata[i].Split(',')[2]);
p.z = double.Parse(showdata[i].Split(',')[3]);
points.Add(p);
}
查看读入数据
textBox1.Clear();
for(int i=0;i<showdata.Count;i++)
{
textBox1.AppendText(showdata[i]+"\r\n");
}
计算
readdata();
Calculate c = new Calculate(testA, testB, H0, points);
double hk1 = c.neica(c.K1, 5);
savedata.Add("K1内插高程:"+hk1.ToString()+"m");
for (int i = 0; i < c.pointss.Count; i++)
{
savedata.Add("点号:"+c.pointss[i].name+" 距离:"+c.pointss[i].d.ToString()+"m");
}
double zduanmianlong = c.longe(c.K0, c.K1) + c.longe(c.K1, c.K2);
savedata.Add("纵断面长度:" + zduanmianlong.ToString()+"m");
List<point> k0k1 = c.neichapoint(c.K1, c.K0, 10);
savedata.Add(" ");
savedata.Add("K0--K1内插点");
for(int i=0;i<k0k1.Count;i++)
{
savedata.Add("p" + i.ToString() + ": " + "(" + k0k1[i].x.ToString() +"," + k0k1[i].y.ToString()+","+k0k1[i].z.ToString()+")");
}
savedata.Add(" ");
savedata.Add("K1--K2内插点");
List<point> k1k2 = c.neichapoint(c.K2,c.K1,10);
for (int i = 0; i < k0k1.Count; i++)
{
savedata.Add("p" + (i+k0k1.Count).ToString() + ": " + "(" + k1k2[i].x.ToString() + "," + k1k2[i].y.ToString() + "," + k1k2[i].z.ToString() + ")");
}
savedata.Add(" ");
double zdum1 = 0;
double zdum2 = 0;
for(int i=0;i<k0k1.Count;i++)
{
if(i==0)
{
zdum1 = zdum1 +Math.Abs( c.trapezoidsize(c.K0.z,k0k1[i].z,c.K0,k0k1[i],c.H0.z));
}
if(i==k0k1.Count-1)
{
zdum1 = zdum1 + Math.Abs(c.trapezoidsize(c.K1.z,k0k1[i].z,c.K1,k0k1[i],c.H0.z));
}
if (i > 0)
{
zdum1 = zdum1 + Math.Abs(c.trapezoidsize(k0k1[i].z, k0k1[i - 1].z, k0k1[i], k0k1[i - 1], c.H0.z));
}
}
for (int i = 0; i < k1k2.Count; i++)
{
if (i == 0)
{
zdum2 = zdum2 + Math.Abs(c.trapezoidsize(c.K1.z, k1k2[i].z, c.K1, k1k2[i], c.H0.z));
}
if (i == k1k2.Count - 1)
{
zdum2 = zdum2 + Math.Abs(c.trapezoidsize(c.K2.z, k1k2[i].z, c.K2, k1k2[i], c.H0.z));
}
if (i > 0)
{
zdum2 = zdum2 + Math.Abs(c.trapezoidsize(k1k2[i].z, k1k2[i - 1].z, k1k2[i], k1k2[i - 1], c.H0.z));
}
}
savedata.Add("纵断面面积:" + (zdum1 + zdum2).ToString() + "㎡");
savedata.Add(" ");
point pointm1 = new point();
pointm1.x = (c.K0.x + c.K1.x) / 2;
pointm1.y = (c.K0.y + c.K1.y) / 2;
pointm1.z = c.neica(pointm1,5);
point pointm2 = new point();
pointm2.x = (c.K1.x + c.K2.x) / 2;
pointm2.y = (c.K1.y + c.K2.y) / 2;
pointm2.z = c.neica(pointm2, 5);
double anglem1 = c.azimuthangle(c.K0,c.K1);
double anglem2 = c.azimuthangle(c.K1,c.K2);
anglem1 = anglem1 + 0.5 * Math.PI;
anglem2 = anglem2 + 0.5 * Math.PI;
if(anglem1>=2*Math.PI)
{
anglem1 -= 2 * Math.PI;
}
if (anglem2 >= 2 * Math.PI)
{
anglem2 -= 2 * Math.PI;
}
double zx1 = Math.Abs(25 * Math.Cos(anglem1));
double zy1 = Math.Abs(25 * Math.Sin(anglem1));
double x1max = pointm1.x + zx1;
double x1min = pointm1.x - zx1;
double y1max = pointm1.y + zy1;
double y1min = pointm1.y - zy1;
double zx2 = Math.Abs(25 * Math.Cos(anglem2));
double zy2 = Math.Abs(25 * Math.Sin(anglem2));
double x2max = pointm2.x + zx2;
double x2min = pointm2.x - zx2;
double y2max = pointm2.y + zy2;
double y2min = pointm2.y - zy2;
point pointm11 = new point();
point pointm12 = new point();
point pointm21 = new point();
point pointm22 = new point();
if (c.azimuthangle(y1max-y1min,x1max-x1min)==anglem1|| c.azimuthangle(y1max - y1min, x1max - x1min) == anglem1-Math.PI|| c.azimuthangle(y1max - y1min, x1max - x1min) == anglem1 + Math.PI)
{
pointm11.x = x1min;
pointm11.y = y1min;
pointm11.z = c.neica(pointm11,5);
pointm12.x = x1max;
pointm12.y = y1max;
pointm12.z = c.neica(pointm12, 5);
}
else
{
pointm11.x = x1min;
pointm11.y = y1max;
pointm11.z = c.neica(pointm11, 5);
pointm12.x = x1max;
pointm12.y = y1min;
pointm12.z = c.neica(pointm12, 5);
}
if (c.azimuthangle(y2max - y2min, x2max - x2min) == anglem2 || c.azimuthangle(y2max - y2min, x2max - x2min) == anglem2 - Math.PI || c.azimuthangle(y2max - y2min, x2max - x2min) == anglem2 + Math.PI)
{
pointm21.x = x2min;
pointm21.y = y2min;
pointm21.z = c.neica(pointm21, 5);
pointm22.x = x2max;
pointm22.y = y2max;
pointm22.z = c.neica(pointm22, 5);
}
else
{
pointm21.x = x2min;
pointm21.y = y2max;
pointm21.z = c.neica(pointm21, 5);
pointm22.x = x2max;
pointm22.y = y2min;
pointm22.z = c.neica(pointm22, 5);
}
savedata.Add("以m1为中点的断面");
List<point> m1m11 = new List<point>();
m1m11 = c.neichapoint(pointm11,pointm1,5);
for(int i=0;i<m1m11.Count;i++)
{
savedata.Add("V"+i.ToString()+":("+m1m11[i].x.ToString()+","+m1m11[i].y.ToString()+","+m1m11[i].z.ToString()+")");
}
savedata.Add("V" + m1m11.Count.ToString() + ":(" + pointm11.x.ToString() + "," + pointm11.y.ToString() + "," + pointm11.z.ToString() + ")");
savedata.Add("m1:"+"("+pointm1.x.ToString()+","+pointm1.y.ToString()+","+pointm1.z.ToString()+")");
List<point> m1m12 = new List<point>();
m1m12 = c.neichapoint(pointm12,pointm1,5);
for (int i = 0; i < m1m12.Count; i++)
{
savedata.Add("V" + (i+m1m11.Count+1).ToString() + ":(" + m1m12[i].x.ToString() + "," + m1m12[i].y.ToString() + "," + m1m12[i].z.ToString() + ")");
}
savedata.Add("V" + (m1m11.Count+m1m12.Count+1).ToString() + ":(" + pointm12.x.ToString() + "," + pointm12.y.ToString() + "," + pointm12.z.ToString() + ")");
savedata.Add(" ");
savedata.Add("以m2为中点的断面");
List<point> m2m21 = new List<point>();
m2m21 = c.neichapoint(pointm21, pointm2, 5);
for (int i = 0; i < m2m21.Count; i++)
{
savedata.Add("V" + i.ToString() + ":(" + m2m21[i].x.ToString() + "," + m2m21[i].y.ToString() + "," + m2m21[i].z.ToString() + ")");
}
savedata.Add("V" + m2m21.Count.ToString() + ":(" + pointm21.x.ToString() + "," + pointm21.y.ToString() + "," + pointm21.z.ToString() + ")");
savedata.Add("m2:" + "(" + pointm2.x.ToString() + "," + pointm2.y.ToString() + "," + pointm2.z.ToString() + ")");
List<point> m2m22 = new List<point>();
m2m22 = c.neichapoint(pointm22, pointm2, 5);
for (int i = 0; i < m2m22.Count; i++)
{
savedata.Add("V" + (i + m2m21.Count + 1).ToString() + ":(" + m2m22[i].x.ToString() + "," + m2m22[i].y.ToString() + "," + m2m22[i].z.ToString() + ")");
}
savedata.Add("V" + (m2m21.Count + m2m22.Count + 1).ToString() + ":(" + pointm22.x.ToString() + "," + pointm22.y.ToString() + "," + pointm22.z.ToString() + ")");
查看计算数据
textBox1.Clear();
for (int i = 0; i < savedata.Count; i++)
{
textBox1.AppendText(savedata[i] + "\r\n");
}
保存数据
string FileName = "";
SaveFileDialog saveFileDialog = new SaveFileDialog();
saveFileDialog.InitialDirectory = Environment.GetFolderPath(Environment.SpecialFolder.Personal);
saveFileDialog.Filter = "文本文件(*.txt)|*.txt|所有文件(*.*)|*.*";
if (saveFileDialog.ShowDialog(this) == DialogResult.OK)
{
FileName = saveFileDialog.FileName;
}
using (StreamWriter wr = new StreamWriter(FileName))
{
for(int i = 0; i < savedata.Count; i++)
{
wr.WriteLine(savedata[i]) ;
}
}
计算结果
点名:P22 x:4536.141 y:3378.766 z:19.502
点名:P28 x:4540.773 y:3379.044 z:18.344
点名:P15 x:4531.567 y:3373.371 z:15.047
点名:P21 x:4540.86 y:3376.395 z:19.223
点名:P27 x:4537.846 y:3387.734 z:18.239
18.4873777206298
x:4565.25107325863 y:3363.12142744759 z:14.0293771596633
x:4556.49014651726 y:3367.94285489518 z:16.789932402761
x:4547.72921977589 y:3372.76428234276 z:17.9498121098473
x:4538.96829303452 y:3377.58570979035 z:19.1022355540576
x:4530.35974055222 y:3382.66386829545 z:18.3802918961682
x:4521.93091862894 y:3388.04484971018 z:18.0538178533896
x:4513.50209670566 y:3393.42583112491 z:17.6654429305891
x:4505.07327478238 y:3398.80681253964 z:13.1100276519559
S:572.858834593616
m0内插点:
x:4551.70878627621 y:3364.86703662932 z:16.8209559303007
x:4549.29807255241 y:3360.48657325863 z:16.7632507590104
x:4546.88735882862 y:3356.10610988795 z:16.7180116121305
x:4544.47664510482 y:3351.72564651726 z:16.7035868692548
x:4542.06593138103 y:3347.34518314658 z:16.6987038662613
x:4556.53021372379 y:3373.62796337069 z:15.6964099271893
x:4558.94092744759 y:3378.00842674137 z:14.2643342279404
x:4561.35164117138 y:3382.38889011206 z:12.4282423222864
x:4563.76235489518 y:3386.76935348274 z:12.2093315447363
x:4566.17306861897 y:3391.14981685342 z:11.4731185127121
s:262.547882660782
m1内插点:
x:4513.34500929264 y:3387.59408903836 z:17.0048210784474
x:4510.65451858527 y:3383.37967807672 z:15.0362018484324
x:4507.96402787791 y:3379.16526711508 z:15.9138590746037
x:4505.27353717054 y:3374.95085615344 z:15.0315032363331
x:4502.58304646318 y:3370.7364451918 z:15.0112641043021
x:4518.72599070736 y:3396.02291096164 z:17.809606600276
x:4521.41648141473 y:3400.23732192328 z:15.7392689535267
x:4524.10697212209 y:3404.45173288492 z:15.0217313498982
x:4526.79746282946 y:3408.66614384656 z:14.6126623098665
x:4529.48795353682 y:3412.8805548082 z:14.4700374598023
s:293.682598217233
总结
最后写崩了,有些内容可以写为函数,有些函数需要改动,主要原因在于对测量学等科目不熟悉,在基本概念上出现了问题。