1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // Projection_Proj4 //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // crs_distance.cpp //
14 // //
15 // Copyright (C) 2015 by //
16 // Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.org //
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Hamburg //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49
50
51 ///////////////////////////////////////////////////////////
52 // //
53 // //
54 // //
55 ///////////////////////////////////////////////////////////
56
57 //---------------------------------------------------------
58 #include "crs_distance.h"
59
60
61 ///////////////////////////////////////////////////////////
62 // //
63 // //
64 // //
65 ///////////////////////////////////////////////////////////
66
67 //---------------------------------------------------------
CCRS_Distance_Calculator(void)68 CCRS_Distance_Calculator::CCRS_Distance_Calculator(void)
69 {}
70
71 //---------------------------------------------------------
CCRS_Distance_Calculator(const CSG_Projection & Projection,double Epsilon)72 CCRS_Distance_Calculator::CCRS_Distance_Calculator(const CSG_Projection &Projection, double Epsilon)
73 {
74 Create(Projection, Epsilon);
75 }
76
77 //---------------------------------------------------------
~CCRS_Distance_Calculator(void)78 CCRS_Distance_Calculator::~CCRS_Distance_Calculator(void)
79 {}
80
81 //---------------------------------------------------------
Create(const CSG_Projection & Projection,double Epsilon)82 bool CCRS_Distance_Calculator::Create(const CSG_Projection &Projection, double Epsilon)
83 {
84 if( !m_ProjToGCS.Set_Source(Projection)
85 || !m_ProjToGCS.Set_Target(CSG_Projection("+proj=longlat +datum=WGS84", SG_PROJ_FMT_Proj4))
86 || !m_Projector.Set_Target(Projection) )
87 {
88 return( false );
89 }
90
91 m_Epsilon = Epsilon;
92
93 return( true );
94 }
95
96
97 ///////////////////////////////////////////////////////////
98 // //
99 ///////////////////////////////////////////////////////////
100
101 //---------------------------------------------------------
Get_Orthodrome(const TSG_Point & A,const TSG_Point & B,CSG_Shape * pLine)102 double CCRS_Distance_Calculator::Get_Orthodrome(const TSG_Point &A, const TSG_Point &B, CSG_Shape *pLine)
103 {
104 static const TSG_Point P0 = { 0.0, 0.0 };
105
106 TSG_Point P = A;
107
108 if( m_ProjToGCS.Get_Projection(P) )
109 {
110 m_Projector.Set_Source(CSG_Projection(
111 CSG_String::Format("+proj=aeqd +R=6371000 +lon_0=%f +lat_0=%f", P.x, P.y), SG_PROJ_FMT_Proj4)
112 );
113
114 m_Projector.Set_Inverse();
115
116 if( m_Projector.Get_Projection(P = B) )
117 {
118 m_Projector.Set_Inverse(false);
119
120 Add_Segment(P0, P, pLine);
121
122 return( SG_Get_Distance(P0, P) );
123 }
124 }
125
126 return( Get_Distance(A, B) );
127 }
128
129 //---------------------------------------------------------
Get_Loxodrome(const TSG_Point & A,const TSG_Point & B,CSG_Shape * pLine)130 double CCRS_Distance_Calculator::Get_Loxodrome(const TSG_Point &A, const TSG_Point &B, CSG_Shape *pLine)
131 {
132 TSG_Point AA, BB;
133
134 m_Projector.Set_Source(CSG_Projection("+proj=merc +datum=WGS84", SG_PROJ_FMT_Proj4));
135
136 m_Projector.Set_Inverse();
137
138 if( m_Projector.Get_Projection(AA = A)
139 && m_Projector.Get_Projection(BB = B) )
140 {
141 double Length = 0.0;
142
143 m_Projector.Set_Inverse(false);
144
145 Add_Segment(AA, BB, pLine, &Length);
146
147 return( Length );
148 }
149
150 return( Get_Distance(A, B) );
151 }
152
153
154 ///////////////////////////////////////////////////////////
155 // //
156 ///////////////////////////////////////////////////////////
157
158 //---------------------------------------------------------
Get_Distance(TSG_Point A,TSG_Point B)159 double CCRS_Distance_Calculator::Get_Distance(TSG_Point A, TSG_Point B)
160 {
161 if( m_ProjToGCS.Get_Projection(A) && m_ProjToGCS.Get_Projection(B) )
162 {
163 return( SG_Get_Distance_Polar(A, B) );
164 }
165
166 return( 0.0 );
167 }
168
169 //---------------------------------------------------------
Add_Segment(const TSG_Point & A,const TSG_Point & B,CSG_Shape * pLine,double * Length)170 void CCRS_Distance_Calculator::Add_Segment(const TSG_Point &A, const TSG_Point &B, CSG_Shape *pLine, double *Length)
171 {
172 if( SG_Get_Distance(A, B) >= m_Epsilon )
173 {
174 TSG_Point C, CC;
175
176 C.x = CC.x = A.x + 0.5 * (B.x - A.x);
177 C.y = CC.y = A.y + 0.5 * (B.y - A.y);
178
179 if( m_Projector.Get_Projection(CC) )
180 {
181 Add_Segment(A, C, pLine);
182
183 pLine->Add_Point(CC);
184
185 Add_Segment(C, B, pLine);
186 }
187 }
188 else if( Length != NULL )
189 {
190 TSG_Point AA, BB;
191
192 if( m_Projector.Get_Projection(AA = A) && m_Projector.Get_Projection(BB = B) )
193 {
194 *Length += Get_Distance(AA, BB);
195 }
196 }
197 }
198
199
200 ///////////////////////////////////////////////////////////
201 // //
202 // //
203 // //
204 ///////////////////////////////////////////////////////////
205
206 //---------------------------------------------------------
CCRS_Distance_Lines(void)207 CCRS_Distance_Lines::CCRS_Distance_Lines(void)
208 {
209 //-----------------------------------------------------
210 Set_Name (_TL("Geographic Distances"));
211
212 Set_Author ("O. Conrad (c) 2015");
213
214 Set_Description (_TW(
215 "Calculates for all segments of the input lines the planar, great elliptic, "
216 "and loxodrome distance and re-projects the latter two to the projection "
217 "of the input lines. "
218 ));
219
220 Set_Description (Get_Description() + "\n" + CSG_CRSProjector::Get_Description());
221
222 //-----------------------------------------------------
223 Parameters.Add_Shapes(
224 NULL , "PLANAR" , _TL("Segments"),
225 _TL(""),
226 PARAMETER_INPUT, SHAPE_TYPE_Line
227 );
228
229 Parameters.Add_Shapes(
230 NULL , "ORTHODROME" , _TL("Great Elliptic"),
231 _TL(""),
232 PARAMETER_OUTPUT, SHAPE_TYPE_Line
233 );
234
235 Parameters.Add_Shapes(
236 NULL , "LOXODROME" , _TL("Loxodrome"),
237 _TL(""),
238 PARAMETER_OUTPUT, SHAPE_TYPE_Line
239 );
240
241 Parameters.Add_Value(
242 NULL , "EPSILON" , _TL("Epsilon"),
243 _TL("defines the maximum resolution [km] for the re-projected distance segments"),
244 PARAMETER_TYPE_Double, 100.0, 1.0, true
245 );
246 }
247
248
249 ///////////////////////////////////////////////////////////
250 // //
251 ///////////////////////////////////////////////////////////
252
253 //---------------------------------------------------------
On_Execute(void)254 bool CCRS_Distance_Lines::On_Execute(void)
255 {
256 //-----------------------------------------------------
257 CSG_Shapes *pPlanars = Parameters("PLANAR" )->asShapes();
258 CSG_Shapes *pOrthodromes = Parameters("ORTHODROME")->asShapes();
259 CSG_Shapes *pLoxodromes = Parameters("LOXODROME" )->asShapes();
260
261 //-----------------------------------------------------
262 CCRS_Distance_Calculator Distance;
263
264 if( !Distance.Create(pPlanars->Get_Projection(), Parameters("EPSILON")->asDouble() * 1000.0) )
265 {
266 Error_Set(_TL("projection initialization failed"));
267
268 return( false );
269 }
270
271 //-----------------------------------------------------
272 pOrthodromes->Create(pPlanars->Get_Type(), CSG_String::Format("%s [%s]", pPlanars->Get_Name(), _TL("Orthodromes")), pPlanars);
273 pOrthodromes->Add_Field("LENGTH_PLAN", SG_DATATYPE_Double);
274 pOrthodromes->Add_Field("LENGTH" , SG_DATATYPE_Double);
275
276 pLoxodromes ->Create(pPlanars->Get_Type(), CSG_String::Format("%s [%s]", pPlanars->Get_Name(), _TL("Loxodromes" )), pPlanars);
277 pLoxodromes ->Add_Field("LENGTH_PLAN", SG_DATATYPE_Double);
278 pLoxodromes ->Add_Field("LENGTH" , SG_DATATYPE_Double);
279
280 //-----------------------------------------------------
281 for(int iShape=0; iShape<pPlanars->Get_Count() && Set_Progress(iShape, pPlanars->Get_Count()); iShape++)
282 {
283 CSG_Shape_Line *pProj = (CSG_Shape_Line *)pPlanars->Get_Shape(iShape);
284
285 for(int iPart=0; iPart<pProj->Get_Part_Count(); iPart++)
286 {
287 if( pProj->Get_Point_Count(iPart) > 1 )
288 {
289 TSG_Point A, B = pProj->Get_Point(0, iPart);
290
291 CSG_Shape_Line *pOrthodrome = (CSG_Shape_Line *)pOrthodromes->Add_Shape(pProj, SHAPE_COPY_ATTR);
292 CSG_Shape_Line *pLoxodrome = (CSG_Shape_Line *)pLoxodromes->Add_Shape(pProj, SHAPE_COPY_ATTR);
293
294 pOrthodrome->Set_Value(pPlanars->Get_Field_Count() + 0, pProj->Get_Length(iPart));
295 pLoxodrome ->Set_Value(pPlanars->Get_Field_Count() + 0, pProj->Get_Length(iPart));
296
297 pOrthodrome->Add_Point(B);
298 pLoxodrome ->Add_Point(B);
299
300 double dOrthodrome = 0.0;
301 double dLoxodrome = 0.0;
302
303 for(int iPoint=1; iPoint<pProj->Get_Point_Count(iPart); iPoint++)
304 {
305 A = B; B = pProj->Get_Point(iPoint, iPart);
306
307 dOrthodrome += Distance.Get_Orthodrome(A, B, pOrthodrome);
308 dLoxodrome += Distance.Get_Loxodrome (A, B, pLoxodrome );
309
310 pOrthodrome->Add_Point(B);
311 pLoxodrome ->Add_Point(B);
312 }
313
314 pOrthodrome->Set_Value(pPlanars->Get_Field_Count() + 1, dOrthodrome);
315 pLoxodrome ->Set_Value(pPlanars->Get_Field_Count() + 1, dLoxodrome );
316 }
317 }
318 }
319
320 //-----------------------------------------------------
321 return( pOrthodromes->Get_Count() > 0 );
322 }
323
324
325 ///////////////////////////////////////////////////////////
326 // //
327 // //
328 // //
329 ///////////////////////////////////////////////////////////
330
331 //---------------------------------------------------------
CCRS_Distance_Points(void)332 CCRS_Distance_Points::CCRS_Distance_Points(void)
333 {
334 Set_Name (_TL("Geographic Distances (Pair of Coordinates)"));
335
336 Set_Author ("O. Conrad (c) 2015");
337
338 Set_Description (_TW(
339 "Calculates for all segments of the input lines the planar, great elliptic, "
340 "and loxodrome distance and re-projects the latter two to the projection "
341 "of the input lines. "
342 ));
343
344 Parameters.Add_Shapes(
345 NULL , "DISTANCES" , _TL("Geographic Distances"),
346 _TL(""),
347 PARAMETER_OUTPUT, SHAPE_TYPE_Line
348 );
349
350 CSG_Parameter *pNode;
351
352 pNode = Parameters.Add_Node(NULL, "NODE_A", _TL("From"), _TL(""));
353 Parameters.Add_Value(pNode, "COORD_X1", _TL("X"), _TL(""), PARAMETER_TYPE_Double, 10.0);
354 Parameters.Add_Value(pNode, "COORD_Y1", _TL("Y"), _TL(""), PARAMETER_TYPE_Double, 53.5);
355
356 pNode = Parameters.Add_Node(NULL, "NODE_B", _TL("To" ), _TL(""));
357 Parameters.Add_Value(pNode, "COORD_X2", _TL("X"), _TL(""), PARAMETER_TYPE_Double, 116.5);
358 Parameters.Add_Value(pNode, "COORD_Y2", _TL("Y"), _TL(""), PARAMETER_TYPE_Double, 6.4);
359
360 Parameters.Add_Value(
361 NULL , "EPSILON" , _TL("Epsilon"),
362 _TL("defines the maximum resolution [km] for the re-projected distance segments"),
363 PARAMETER_TYPE_Double, 100.0, 1.0, true
364 );
365 }
366
367
368 ///////////////////////////////////////////////////////////
369 // //
370 ///////////////////////////////////////////////////////////
371
372 //---------------------------------------------------------
On_Execute(void)373 bool CCRS_Distance_Points::On_Execute(void)
374 {
375 //-----------------------------------------------------
376 CSG_Projection Projection;
377
378 if( !Get_Projection(Projection) )
379 {
380 return( false );
381 }
382
383 //-----------------------------------------------------
384 CCRS_Distance_Calculator Distance;
385
386 if( !Distance.Create(Projection, Parameters("EPSILON")->asDouble() * 1000.0) )
387 {
388 Error_Set(_TL("projection initialization failed"));
389
390 return( false );
391 }
392
393 //-----------------------------------------------------
394 TSG_Point A, B;
395
396 A.x = Parameters("COORD_X1")->asDouble();
397 A.y = Parameters("COORD_Y1")->asDouble();
398
399 B.x = Parameters("COORD_X2")->asDouble();
400 B.y = Parameters("COORD_Y2")->asDouble();
401
402 //-----------------------------------------------------
403 CSG_Shapes *pDistances = Parameters("DISTANCES")->asShapes();
404
405 pDistances ->Create(SHAPE_TYPE_Line, CSG_String::Format("%s", _TL("Geographic Distances")));
406
407 pDistances ->Add_Field("TYPE" , SG_DATATYPE_String);
408 pDistances ->Add_Field("LENGTH", SG_DATATYPE_Double);
409
410 pDistances ->Get_Projection().Create(Projection);
411
412 //-----------------------------------------------------
413 CSG_Shape *pPlanar = pDistances->Add_Shape(); pPlanar ->Set_Value(0, "Planar" );
414 CSG_Shape *pOrthodrome = pDistances->Add_Shape(); pOrthodrome->Set_Value(0, "Orthodrome");
415 CSG_Shape *pLoxodrome = pDistances->Add_Shape(); pLoxodrome ->Set_Value(0, "Loxodrome" );
416
417 pPlanar ->Add_Point(A);
418 pOrthodrome->Add_Point(A);
419 pLoxodrome ->Add_Point(A);
420
421 pPlanar ->Set_Value(1, SG_Get_Distance (A, B ));
422 pOrthodrome->Set_Value(1, Distance.Get_Orthodrome(A, B, pOrthodrome));
423 pLoxodrome ->Set_Value(1, Distance.Get_Loxodrome (A, B, pLoxodrome ));
424
425 pPlanar ->Add_Point(B);
426 pOrthodrome->Add_Point(B);
427 pLoxodrome ->Add_Point(B);
428
429 //-----------------------------------------------------
430 return( true );
431 }
432
433
434 ///////////////////////////////////////////////////////////
435 // //
436 // //
437 // //
438 ///////////////////////////////////////////////////////////
439
440 //---------------------------------------------------------
CCRS_Distance_Interactive(void)441 CCRS_Distance_Interactive::CCRS_Distance_Interactive(void)
442 {
443 Set_Name (_TL("Geographic Distances"));
444
445 Set_Author ("O. Conrad (c) 2015");
446
447 Set_Description (_TW(
448 "Calculates for all segments of the input lines the planar, great elliptic, "
449 "and loxodrome distance and re-projects the latter two to the projection "
450 "of the input lines. "
451 ));
452
453 Parameters.Add_Shapes(
454 NULL , "DISTANCES" , _TL("Geographic Distances"),
455 _TL(""),
456 PARAMETER_OUTPUT, SHAPE_TYPE_Line
457 );
458
459 Parameters.Add_Value(
460 NULL , "EPSILON" , _TL("Epsilon"),
461 _TL("defines the maximum resolution [km] for the re-projected distance segments"),
462 PARAMETER_TYPE_Double, 100.0, 1.0, true
463 );
464
465 Set_Drag_Mode(TOOL_INTERACTIVE_DRAG_LINE);
466 }
467
468 //---------------------------------------------------------
On_Execute(void)469 bool CCRS_Distance_Interactive::On_Execute(void)
470 {
471 CCRS_Picker CRS;
472
473 if( !Dlg_Parameters(CRS.Get_Parameters(), CRS.Get_Name()) )
474 {
475 m_Projection.Destroy();
476
477 return( false );
478 }
479
480 return( m_Projection.Create(CRS.Get_Parameters()->Get_Parameter("CRS_PROJ4")->asString(), SG_PROJ_FMT_Proj4) );
481 }
482
483 //---------------------------------------------------------
On_Execute_Position(CSG_Point ptWorld,TSG_Tool_Interactive_Mode Mode)484 bool CCRS_Distance_Interactive::On_Execute_Position(CSG_Point ptWorld, TSG_Tool_Interactive_Mode Mode)
485 {
486 if( Mode == TOOL_INTERACTIVE_LDOWN )
487 {
488 m_Down = ptWorld;
489 }
490 else if( Mode == TOOL_INTERACTIVE_LUP )
491 {
492 if( m_Down != ptWorld )
493 {
494 CCRS_Distance_Points Distance;
495
496 Distance.Set_Parameter("DISTANCES", Parameters("DISTANCES")->asShapes());
497 Distance.Set_Parameter("EPSILON" , Parameters("EPSILON" )->asDouble());
498 Distance.Set_Parameter("CRS_PROJ4", m_Projection.Get_Proj4());
499 Distance.Set_Parameter("COORD_X1" , m_Down .Get_X());
500 Distance.Set_Parameter("COORD_Y1" , m_Down .Get_Y());
501 Distance.Set_Parameter("COORD_X2" , ptWorld.Get_X());
502 Distance.Set_Parameter("COORD_Y2" , ptWorld.Get_Y());
503
504 Distance.Execute();
505
506 DataObject_Update(Parameters("DISTANCES")->asShapes());
507 }
508 }
509
510 return( true );
511 }
512
513
514 ///////////////////////////////////////////////////////////
515 // //
516 // //
517 // //
518 ///////////////////////////////////////////////////////////
519
520 //---------------------------------------------------------
521