好几天没更博了,yogurt最近忙得飞起啊,没办法,相信付出总是会有收获的!每次更博的时候都是yogurt最开心的时候,啦啦啦~~~好了,废话不多说,赶紧更完写作业 去了~~~~(>_<)~~~~
今天yogurt要给大家分享的是我前几周刚学会的地图投影!就是把.gen的矢量形式的地图文件读入,并通过数学运算,实现墨卡托投影和兰伯特投影的效果~~(可以利用ArcGIS的DtaInteroperability的快速导入工具进行查看投影后效果哦~~开心到飞起,感觉自己完成了什么了不起的大事情呢,哈哈哈)
首先yogurt给大家简单普及一下坐标系这个概念!要知道我们有三种坐标!三维的地心坐标系(XYZ)、三维椭球体的地理坐标系(经纬度)以及二维平面上的投影坐标系。它们之间的关系就如图:所以我们不能直接将一张地图(经过投影了)直接转换到另一种坐标系下的地图,而是需要复杂的变化操作。先从投影坐标系变回地理坐标系,再由地理坐标系变到大地坐标系,最后通过三参数法或者七参数法进行坐标变换,然后换算到另一种椭球体对应的地理坐标系中,最后进行投影得到另一个坐标系下的地图。
因此,不同的投影坐标系的投影方法对应的数学运算是不一样的。然而yogurt这个小程序实现了从经纬度利用墨卡托投影(正轴等角圆柱投影)和兰伯特投影(等角圆锥投影)投影到二维平面单位转化为米或者千米。
=================================下课了================================
好了,话不多说,Yogur上代码了。注意:我用来存储.gen数据的二维数组我简单解释一下:第一维存放线号(第0号位置存放总的线数目,第i号位置存放i),第二维存放点号(第0号位置存放该线对应的总的点数目,第j号存放该线上第j个点的坐标),如:a[7][0].x:代表第7条线对应的点的数目;a[3][4]:代表第3条线的第4个点的坐标。
1 1 // 投影.cpp : 定义控制台应用程序的入口点。
2 2 //
3 3
4 4 #include "stdafx.h"
5 5 #include<stdlib.h>
6 6 #include<string.h>
7 7 #include<math.h>
8 8 #define PI 3.1415926535898
9 9
10 10 typedef struct POINT
11 11 {
12 12 double x, y;
13 13 }point;
14 14
15 15 point ** read(char * InFILEname)
16 16 {
17 17 int lines_p[500];//用于记录每条线有的点数
18 18 FILE * fp = fopen(InFILEname, "r");
19 19 if (fp == NULL)
20 20 {
21 21 printf("Can not open it.");
22 22 return 0;
23 23 }
24 24 char s[100] = { 0 };
25 25 int i = 0;//记录线数
26 26 int j = 0;//记录点数
27 27 fscanf(fp, "%s", s); //s不是线编号就是文件读取结束的标志“END”
28 28 while (s[0] != 'E')
29 29 {
30 30 i++;
31 31 j = 0;
32 32 fscanf(fp, "%s", s); //s不是点对就是线读取结束的标志“END”
33 33 while (s[0] != 'E')
34 34 {
35 35 j++;
36 36 fscanf(fp, "%s", s);
37 37 }
38 38 lines_p[i] = j;
39 39 fscanf(fp, "%s", s);
40 40 }
41 41 lines_p[0] = i; //记录线的数目
42 42 rewind(fp);
43 43
44 44 //分配空间
45 45 point **a = (point **)malloc(sizeof(point *)*(lines_p[0] + 1));
46 46 a[0] = (point *)malloc(sizeof(point));
47 47
48 48 for (int i = 1; i <= lines_p[0]; i++) //为每条线开辟足够的点空间
49 49 {
50 50
51 51 a[i] = (point *)malloc((lines_p[i] + 1)* sizeof(point));
52 52
53 53 }
54 54
55 55 //0位置赋值(线数和点数)
56 56 for (int i = 0; i <= lines_p[0]; i++)
57 57 {
58 58 a[i][0].x = a[i][0].y = lines_p[i];//第一个位置存放线数或者点数
59 59 }
60 60
61 61 for (int i = 1; i <= a[0][0].x; i++)//第i条线
62 62 {
63 63 fscanf(fp, "%s", s);//过滤掉线号
64 64 for (int j = 1; j <= a[i][0].x; j++)//第j个点
65 65 fscanf(fp, "%lf,%lf", &a[i][j].x, &a[i][j].y);
66 66 fscanf(fp, "%s", s);//过滤掉END标志
67 67 }
68 68
69 69 return a;
70 70 }
71 71
72 72 void LambertPro(point ** p, char * OutFILEname)
73 73 {
74 74 FILE *fp = fopen(OutFILEname, "w");
75 75 if (fp == NULL)
76 76 printf("Can not open it.");
77 77 for (int i = 1; i <= p[0][0].x; i++) //第i条线
78 78 {
79 79 fprintf(fp, "%d\n", i);
80 80 for (int j = 1; j<=p[i][0].x; j++) //第j个点
81 81 {
82 82 double lon = p[i][j].x*PI / 180;
83 83 double lat = p[i][j].y*PI / 180; //被投影的点的经纬度
84 84
85 85 double a = 6378245; //长半轴
86 86 double b = 6356863.0188; //短半轴
87 87 double e = sqrt(a*a - b*b) / a; //第一偏心率
88 88 double originLon = 105 * PI / 180;
89 89 double originLat = 0; //原始经纬度
90 90 double SP1 = 20 * PI / 180;
91 91 double SP2 = 40 * PI / 180; //第一、二标准纬线
92 92 double Ef = 609600;
93 93 double Nf = 0; //假东和假北
94 94
95 95 double t = tan(PI / 4 - lat / 2) / pow(((1 - e*sin(lat)) / (1 + e*sin(lat))), e / 2);
96 96 double t1 = tan(PI / 4 - SP1 / 2) / pow(((1 - e*sin(SP1)) / (1 + e*sin(SP1))), e / 2);
97 97 double t2 = tan(PI / 4 - SP2 / 2) / pow(((1 - e*sin(SP2)) / (1 + e*sin(SP2))), e / 2);
98 98 double tF = tan(PI / 4 - originLat / 2) / pow(((1 - e*sin(originLat)) / (1 + e*sin(originLat))), e / 2);
99 99 double m1 = cos(SP1) / pow((1 - e*e*sin(SP1)*sin(SP1)), 0.5);
100 100 double m2 = cos(SP2) / pow((1 - e*e*sin(SP2)*sin(SP2)), 0.5);
101 101 double n = (log(m1) - log(m2)) / (log(t1) - log(t2));
102 102 double F = m1 / (n*pow(t1, n));
103 103 double A = n*(lon - originLon);
104 104 double r = a*F*pow(t, n);
105 105 double rF = a*F*pow(tF, n);
106 106
107 107 double E = Ef + r*sin(A);
108 108 double N = Nf + rF - r*cos(A);
109 109 fprintf(fp, "%lf,%lf\n", E, N);
110 110 }
111 111 fprintf(fp, "END\n");
112 112 }
113 113 fprintf(fp, "END\n");
114 114 fclose(fp);
115 115 return;
116 116 }
117 117
118 118 void MercatoPro(point ** p, char * OutFILEname)
119 119 {
120 120 FILE *fp = fopen(OutFILEname, "w");
121 121 if (fp == NULL)
122 122 printf("Can not open it.");
123 123 for (int i = 1; i <= p[0][0].x; i++) //第i条线
124 124 {
125 125 fprintf(fp, "%d\n", i);
126 126 for (int j = 1; j <= p[i][0].x; j++) //第j个点
127 127 {
128 128 double lon = p[i][j].x*PI / 180;
129 129 double lat = p[i][j].y*PI / 180; //被投影的点的经纬度
130 130
131 131 double a = 6378245; //长半轴
132 132 double b = 6356863.0188; //短半轴
133 133 double e = sqrt(a*a - b*b) / a; //第一偏心率
134 134 double e2 = sqrt(a*a - b*b) / b; //第二偏心率
135 135 double L0 = 0; //原始经度
136 136 double B0 = 42 * PI / 180; //标准纬度
137 137
138 138 double K = (a*a / b) / sqrt(1 + e2*e2*cos(B0)*cos(B0)) * cos(B0);
139 139
140 140 double N = K*log(tan(PI / 4 + lat / 2))*(pow((1 - e*sin(lat)) / (1 + e*sin(lat)), e / 2));
141 141 double E = K*(lon - L0);
142 142 fprintf(fp, "%lf,%lf\n", E, N);
143 143 }
144 144 fprintf(fp, "END\n");
145 145 }
146 146 fprintf(fp, "END\n");
147 147 fclose(fp);
148 148 return;
149 149 }
150 150
151 151 int _tmain(int argc, _TCHAR* argv[])
152 152 {
153 153 point ** a = read("CHINA_Arc.gen");
154 154 char outL[15] = "LambertPro.gen";
155 155 char outM[15] = "MercatoPro.gen";
156 156 LambertPro(a, outL);
157 157 MercatoPro(a, outM);
158 158 return 0;
159 159 }
160 View Code
这样就大功告成啦!我把结构的导入ArcGIS里面试一试哈,见证奇迹的时候到了~~
原文件:
投影后: