pktools 2.6.7
Processing Kernel for geospatial data
Egcs.cc
1/**********************************************************************
2Egcs.cc: Conversions from and to european grid coding system
3Copyright (C) 2008-2012 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along with pktools. If not, see <http://www.gnu.org/licenses/>.
19***********************************************************************/
20#include "Egcs.h"
21#include <iostream>
22#include <sstream>
23#include <iomanip>
24#include <assert.h>
25
26Egcs::Egcs(){
27}
28
29// Egcs::Egcs(unsigned short level)
30// : m_level(level){
31// }
32
33
34Egcs::Egcs(unsigned short level)
35{
36 m_level=level;
37}
38
39unsigned short Egcs::cell2level(const std::string& cellCode) const{
40 size_t pos=cellCode.find("-");
41 std::string TILE=cellCode.substr(0,pos);
42 unsigned short level=0;
43 int base_level=19-(TILE.size()/2-1)*3;
44 int quad_level=0;
45 if(pos!=std::string::npos)
46 quad_level=cellCode.size()-pos-1;
47 if(quad_level>1)
48 level=base_level-2;
49 else if(quad_level>0)
50 level=base_level-1;
51 else
52 level=base_level;
53 return level;
54}
55
56Egcs::~Egcs(){
57}
58
59unsigned short Egcs::res2level(double resolution) const{
60 double base=pow(10,log(resolution*4.0)/log(10.0));
61 double diff=base/(2*resolution);
62 return 0.5+(log(base)/log(10.0)*3+1-diff);
63}
64
65double Egcs::getResolution() const{
66 unsigned short exponent=(m_level+1)/3;
67 double base=pow(10.0,exponent);
68 if((m_level)%3==2)
69 return(base/4);
70 else if((m_level)%3==0)
71 return(base/2);
72 else
73 return(base);
74}
75
76void Egcs::force2grid(double& ulx, double& uly, double& lrx, double &lry) const{
77 double dx=getResolution();
78 double dy=dx;
79 //ulx
80 ulx=floor(ulx);
81 ulx-=static_cast<int>(ulx)%(static_cast<int>(dx));
82 //uly
83 uly=ceil(uly);
84 if(static_cast<int>(uly)%static_cast<int>(dy))
85 uly+=dy;
86 uly-=static_cast<int>(uly)%(static_cast<int>(dy));
87 //lrx
88 lrx=ceil(lrx);
89 if(static_cast<int>(lrx)%static_cast<int>(dx))
90 lrx+=dx;
91 lrx-=static_cast<int>(lrx)%(static_cast<int>(dx));
92 //lry
93 lry=floor(lry);
94 lry-=static_cast<int>(lry)%(static_cast<int>(dy));
95}
96
97void Egcs::cell2bb(const std::string& cellCode, int &ulx, int &uly, int &lrx, int &lry) const{
98 size_t pos=cellCode.find("-");
99 std::string TILE=cellCode.substr(0,pos);
100 std::string TILEX=TILE.substr(0,TILE.size()/2);
101 std::string TILEY=TILE.substr(TILE.size()/2);
102 std::string QUAD;
103 char QUAD1,QUAD2;
104 std::istringstream stilex(TILEX);
105 std::istringstream stiley(TILEY);
106 int llx,lly;
107 stilex >> llx;
108 stiley >> lly;
109 llx*=getBaseSize();
110 lly*=getBaseSize();
111 switch((19-m_level)%3){
112 case(0)://there should be no QUAD
113 assert(pos==std::string::npos);
114 break;
115 case(2)://there is a QUAD2
116 assert(pos+1!=std::string::npos);
117 QUAD=cellCode.substr(pos+1);
118 QUAD2=QUAD.substr(1,1).c_str()[0];
119 switch(QUAD2){
120 case('A'):
121 break;
122 case('C'):
123 lly+=getBaseSize()/4;
124 break;
125 case('D'):
126 lly+=getBaseSize()/4;
127 case('B'):
128 llx+=getBaseSize()/4;
129 break;
130 }
131 case(1)://QUAD1: deliberate fall through from case(2)!
132 if(!QUAD.size()){
133 assert(pos+1!=std::string::npos);
134 QUAD=cellCode.substr(pos+1);
135 }
136 QUAD1=QUAD.substr(0,1).c_str()[0];
137 switch(QUAD1){
138 case('A'):
139 break;
140 case('C'):
141 lly+=getBaseSize()/2;
142 break;
143 case('D'):
144 lly+=getBaseSize()/2;
145 case('B'):
146 llx+=getBaseSize()/2;
147 break;
148 }
149 break;
150 }
151 ulx=llx;
152 uly=static_cast<int>(lly+getSize());
153 lrx=static_cast<int>(llx+getSize());
154 lry=lly;
155}
156
157void Egcs::cell2mid(const std::string& cellCode, double& midX, double& midY) const{
158 int ulx,uly,lrx,lry;
159 cell2bb(cellCode,ulx,uly,lrx,lry);
160 midX=(ulx+lrx)/2.0;
161 midY=(lry+uly)/2.0;
162}
163
164std::string Egcs::geo2cell(double geoX, double geoY) const{
165 int ndgts=7-(m_level+1)/3;
166 double xcel=static_cast<int>(geoX)/getBaseSize();
167 double ycel=static_cast<int>(geoY)/getBaseSize();
168 std::ostringstream osx;
169 std::ostringstream osy, osxy;
170 // osx << setprecision(ndgts) << xcel;
171 // osy << setprecision(ndgts) << ycel;
172 // osx << setprecision(0) << fixed << geoX;
173 // osy << setprecision(0) << fixed << geoY;
174 osx << std::fixed << geoX;
175 osy << std::fixed << geoY;
176 std::string quad1="";
177 std::string quad2="";
178 switch((19-m_level)%3){
179 case(2):
180 if(static_cast<int>(geoX)%(getBaseSize()/2)>=getBaseSize()/2.0){
181 if(static_cast<int>(geoY)%(getBaseSize()/2)>=getBaseSize()/2.0)
182 quad2="D";
183 else
184 quad2="B";
185 }
186 else if(static_cast<int>(geoY)%getBaseSize()>=getBaseSize()/4.0)
187 quad2="C";
188 else
189 quad2="A";
190 case(1)://deliberate fall through!
191 if(static_cast<int>(geoX)%getBaseSize()>=getBaseSize()/2.0){
192 if(static_cast<int>(geoY)%getBaseSize()>=getBaseSize()/2.0)
193 quad1="-D";
194 else
195 quad1="-B";
196 }
197 else if(static_cast<int>(geoY)%getBaseSize()>=getBaseSize()/2.0)
198 quad1="-C";
199 else
200 quad1="-A";
201 break;
202 default:
203 break;
204 }
205 osxy << osx.str().substr(0,ndgts) << osy.str().substr(0,ndgts) << quad1 << quad2;
206 return osxy.str();
207}
208