1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // Projection_Proj4 //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // gcs_graticule.cpp //
14 // //
15 // Copyright (C) 2014 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 #include "gcs_graticule.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
59 #define AXIS_LEFT 1
60 #define AXIS_RIGHT 2
61 #define AXIS_BOTTOM 3
62 #define AXIS_TOP 4
63
64 //---------------------------------------------------------
65 enum
66 {
67 DEG_PREC_AUTO,
68 DEG_PREC_FULL,
69 DEG_PREC_SEC,
70 DEG_PREC_MIN,
71 DEG_PREC_DEG
72 };
73
74
75 ///////////////////////////////////////////////////////////
76 // //
77 // //
78 // //
79 ///////////////////////////////////////////////////////////
80
81 //---------------------------------------------------------
CGCS_Graticule(void)82 CGCS_Graticule::CGCS_Graticule(void)
83 {
84 Set_Name (_TL("Latitude/Longitude Graticule"));
85
86 Set_Author ("O.Conrad (c) 2014");
87
88 Set_Description (_TW(
89 "Creates a longitude/latitude graticule for the extent and projection of the input shapes layer. "
90 ));
91
92 Set_Description (Get_Description() + "\n" + CSG_CRSProjector::Get_Description());
93
94 //-----------------------------------------------------
95 Parameters.Add_Shapes("",
96 "GRATICULE" , _TL("Graticule"),
97 _TL(""),
98 PARAMETER_OUTPUT, SHAPE_TYPE_Line
99 );
100
101 Parameters.Add_Shapes("",
102 "COORDS" , _TL("Frame Coordinates"),
103 _TL(""),
104 PARAMETER_OUTPUT_OPTIONAL, SHAPE_TYPE_Point
105 );
106
107 Parameters.Add_Node("",
108 "NODE_GRID" , _TL("Graticule"),
109 _TL("")
110 );
111
112 Parameters.Add_Node("NODE_GRID", "NODE_X", _TL("X Range"), _TL(""));
113 Parameters.Add_Double("NODE_X" , "XMIN" , _TL("Minimum"), _TL(""));
114 Parameters.Add_Double("NODE_X" , "XMAX" , _TL("Maximum"), _TL(""));
115
116 Parameters.Add_Node("NODE_GRID", "NODE_Y", _TL("Y Range"), _TL(""));
117 Parameters.Add_Double("NODE_Y" , "YMIN" , _TL("Minimum"), _TL(""));
118 Parameters.Add_Double("NODE_Y" , "YMAX" , _TL("Maximum"), _TL(""));
119
120 Parameters.Add_Choice("NODE_GRID",
121 "INTERVAL" , _TL("Interval"),
122 _TL(""),
123 CSG_String::Format("%s|%s",
124 _TL("fixed interval"),
125 _TL("fitted interval")
126 ), 0
127 );
128
129 Parameters.Add_Double("NODE_GRID",
130 "FIXED" , _TL("Fixed Interval (Degree)"),
131 _TL(""),
132 1., 0., true, 20.
133 );
134
135 Parameters.Add_Int("NODE_GRID",
136 "FITTED" , _TL("Number of Intervals"),
137 _TL(""),
138 10, 1, true
139 );
140
141 Parameters.Add_Double("NODE_GRID",
142 "RESOLUTION", _TL("Minimum Resolution (Degree)"),
143 _TL(""),
144 0.5, 0., true
145 );
146 }
147
148
149 ///////////////////////////////////////////////////////////
150 // //
151 ///////////////////////////////////////////////////////////
152
153 //---------------------------------------------------------
On_Parameter_Changed(CSG_Parameters * pParameters,CSG_Parameter * pParameter)154 int CGCS_Graticule::On_Parameter_Changed(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
155 {
156 if( pParameter->Cmp_Identifier("CRS_GRID" )
157 || pParameter->Cmp_Identifier("CRS_SHAPES") )
158 {
159 CSG_Rect r(pParameter->Cmp_Identifier("CRS_GRID")
160 ? pParameter->asParameters()->Get_Parameter("PICK")->asGrid ()->Get_Extent()
161 : pParameter->asParameters()->Get_Parameter("PICK")->asShapes()->Get_Extent()
162 );
163
164 if( r.Get_XRange() > 0. && r.Get_YRange() > 0. )
165 {
166 pParameters->Set_Parameter("XMIN", r.Get_XMin());
167 pParameters->Set_Parameter("XMAX", r.Get_XMax());
168 pParameters->Set_Parameter("YMIN", r.Get_YMin());
169 pParameters->Set_Parameter("YMAX", r.Get_YMax());
170 }
171 }
172
173 return( CCRS_Base::On_Parameter_Changed(pParameters, pParameter) );
174 }
175
176 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)177 int CGCS_Graticule::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
178 {
179 if( pParameter->Cmp_Identifier("INTERVAL") )
180 {
181 pParameters->Set_Enabled("FIXED" , pParameter->asInt() == 0);
182 pParameters->Set_Enabled("FITTED", pParameter->asInt() == 1);
183 }
184
185 return( CCRS_Base::On_Parameters_Enable(pParameters, pParameter) );
186 }
187
188
189 ///////////////////////////////////////////////////////////
190 // //
191 ///////////////////////////////////////////////////////////
192
193 //---------------------------------------------------------
On_Execute(void)194 bool CGCS_Graticule::On_Execute(void)
195 {
196 //-----------------------------------------------------
197 CSG_Projection Projection;
198
199 if( !Get_Projection(Projection) )
200 {
201 return( false );
202 }
203
204 //-----------------------------------------------------
205 m_Projector.Set_Source(CSG_Projection("+proj=longlat +ellps=WGS84 +datum=WGS84", SG_PROJ_FMT_Proj4));
206
207 if( !m_Projector.Set_Target(Projection) )
208 {
209 m_Projector.Destroy();
210
211 return( false );
212 }
213
214 //-----------------------------------------------------
215 CSG_Rect Extent(
216 Parameters("XMIN")->asDouble(),
217 Parameters("YMIN")->asDouble(),
218 Parameters("XMAX")->asDouble(),
219 Parameters("YMAX")->asDouble()
220 );
221
222 if( !Get_Graticule(Extent) )
223 {
224 m_Projector.Destroy();
225
226 return( false );
227 }
228
229 //-----------------------------------------------------
230 m_Projector.Destroy();
231
232 return( true );
233 }
234
235
236 ///////////////////////////////////////////////////////////
237 // //
238 ///////////////////////////////////////////////////////////
239
240 //---------------------------------------------------------
Get_Graticule(const CSG_Rect & Extent)241 bool CGCS_Graticule::Get_Graticule(const CSG_Rect &Extent)
242 {
243 double x, y, Interval;
244 CSG_Rect r;
245
246 if( !Get_Extent(Extent, r) || (Interval = Get_Interval(r)) <= 0. )
247 {
248 return( false );
249 }
250
251 //-----------------------------------------------------
252 r.m_rect.xMin = Interval * floor(r.Get_XMin() / Interval);
253 r.m_rect.xMax = Interval * ceil (r.Get_XMax() / Interval);
254 r.m_rect.yMin = Interval * floor(r.Get_YMin() / Interval);
255 r.m_rect.yMax = Interval * ceil (r.Get_YMax() / Interval);
256
257 r.Inflate(Interval, false);
258
259 if( r.Get_XMin() < -180. ) r.m_rect.xMin = -180.;
260 if( r.Get_XMax() > 180. ) r.m_rect.xMax = 180.;
261 if( r.Get_YMin() < -90. ) r.m_rect.yMin = -90.;
262 if( r.Get_YMax() > 90. ) r.m_rect.yMax = 90.;
263
264 //-----------------------------------------------------
265 double Resolution = Parameters("RESOLUTION")->asDouble(); if( Resolution <= 0. ) Resolution = Interval;
266
267 if( Interval > Resolution )
268 {
269 Resolution = Interval / ceil(Interval / Resolution);
270 }
271
272 //-----------------------------------------------------
273 CSG_Shapes *pGraticule = Parameters("GRATICULE")->asShapes();
274
275 pGraticule->Create(SHAPE_TYPE_Line);
276 pGraticule->Set_Name(_TL("Graticule"));
277
278 pGraticule->Add_Field("TYPE" , SG_DATATYPE_String);
279 pGraticule->Add_Field("LABEL" , SG_DATATYPE_String);
280 pGraticule->Add_Field("DEGREE", SG_DATATYPE_Double);
281
282 //-----------------------------------------------------
283 CSG_Shapes *pCoordinates = Parameters("COORDS")->asShapes();
284
285 if( pCoordinates )
286 {
287 pCoordinates->Create(SHAPE_TYPE_Point);
288 pCoordinates->Set_Name(_TL("Coordinates"));
289
290 pCoordinates->Add_Field("TYPE" , SG_DATATYPE_String);
291 pCoordinates->Add_Field("LABEL", SG_DATATYPE_String);
292 }
293
294 //-----------------------------------------------------
295 CSG_Shapes Clip(SHAPE_TYPE_Polygon);
296 CSG_Shape *pClip = Clip.Add_Shape();
297
298 pClip->Add_Point(Extent.Get_XMin(), Extent.Get_YMin());
299 pClip->Add_Point(Extent.Get_XMin(), Extent.Get_YMax());
300 pClip->Add_Point(Extent.Get_XMax(), Extent.Get_YMax());
301 pClip->Add_Point(Extent.Get_XMax(), Extent.Get_YMin());
302 pClip->Add_Point(Extent.Get_XMin(), Extent.Get_YMin());
303
304 //-----------------------------------------------------
305 for(y=r.Get_YMin(); y<=r.Get_YMax(); y+=Interval)
306 {
307 CSG_Shape *pLine = pGraticule->Add_Shape();
308
309 pLine->Set_Value(0, "LAT");
310 pLine->Set_Value(1, Get_Degree(y, DEG_PREC_DEG));
311 pLine->Set_Value(2, y);
312
313 for(x=r.Get_XMin(); x<=r.Get_XMax(); x+=Interval)
314 {
315 CSG_Point p(x, y); m_Projector.Get_Projection(p); pLine->Add_Point(p);
316
317 if( Resolution < Interval && x < r.Get_XMax() )
318 {
319 for(double i=x+Resolution; i<x+Interval; i+=Resolution)
320 {
321 CSG_Point p(i, y); m_Projector.Get_Projection(p); pLine->Add_Point(p);
322 }
323 }
324 }
325
326 Get_Coordinate(Extent, pCoordinates, pLine, AXIS_LEFT);
327 Get_Coordinate(Extent, pCoordinates, pLine, AXIS_RIGHT);
328
329 if( !SG_Polygon_Intersection(pLine, pClip) )
330 {
331 pGraticule->Del_Shape(pLine);
332 }
333 }
334
335 //-----------------------------------------------------
336 for(x=r.Get_XMin(); x<=r.Get_XMax(); x+=Interval)
337 {
338 CSG_Shape *pLine = pGraticule->Add_Shape();
339
340 pLine->Set_Value(0, "LON");
341 pLine->Set_Value(1, Get_Degree(x, DEG_PREC_DEG));
342 pLine->Set_Value(2, x);
343
344 for(y=r.Get_YMin(); y<=r.Get_YMax(); y+=Interval)
345 {
346 CSG_Point p(x, y); m_Projector.Get_Projection(p); pLine->Add_Point(p);
347
348 if( Resolution < Interval && y < r.Get_YMax() )
349 {
350 for(double i=y+Resolution; i<y+Interval; i+=Resolution)
351 {
352 CSG_Point p(x, i); m_Projector.Get_Projection(p); pLine->Add_Point(p);
353 }
354 }
355 }
356
357 Get_Coordinate(Extent, pCoordinates, pLine, AXIS_BOTTOM);
358 Get_Coordinate(Extent, pCoordinates, pLine, AXIS_TOP);
359
360 if( !SG_Polygon_Intersection(pLine, pClip) )
361 {
362 pGraticule->Del_Shape(pLine);
363 }
364 }
365
366 //-----------------------------------------------------
367 return( true );
368 }
369
370
371 ///////////////////////////////////////////////////////////
372 // //
373 ///////////////////////////////////////////////////////////
374
375 //---------------------------------------------------------
Get_Coordinate(const CSG_Rect & Extent,CSG_Shapes * pCoordinates,CSG_Shape * pLine,int Axis)376 bool CGCS_Graticule::Get_Coordinate(const CSG_Rect &Extent, CSG_Shapes *pCoordinates, CSG_Shape *pLine, int Axis)
377 {
378 if( !pCoordinates || !Extent.Intersects(pLine->Get_Extent()) || pLine->Get_Point_Count(0) < 2 )
379 {
380 return( false );
381 }
382
383 TSG_Point A[2], B[2], C;
384
385 switch( Axis )
386 {
387 case AXIS_LEFT : A[0].x = A[1].x = Extent.Get_XMin(); A[0].y = Extent.Get_YMin(); A[1].y = Extent.Get_YMax(); break;
388 case AXIS_RIGHT : A[0].x = A[1].x = Extent.Get_XMax(); A[0].y = Extent.Get_YMin(); A[1].y = Extent.Get_YMax(); break;
389 case AXIS_BOTTOM: A[0].y = A[1].y = Extent.Get_YMin(); A[0].x = Extent.Get_XMin(); A[1].x = Extent.Get_XMax(); break;
390 case AXIS_TOP : A[0].y = A[1].y = Extent.Get_YMax(); A[0].x = Extent.Get_XMin(); A[1].x = Extent.Get_XMax(); break;
391
392 default:
393 return( false );
394 }
395
396 //-----------------------------------------------------
397 B[1] = pLine->Get_Point(0);
398
399 for(int i=1; i<pLine->Get_Point_Count(); i++)
400 {
401 B[0] = B[1];
402 B[1] = pLine->Get_Point(i);
403
404 if( SG_Get_Crossing(C, A[0], A[1], B[0], B[1], true) )
405 {
406 CSG_Shape *pPoint = pCoordinates->Add_Shape();
407 pPoint->Add_Point(C);
408 pPoint->Set_Value(0, CSG_String(pLine->asString(0)) + (Axis == AXIS_LEFT || Axis == AXIS_BOTTOM ? "_MIN" : "_MAX"));
409 pPoint->Set_Value(1, pLine->asString(1));
410
411 return( true );
412 }
413 }
414
415 //-----------------------------------------------------
416 switch( Axis )
417 {
418 case AXIS_LEFT : C = pLine->Get_Point(0, 0, true ); break;
419 case AXIS_RIGHT : C = pLine->Get_Point(0, 0, false); break;
420 case AXIS_BOTTOM: C = pLine->Get_Point(0, 0, true ); break;
421 case AXIS_TOP : C = pLine->Get_Point(0, 0, false); break;
422 }
423
424 if( Extent.Contains(C) )
425 {
426 CSG_Shape *pPoint = pCoordinates->Add_Shape();
427 pPoint->Add_Point(C);
428 pPoint->Set_Value(0, CSG_String(pLine->asString(0)) + (Axis == AXIS_LEFT || Axis == AXIS_BOTTOM ? "_MIN" : "_MAX"));
429 pPoint->Set_Value(1, pLine->asString(1));
430
431 return( true );
432 }
433
434 //-----------------------------------------------------
435 return( false );
436 }
437
438
439 ///////////////////////////////////////////////////////////
440 // //
441 ///////////////////////////////////////////////////////////
442
443 //---------------------------------------------------------
Get_Interval(const CSG_Rect & Extent)444 double CGCS_Graticule::Get_Interval(const CSG_Rect &Extent)
445 {
446 if( Parameters("INTERVAL")->asInt() == 0 )
447 {
448 return( Parameters("FIXED")->asDouble() );
449 }
450
451 double Interval = Extent.Get_XRange() > Extent.Get_YRange() ? Extent.Get_XRange() : Extent.Get_YRange();
452
453 if( Interval > 360. )
454 {
455 Interval = 360.;
456 }
457
458 Interval = Interval / Parameters("FITTED")->asInt();
459
460 double d = pow(10., (int)(log10(Interval)) - (Interval < 1. ? 1. : 0.));
461
462 Interval = (int)(Interval / d) * d;
463
464 return( Interval );
465 }
466
467
468 ///////////////////////////////////////////////////////////
469 // //
470 ///////////////////////////////////////////////////////////
471
472 //---------------------------------------------------------
Get_Extent(const CSG_Rect & Extent,CSG_Rect & r)473 bool CGCS_Graticule::Get_Extent(const CSG_Rect &Extent, CSG_Rect &r)
474 {
475 if( m_Projector.Set_Inverse() )
476 {
477 double x, y, d;
478
479 CSG_Point p(Extent.Get_XMin(), Extent.Get_YMin());
480
481 m_Projector.Get_Projection(p); r.Assign(p, p);
482
483 d = Extent.Get_XRange() / 10.;
484
485 for(y=Extent.Get_YMin(), x=Extent.Get_XMin(); x<=Extent.Get_XMax(); x+=d)
486 {
487 p.Assign(x, y); m_Projector.Get_Projection(p); r.Union(p);
488 }
489
490 for(y=Extent.Get_YMax(), x=Extent.Get_XMin(); x<=Extent.Get_XMax(); x+=d)
491 {
492 p.Assign(x, y); m_Projector.Get_Projection(p); r.Union(p);
493 }
494
495 d = Extent.Get_YRange() / 10.;
496
497 for(x=Extent.Get_XMin(), y=Extent.Get_YMin(); y<=Extent.Get_YMax(); y+=d)
498 {
499 p.Assign(x, y); m_Projector.Get_Projection(p); r.Union(p);
500 }
501
502 for(x=Extent.Get_XMax(), y=Extent.Get_YMin(); y<=Extent.Get_YMax(); y+=d)
503 {
504 p.Assign(x, y); m_Projector.Get_Projection(p); r.Union(p);
505 }
506
507 m_Projector.Set_Inverse(false);
508
509 if( r.Get_XMin() < -180. ) r.m_rect.xMin = -180.; else if( r.Get_XMax() > 180. ) r.m_rect.xMax = 180.;
510 if( r.Get_YMin() < -90. ) r.m_rect.yMin = -90.; else if( r.Get_YMax() > 90. ) r.m_rect.yMax = 90.;
511
512 return( r.Get_XRange() > 0. && r.Get_YRange() > 0. );
513 }
514
515 return( false );
516 }
517
518
519 ///////////////////////////////////////////////////////////
520 // //
521 ///////////////////////////////////////////////////////////
522
523 //---------------------------------------------------------
Get_Degree(double Value,int Precision)524 CSG_String CGCS_Graticule::Get_Degree(double Value, int Precision)
525 {
526 if( Precision == DEG_PREC_DEG )
527 {
528 return( SG_Get_String(Value, -12) + "\xb0" );
529 }
530
531 SG_Char c;
532 int d, h;
533 double s;
534 CSG_String String;
535
536 if( Value < 0. )
537 {
538 Value = -Value;
539 c = '-';
540 }
541 else
542 {
543 c = '+';
544 }
545
546 Value = fmod(Value, 360.);
547 d = (int)Value;
548 Value = 60. * (Value - d);
549 h = (int)Value;
550 Value = 60. * (Value - h);
551 s = Value;
552
553 if( s > 0. || Precision == DEG_PREC_FULL )
554 {
555 String.Printf("%c%d\xb0%02d'%02.*f''", c, d, h, SG_Get_Significant_Decimals(s), s);
556 }
557 else if( h > 0 || Precision == DEG_PREC_MIN )
558 {
559 String.Printf("%c%d\xb0%02d'" , c, d, h);
560 }
561 else
562 {
563 String.Printf("%c%d\xb0" , c, d);
564 }
565
566 return( String );
567 }
568
569
570 ///////////////////////////////////////////////////////////
571 // //
572 // //
573 // //
574 ///////////////////////////////////////////////////////////
575
576 //---------------------------------------------------------
577