.gen地圖文件的投影編程實現(以墨卡託投影和蘭伯特投影爲例)

      好幾天沒更博了,yogurt最近忙得飛起啊,沒辦法,相信付出老是會有收穫的!每次更博的時候都是yogurt最開心的時候,啦啦啦~~~好了,廢話很少說,趕忙更完寫做業 去了~~~~(>_<)~~~~面試

      今天yogurt要給你們分享的是我前幾周剛學會的地圖投影!就是把.gen的矢量形式的地圖文件讀入,並經過數學運算,實現墨卡託投影和蘭伯特投影的效果~~(能夠利用ArcGIS的DtaInteroperability的快速導入工具進行查看投影后效果哦~~開心到飛起,感受本身完成了什麼了不得的大事情呢,哈哈哈)小程序

=================================yogurt小課堂開課了=============================數組

    首先yogurt給你們簡單普及一下座標系這個概念!要知道咱們有三種座標!三維的地心座標系(XYZ)、三維橢球體的地理座標系(經緯度)以及二維平面上的投影座標系。它們之間的關係就如圖:因此咱們不能直接將一張地圖(通過投影了)直接轉換到另外一種座標系下的地圖,而是須要複雜的變化操做。先從投影座標系變回地理座標系,再由地理座標系變到大地座標系,最後經過三參數法或者七參數法進行座標變換,而後換算到另外一種橢球體對應的地理座標系中,最後進行投影獲得另外一個座標系下的地圖。ide

 

    所以,不一樣的投影座標系的投影方法對應的數學運算是不同的。然而yogurt這個小程序實現了從經緯度利用墨卡託投影(正軸等角圓柱投影)和蘭伯特投影(等角圓錐投影)投影到二維平面單位轉化爲米或者公里。工具

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

好了,話很少說,Yogur上代碼了。注意:我用來存儲.gen數據的二維數組我簡單解釋一下:第一維存放線號(第0號位置存放總的線數目,第i號位置存放i),第二維存放點號(第0號位置存放該線對應的總的點數目,第j號存放該線上第j個點的座標),如:a[7][0].x:表明第7條線對應的點的數目;a[3][4]:表明第3條線的第4個點的座標。code

 

  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;            //第1、二標準緯線
 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
View Code

這樣就大功告成啦!我把結構的導入ArcGIS裏面試一試哈,見證奇蹟的時候到了~~blog

原文件:數學

投影后:string

相關文章
相關標籤/搜索