? ????好幾天沒(méi)更博了，yogurt最近忙得飛起啊，沒(méi)辦法，相信付出總是會(huì)有收獲的！每次更博的時(shí)候都是yogurt最開(kāi)心的時(shí)候，啦啦啦~~~好了，廢話不多說(shuō)，趕緊更完寫(xiě)作業(yè) 去了~~~~(>_<)~~~~

??????今天yogurt要給大家分享的是我前幾周剛學(xué)會(huì)的地圖投影！就是把.gen的矢量形式的地圖文件讀入，并通過(guò)數(shù)學(xué)運(yùn)算，實(shí)現(xiàn)墨卡托投影和蘭伯特投影的效果~~（可以利用ArcGIS的DtaInteroperability的快速導(dǎo)入工具進(jìn)行查看投影后效果哦~~開(kāi)心到飛起，感覺(jué)自己完成了什么了不起的大事情呢，哈哈哈）

=================================yogurt小課堂開(kāi)課了=============================

????首先yogurt給大家簡(jiǎn)單普及一下坐標(biāo)系這個(gè)概念！要知道我們有三種坐標(biāo)！三維的地心坐標(biāo)系（XYZ）、三維橢球體的地理坐標(biāo)系（經(jīng)緯度）以及二維平面上的投影坐標(biāo)系。它們之間的關(guān)系就如圖：所以我們不能直接將一張地圖（經(jīng)過(guò)投影了）直接轉(zhuǎn)換到另一種坐標(biāo)系下的地圖，而是需要復(fù)雜的變化操作。先從投影坐標(biāo)系變回地理坐標(biāo)系，再由地理坐標(biāo)系變到大地坐標(biāo)系，最后通過(guò)三參數(shù)法或者七參數(shù)法進(jìn)行坐標(biāo)變換，然后換算到另一種橢球體對(duì)應(yīng)的地理坐標(biāo)系中，最后進(jìn)行投影得到另一個(gè)坐標(biāo)系下的地圖。

?

????因此，不同的投影坐標(biāo)系的投影方法對(duì)應(yīng)的數(shù)學(xué)運(yùn)算是不一樣的。然而yogurt這個(gè)小程序?qū)崿F(xiàn)了從經(jīng)緯度利用墨卡托投影（正軸等角圓柱投影）和蘭伯特投影（等角圓錐投影）投影到二維平面單位轉(zhuǎn)化為米或者千米。

=================================下課了================================

好了，話不多說(shuō)，Yogur上代碼了。注意：我用來(lái)存儲(chǔ).gen數(shù)據(jù)的二維數(shù)組我簡(jiǎn)單解釋一下：第一維存放線號(hào)（第0號(hào)位置存放總的線數(shù)目，第i號(hào)位置存放i），第二維存放點(diǎn)號(hào)（第0號(hào)位置存放該線對(duì)應(yīng)的總的點(diǎn)數(shù)目，第j號(hào)存放該線上第j個(gè)點(diǎn)的坐標(biāo)），如：a[7][0].x:代表第7條線對(duì)應(yīng)的點(diǎn)的數(shù)目；a[3][4]:代表第3條線的第4個(gè)點(diǎn)的坐標(biāo)。

?

1 1 // 投影.cpp : 定義控制臺(tái)應(yīng)用程序的入口點(diǎn)。 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];//用于記錄每條線有的點(diǎn)數(shù) 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;//記錄線數(shù) 26 26 int j = 0;//記錄點(diǎn)數(shù) 27 27 fscanf(fp, "%s", s); //s不是線編號(hào)就是文件讀取結(jié)束的標(biāo)志“END” 28 28 while (s[0] != 'E') 29 29 { 30 30 i++; 31 31 j = 0; 32 32 fscanf(fp, "%s", s); //s不是點(diǎn)對(duì)就是線讀取結(jié)束的標(biāo)志“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; //記錄線的數(shù)目 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++) //為每條線開(kāi)辟足夠的點(diǎn)空間 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位置賦值（線數(shù)和點(diǎn)數(shù)） 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];//第一個(gè)位置存放線數(shù)或者點(diǎn)數(shù) 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);//過(guò)濾掉線號(hào) 64 64 for (int j = 1; j <= a[i][0].x; j++)//第j個(gè)點(diǎn) 65 65 fscanf(fp, "%lf,%lf", &a[i][j].x, &a[i][j].y); 66 66 fscanf(fp, "%s", s);//過(guò)濾掉END標(biāo)志 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個(gè)點(diǎn) 81 81 { 82 82 double lon = p[i][j].x*PI / 180; 83 83 double lat = p[i][j].y*PI / 180; //被投影的點(diǎn)的經(jīng)緯度 84 84 85 85 double a = 6378245; //長(zhǎng)半軸 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; //原始經(jīng)緯度 90 90 double SP1 = 20 * PI / 180; 91 91 double SP2 = 40 * PI / 180; //第一、二標(biāo)準(zhǔn)緯線 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個(gè)點(diǎn) 127 127 { 128 128 double lon = p[i][j].x*PI / 180; 129 129 double lat = p[i][j].y*PI / 180; //被投影的點(diǎn)的經(jīng)緯度 130 130 131 131 double a = 6378245; //長(zhǎng)半軸 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; //原始經(jīng)度 136 136 double B0 = 42 * PI / 180; //標(biāo)準(zhǔn)緯度 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 View Code

這樣就大功告成啦！我把結(jié)構(gòu)的導(dǎo)入ArcGIS里面試一試哈，見(jiàn)證奇跡的時(shí)候到了~~

原文件：

投影后：

?

轉(zhuǎn)載于:https://www.cnblogs.com/to-sunshine/p/6048438.html

總結(jié)

以上是生活随笔為你收集整理的.gen地图文件的投影编程实现（以墨卡托投影和兰伯特投影为例）的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。