1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // Grid_Gridding //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // Shapes2Grid.cpp //
14 // //
15 // Copyright (C) 2003 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 Goettingen //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49 #include "Shapes2Grid.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
59 #define X_WORLD_TO_GRID(X) (((X) - m_pGrid->Get_XMin()) / m_pGrid->Get_Cellsize())
60 #define Y_WORLD_TO_GRID(Y) (((Y) - m_pGrid->Get_YMin()) / m_pGrid->Get_Cellsize())
61
62 //---------------------------------------------------------
63 #define OUTPUT_NODATA -2
64 #define OUTPUT_INDEX -1
65
66
67 ///////////////////////////////////////////////////////////
68 // //
69 // //
70 // //
71 ///////////////////////////////////////////////////////////
72
73 //---------------------------------------------------------
CShapes2Grid(void)74 CShapes2Grid::CShapes2Grid(void)
75 {
76 Set_Name (_TL("Shapes to Grid"));
77
78 Set_Author ("O.Conrad (c) 2003");
79
80 Set_Description (_TW(
81 "Gridding of a shapes layer. If some shapes are selected, only these will be gridded."
82 ));
83
84 //-----------------------------------------------------
85 Parameters.Add_Shapes("",
86 "INPUT" , _TL("Shapes"),
87 _TL(""),
88 PARAMETER_INPUT
89 );
90
91 Parameters.Add_Table_Field("INPUT",
92 "FIELD" , _TL("Attribute"),
93 _TL("")
94 );
95
96 Parameters.Add_Choice("",
97 "OUTPUT" , _TL("Output Values"),
98 _TL(""),
99 CSG_String::Format("%s|%s|%s",
100 _TL("data / no-data"),
101 _TL("index number"),
102 _TL("attribute")
103 ), 2
104 );
105
106 Parameters.Add_Choice("",
107 "MULTIPLE" , _TL("Method for Multiple Values"),
108 _TL(""),
109 CSG_String::Format("%s|%s|%s|%s|%s",
110 _TL("first"),
111 _TL("last"),
112 _TL("minimum"),
113 _TL("maximum"),
114 _TL("mean")
115 ), 1
116 );
117
118 Parameters.Add_Choice("",
119 "LINE_TYPE" , _TL("Lines"),
120 _TL(""),
121 CSG_String::Format("%s|%s",
122 _TL("thin"),
123 _TL("thick")
124 ), 1
125 );
126
127 Parameters.Add_Choice("",
128 "POLY_TYPE" , _TL("Polygon"),
129 _TL(""),
130 CSG_String::Format("%s|%s",
131 _TL("node"),
132 _TL("cell")
133 ), 1
134 );
135
136 Parameters.Add_Choice("",
137 "GRID_TYPE" , _TL("Data Type"),
138 _TL(""),
139 CSG_String::Format("%s|%s|%s|%s|%s|%s|%s|%s|%s|%s",
140 _TL("1 bit"),
141 _TL("1 byte unsigned integer"),
142 _TL("1 byte signed integer"),
143 _TL("2 byte unsigned integer"),
144 _TL("2 byte signed integer"),
145 _TL("4 byte unsigned integer"),
146 _TL("4 byte signed integer"),
147 _TL("4 byte floating point"),
148 _TL("8 byte floating point"),
149 _TL("same as attribute")
150 ), 9
151 );
152
153 //-----------------------------------------------------
154 m_Grid_Target.Create(&Parameters, false, "", "TARGET_");
155
156 m_Grid_Target.Add_Grid("GRID" , _TL("Grid" ), false);
157 m_Grid_Target.Add_Grid("COUNT", _TL("Number of Values"), true);
158 }
159
160
161 ///////////////////////////////////////////////////////////
162 // //
163 ///////////////////////////////////////////////////////////
164
165 //---------------------------------------------------------
On_Parameter_Changed(CSG_Parameters * pParameters,CSG_Parameter * pParameter)166 int CShapes2Grid::On_Parameter_Changed(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
167 {
168 if( pParameter->Cmp_Identifier("INPUT") )
169 {
170 m_Grid_Target.Set_User_Defined(pParameters, pParameter->asShapes());
171 }
172
173 m_Grid_Target.On_Parameter_Changed(pParameters, pParameter);
174
175 return( CSG_Tool::On_Parameter_Changed(pParameters, pParameter) );
176 }
177
178 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)179 int CShapes2Grid::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
180 {
181 if( pParameter->Cmp_Identifier("INPUT") )
182 {
183 pParameters->Set_Enabled("LINE_TYPE", pParameter->asShapes() && pParameter->asShapes()->Get_Type() == SHAPE_TYPE_Line);
184 pParameters->Set_Enabled("POLY_TYPE", pParameter->asShapes() && pParameter->asShapes()->Get_Type() == SHAPE_TYPE_Polygon);
185 }
186
187 if( pParameter->Cmp_Identifier("OUTPUT") )
188 {
189 pParameters->Set_Enabled("FIELD" , pParameter->asInt() == 2);
190 pParameters->Set_Enabled("MULTIPLE" , pParameter->asInt() != 0);
191 pParameters->Set_Enabled("GRID_TYPE", pParameter->asInt() == 2);
192 }
193
194 m_Grid_Target.On_Parameters_Enable(pParameters, pParameter);
195
196 return( CSG_Tool::On_Parameters_Enable(pParameters, pParameter) );
197 }
198
199
200 ///////////////////////////////////////////////////////////
201 // //
202 ///////////////////////////////////////////////////////////
203
204 //---------------------------------------------------------
Get_Data_Type(int Field)205 TSG_Data_Type CShapes2Grid::Get_Data_Type(int Field)
206 {
207 CSG_Shapes *pShapes = Parameters("INPUT")->asShapes();
208
209 if( Field >= 0 && (Field >= pShapes->Get_Field_Count() || !SG_Data_Type_is_Numeric(pShapes->Get_Field_Type(Field))) )
210 {
211 Field = OUTPUT_INDEX; // index number
212 }
213
214 if( Field == OUTPUT_NODATA ) // data / no-data
215 {
216 return( SG_DATATYPE_Byte );
217 }
218
219 if( Field <= OUTPUT_INDEX ) // index number
220 {
221 return( pShapes->Get_Count() < 65535 ? SG_DATATYPE_Word : SG_DATATYPE_DWord );
222 }
223
224 // if( Field >= 0 ) // attribute
225 {
226 switch( Parameters("GRID_TYPE")->asInt() )
227 {
228 case 0: return( SG_DATATYPE_Bit );
229 case 1: return( SG_DATATYPE_Byte );
230 case 2: return( SG_DATATYPE_Char );
231 case 3: return( SG_DATATYPE_Word );
232 case 4: return( SG_DATATYPE_Short );
233 case 5: return( SG_DATATYPE_DWord );
234 case 6: return( SG_DATATYPE_Int );
235 case 7: return( SG_DATATYPE_Float );
236 case 8: return( SG_DATATYPE_Double );
237 }
238
239 return( pShapes->Get_Field_Type(Field) );
240 }
241 }
242
243
244 ///////////////////////////////////////////////////////////
245 // //
246 ///////////////////////////////////////////////////////////
247
248 //---------------------------------------------------------
On_Execute(void)249 bool CShapes2Grid::On_Execute(void)
250 {
251 CSG_Shapes *pShapes = Parameters("INPUT")->asShapes();
252
253 m_Multiple = Parameters("MULTIPLE")->asInt();
254
255 //-----------------------------------------------------
256 bool bFat;
257
258 switch( pShapes->Get_Type() )
259 {
260 default : bFat = false; break;
261 case SHAPE_TYPE_Line : bFat = Parameters("LINE_TYPE")->asInt() == 1; break;
262 case SHAPE_TYPE_Polygon: bFat = Parameters("POLY_TYPE")->asInt() == 1; break;
263 }
264
265 //-----------------------------------------------------
266 int Field;
267
268 switch( Parameters("OUTPUT")->asInt() )
269 {
270 case 0: Field = OUTPUT_NODATA; break; // data / no-data
271 case 1: Field = OUTPUT_INDEX ; break; // index number
272 default: Field = Parameters("FIELD")->asInt(); // attribute
273 if( Field < 0 || !SG_Data_Type_is_Numeric(pShapes->Get_Field_Type(Field)) )
274 {
275 Message_Add(_TL("WARNING: selected attribute is not numeric."));
276 }
277 break;
278 }
279
280 //-----------------------------------------------------
281 if( (m_pGrid = m_Grid_Target.Get_Grid("GRID", Get_Data_Type(Field))) == NULL )
282 {
283 return( false );
284 }
285
286 if( !pShapes->Get_Extent().Intersects(m_pGrid->Get_Extent()) )
287 {
288 Error_Set(_TL("Polygons' and target grid's extent do not intersect."));
289
290 return( false );
291 }
292
293 if( Field < 0 )
294 {
295 m_pGrid->Set_NoData_Value(0.);
296 }
297
298 m_pGrid->Fmt_Name("%s [%s]", pShapes->Get_Name(), Field < 0 ? _TL("ID") : pShapes->Get_Field_Name(Field));
299 m_pGrid->Assign_NoData();
300
301 //-------------------------------------------------
302 CSG_Grid Count;
303
304 m_pCount = m_Grid_Target.Get_Grid("COUNT", pShapes->Get_Count() < 256 ? SG_DATATYPE_Byte : SG_DATATYPE_Word);
305
306 if( m_pCount == NULL )
307 {
308 Count.Create(m_pGrid->Get_System(), SG_DATATYPE_Word);
309
310 m_pCount = &Count;
311 }
312
313 m_pCount->Fmt_Name("%s [%s]", pShapes->Get_Name(), _TL("Count"));
314 m_pCount->Set_NoData_Value(0.);
315 m_pCount->Assign(0.);
316
317 //-----------------------------------------------------
318 for(int i=0; i<pShapes->Get_Count() && Set_Progress(i, pShapes->Get_Count()); i++)
319 {
320 CSG_Shape *pShape = pShapes->Get_Shape(i);
321
322 m_Cells_On_Shape.clear();
323
324 if( pShapes->Get_Selection_Count() <= 0 || pShape->is_Selected() )
325 {
326 if( Field < 0 || !pShape->is_NoData(Field) )
327 {
328 if( pShape->Intersects(m_pGrid->Get_Extent()) )
329 {
330 double Value = Field >= 0 ? pShape->asDouble(Field) : Field == OUTPUT_INDEX ? i + 1 : 1;
331
332 switch( pShapes->Get_Type() )
333 {
334 case SHAPE_TYPE_Point: case SHAPE_TYPE_Points:
335 Set_Points (pShape, Value);
336 break;
337
338 case SHAPE_TYPE_Line:
339 Set_Line (pShape, bFat, Value);
340 break;
341
342 case SHAPE_TYPE_Polygon:
343 Set_Polygon (pShape, bFat, Value);
344 break;
345 }
346 }
347 }
348 }
349 }
350
351 //-----------------------------------------------------
352 if( m_Multiple == 4 ) // mean
353 {
354 for(int y=0; y<m_pGrid->Get_NY() && Set_Progress(y, m_pGrid->Get_NY()); y++)
355 {
356 for(int x=0; x<m_pGrid->Get_NX(); x++)
357 {
358 if( m_pCount->asInt(x, y) > 1 )
359 {
360 m_pGrid->Mul_Value(x, y, 1.0 / m_pCount->asDouble(x, y));
361 }
362 }
363 }
364 }
365
366 //-----------------------------------------------------
367 return( true );
368 }
369
370
371 ///////////////////////////////////////////////////////////
372 // //
373 ///////////////////////////////////////////////////////////
374
375 //---------------------------------------------------------
Set_Value(int x,int y,double Value,bool bCheckDuplicates)376 inline void CShapes2Grid::Set_Value(int x, int y, double Value, bool bCheckDuplicates)
377 {
378 if( bCheckDuplicates )
379 {
380 sLong n = y * m_pGrid->Get_NX() + x;
381
382 if( !m_Cells_On_Shape.insert(n).second )
383 {
384 return; // this cell has already been rendered for this shape
385 }
386 }
387
388 if( m_pGrid->is_InGrid(x, y, false) )
389 {
390 if( m_pCount->asInt(x, y) == 0 )
391 {
392 m_pGrid->Set_Value(x, y, Value);
393 }
394 else switch( m_Multiple )
395 {
396 default: // first
397 break;
398
399 case 1: // last
400 m_pGrid->Set_Value(x, y, Value);
401 break;
402
403 case 2: // minimum
404 if( m_pGrid->asDouble(x, y) > Value )
405 {
406 m_pGrid->Set_Value(x, y, Value);
407 }
408 break;
409
410 case 3: // maximum
411 if( m_pGrid->asDouble(x, y) < Value )
412 {
413 m_pGrid->Set_Value(x, y, Value);
414 }
415 break;
416
417 case 4: // mean
418 m_pGrid->Add_Value(x, y, Value);
419 break;
420 }
421
422 m_pCount->Add_Value(x, y, 1);
423 }
424 }
425
426
427 ///////////////////////////////////////////////////////////
428 // //
429 ///////////////////////////////////////////////////////////
430
431 //---------------------------------------------------------
Set_Points(CSG_Shape * pShape,double Value)432 void CShapes2Grid::Set_Points(CSG_Shape *pShape, double Value)
433 {
434 for(int iPart=0; iPart<pShape->Get_Part_Count(); iPart++)
435 {
436 for(int iPoint=0; iPoint<pShape->Get_Point_Count(iPart); iPoint++)
437 {
438 TSG_Point p = pShape->Get_Point(iPoint, iPart);
439
440 Set_Value(
441 (int)(0.5 + X_WORLD_TO_GRID(p.x)),
442 (int)(0.5 + Y_WORLD_TO_GRID(p.y)), Value,
443 false
444 );
445 }
446 }
447 }
448
449
450 ///////////////////////////////////////////////////////////
451 // //
452 ///////////////////////////////////////////////////////////
453
454 //---------------------------------------------------------
Set_Line(CSG_Shape * pShape,bool bFat,double Value)455 void CShapes2Grid::Set_Line(CSG_Shape *pShape, bool bFat, double Value)
456 {
457 for(int iPart=0; iPart<pShape->Get_Part_Count(); iPart++)
458 {
459 if( !((CSG_Shape_Line *)pShape)->Get_Part(iPart)->Get_Extent().Intersects(m_pGrid->Get_Extent()) )
460 {
461 continue;
462 }
463
464 int iPoint = pShape->Get_Type() == SHAPE_TYPE_Polygon ? 0 : 1;
465
466 TSG_Point A = pShape->Get_Point(0, iPart, iPoint != 0);
467
468 A.x = X_WORLD_TO_GRID(A.x);
469 A.y = Y_WORLD_TO_GRID(A.y);
470
471 for( ; iPoint<pShape->Get_Point_Count(iPart); iPoint++)
472 {
473 TSG_Point B = A; A = pShape->Get_Point(iPoint, iPart);
474
475 A.x = X_WORLD_TO_GRID(A.x);
476 A.y = Y_WORLD_TO_GRID(A.y);
477
478 if( bFat )
479 {
480 Set_Line_Fat(A, B, Value);
481 }
482 else
483 {
484 Set_Line_Thin(A, B, Value);
485 }
486 }
487 }
488 }
489
490 //---------------------------------------------------------
Set_Line_Thin(TSG_Point a,TSG_Point b,double Value)491 void CShapes2Grid::Set_Line_Thin(TSG_Point a, TSG_Point b, double Value)
492 {
493 TSG_Point_Int A; A.x = (int)(a.x += 0.5); A.y = (int)(a.y += 0.5);
494 TSG_Point_Int B; B.x = (int)(b.x += 0.5); B.y = (int)(b.y += 0.5);
495
496 //-----------------------------------------------------
497 if( A.x != B.x || A.y != B.y )
498 {
499 double dx = b.x - a.x;
500 double dy = b.y - a.y;
501
502 if( fabs(dx) > fabs(dy) )
503 {
504 int sig = dx < 0 ? -1 : 1;
505 dx = fabs(dx);
506 dy /= dx;
507
508 for(int ix=0; ix<=dx; ix++, a.x+=sig, a.y+=dy)
509 {
510 Set_Value((int)a.x, (int)a.y, Value);
511 }
512 }
513 else if( fabs(dy) >= fabs(dx) && dy != 0 )
514 {
515 int sig = dy < 0 ? -1 : 1;
516 dy = fabs(dy);
517 dx /= dy;
518
519 for(int iy=0; iy<=dy; iy++, a.x+=dx, a.y+=sig)
520 {
521 Set_Value((int)a.x, (int)a.y, Value);
522 }
523 }
524 }
525 else
526 {
527 Set_Value(A.x, A.y, Value);
528 }
529 }
530
531 //---------------------------------------------------------
Set_Line_Fat(TSG_Point a,TSG_Point b,double Value)532 void CShapes2Grid::Set_Line_Fat(TSG_Point a, TSG_Point b, double Value)
533 {
534 TSG_Point_Int A; A.x = (int)(a.x += 0.5); A.y = (int)(a.y += 0.5);
535 TSG_Point_Int B; B.x = (int)(b.x += 0.5); B.y = (int)(b.y += 0.5);
536
537 Set_Value(A.x, A.y, Value);
538
539 //-----------------------------------------------------
540 if( A.x != B.x || A.y != B.y )
541 {
542 double dx = b.x - a.x;
543 double dy = b.y - a.y;
544
545 a.x = a.x > 0. ? a.x - (int)a.x : 1. + (a.x - (int)a.x);
546 a.y = a.y > 0. ? a.y - (int)a.y : 1. + (a.y - (int)a.y);
547
548 //-------------------------------------------------
549 if( fabs(dx) > fabs(dy) )
550 {
551 int ix = dx > 0. ? 1 : -1;
552 int iy = dy > 0. ? 1 : -1;
553
554 double d, e;
555
556 d = fabs(dy / dx);
557 dx = ix < 0 ? a.x : 1. - a.x;
558 e = iy > 0 ? a.y : 1. - a.y;
559 e += d * dx;
560
561 while( e > 1. )
562 {
563 Set_Value(A.x, A.y += iy, Value); e--;
564 }
565
566 while( A.x != B.x )
567 {
568 Set_Value(A.x += ix, A.y, Value); e += d;
569
570 if( A.x != B.x )
571 {
572 while( e > 1. )
573 {
574 Set_Value(A.x, A.y += iy, Value); e--;
575 }
576 }
577 }
578
579 if( A.y != B.y )
580 {
581 iy = A.y < B.y ? 1 : -1;
582
583 while( A.y != B.y )
584 {
585 Set_Value(A.x, A.y += iy, Value);
586 }
587 }
588 }
589
590 //-------------------------------------------------
591 else // if( fabs(dy) > fabs(dx) )
592 {
593 int ix = dx > 0.0 ? 1 : -1;
594 int iy = dy > 0.0 ? 1 : -1;
595
596 double d, e;
597
598 d = fabs(dx / dy);
599 dy = iy < 0 ? a.y : 1.0 - a.y;
600 e = ix > 0 ? a.x : 1.0 - a.x;
601 e += d * dy;
602
603 while( e > 1. )
604 {
605 Set_Value(A.x += ix, A.y, Value); e--;
606 }
607
608 while( A.y != B.y )
609 {
610 Set_Value(A.x, A.y += iy, Value); e += d;
611
612 if( A.y != B.y )
613 {
614 while( e > 1. )
615 {
616 Set_Value(A.x += ix, A.y, Value); e--;
617 }
618 }
619 }
620
621 if( A.x != B.x )
622 {
623 ix = A.x < B.x ? 1 : -1;
624
625 while( A.x != B.x )
626 {
627 Set_Value(A.x += ix, A.y, Value);
628 }
629 }
630 }
631 }
632 }
633
634
635 ///////////////////////////////////////////////////////////
636 // //
637 ///////////////////////////////////////////////////////////
638
639 //---------------------------------------------------------
Set_Polygon(CSG_Shape * pShape,bool bFat,double Value)640 void CShapes2Grid::Set_Polygon(CSG_Shape *pShape, bool bFat, double Value)
641 {
642 Set_Polygon((CSG_Shape_Polygon *)pShape, Value);
643
644 if( bFat ) // all cells intersected have to be marked
645 {
646 Set_Line(pShape, true, Value); // thick, each cell crossed by polygon boundary will be marked additionally
647 }
648 }
649
650 //---------------------------------------------------------
Set_Polygon(CSG_Shape_Polygon * pPolygon,double Value)651 void CShapes2Grid::Set_Polygon(CSG_Shape_Polygon *pPolygon, double Value)
652 {
653 //-----------------------------------------------------
654 bool *bCrossing = (bool *)SG_Malloc(m_pGrid->Get_NX() * sizeof(bool));
655
656 CSG_Rect Extent(pPolygon->Get_Extent());
657
658 int xStart = (int)((Extent.Get_XMin() - m_pGrid->Get_XMin()) / m_pGrid->Get_Cellsize()) - 1; if( xStart < 0 ) xStart = 0;
659 int xStop = (int)((Extent.Get_XMax() - m_pGrid->Get_XMin()) / m_pGrid->Get_Cellsize()) + 1; if( xStop >= m_pGrid->Get_NX() ) xStop = m_pGrid->Get_NX() - 1;
660
661 TSG_Point A, B;
662
663 A.y = m_pGrid->Get_YMin();
664 A.x = m_pGrid->Get_XMin() - 1.;
665 B.x = m_pGrid->Get_XMax() + 1.;
666
667 //-----------------------------------------------------
668 for(int y=0; y<m_pGrid->Get_NY(); y++, A.y+=m_pGrid->Get_Cellsize())
669 {
670 if( A.y >= Extent.m_rect.yMin && A.y <= Extent.m_rect.yMax )
671 {
672 B.y = A.y;
673
674 memset(bCrossing, 0, m_pGrid->Get_NX() * sizeof(bool));
675
676 for(int iPart=0; iPart<pPolygon->Get_Part_Count(); iPart++)
677 {
678 if( pPolygon->Get_Part(iPart)->Get_Extent().Intersects(m_pGrid->Get_Extent(true)) )
679 {
680 TSG_Point a, b, c;
681
682 b = pPolygon->Get_Point(pPolygon->Get_Point_Count(iPart) - 1, iPart);
683
684 for(int iPoint=0; iPoint<pPolygon->Get_Point_Count(iPart); iPoint++)
685 {
686 a = b;
687 b = pPolygon->Get_Point(iPoint, iPart);
688
689 if( ((a.y <= A.y && A.y < b.y)
690 || (a.y > A.y && A.y >= b.y)) )
691 {
692 SG_Get_Crossing(c, a, b, A, B, false);
693
694 int x = (int)(1.0 + X_WORLD_TO_GRID(c.x));
695
696 if( x < 0 )
697 {
698 x = 0;
699 }
700 else if( x >= m_pGrid->Get_NX() )
701 {
702 continue;
703 }
704
705 bCrossing[x] = !bCrossing[x];
706 }
707 }
708 }
709 }
710
711 //---------------------------------------------
712 bool bFill = false;
713
714 for(int x=xStart; x<=xStop; x++)
715 {
716 if( bCrossing[x] )
717 {
718 bFill = !bFill;
719 }
720
721 if( bFill )
722 {
723 Set_Value(x, y, Value);
724 }
725 }
726 }
727 }
728
729 //-----------------------------------------------------
730 SG_Free(bCrossing);
731 }
732
733
734 ///////////////////////////////////////////////////////////
735 // //
736 // //
737 // //
738 ///////////////////////////////////////////////////////////
739
740 //---------------------------------------------------------
CPolygons2Grid(void)741 CPolygons2Grid::CPolygons2Grid(void)
742 {
743 Set_Name (_TL("Polygons to Grid"));
744
745 Set_Author ("O.Conrad (c) 2018");
746
747 Set_Description (_TW(
748 "Gridding of polygons. If any polygons are selected, only these will be gridded."
749 ));
750
751 //-----------------------------------------------------
752 Parameters.Add_Shapes("",
753 "POLYGONS" , _TL("Polygons"),
754 _TL(""),
755 PARAMETER_INPUT, SHAPE_TYPE_Polygon
756 );
757
758 Parameters.Add_Table_Field("POLYGONS",
759 "FIELD" , _TL("Attribute"),
760 _TL("")
761 );
762
763 Parameters.Add_Choice("",
764 "OUTPUT" , _TL("Output Values"),
765 _TL(""),
766 CSG_String::Format("%s|%s",
767 _TL("index number"),
768 _TL("attribute")
769 ), 2
770 );
771
772 Parameters.Add_Choice("",
773 "MULTIPLE" , _TL("Multiple Polygons"),
774 _TL("Output value for cells that intersect with more than one polygon."),
775 CSG_String::Format("%s|%s|%s",
776 _TL("minimum coverage"),
777 _TL("maximum coverage"),
778 _TL("average proportional to area coverage")
779 ), 1
780 );
781
782 Parameters.Add_Choice("",
783 "GRID_TYPE" , _TL("Data Type"),
784 _TL(""),
785 CSG_String::Format("%s|%s|%s|%s|%s|%s|%s|%s|%s|%s",
786 _TL("1 bit"),
787 _TL("1 byte unsigned integer"),
788 _TL("1 byte signed integer"),
789 _TL("2 byte unsigned integer"),
790 _TL("2 byte signed integer"),
791 _TL("4 byte unsigned integer"),
792 _TL("4 byte signed integer"),
793 _TL("4 byte floating point"),
794 _TL("8 byte floating point"),
795 _TL("same as attribute")
796 ), 9
797 );
798
799 //-----------------------------------------------------
800 m_Grid_Target.Create(&Parameters, false, "", "TARGET_");
801
802 m_Grid_Target.Add_Grid("GRID" , _TL("Grid" ), false);
803 m_Grid_Target.Add_Grid("COVERAGE", _TL("Coverage"), true);
804 }
805
806
807 ///////////////////////////////////////////////////////////
808 // //
809 ///////////////////////////////////////////////////////////
810
811 //---------------------------------------------------------
On_Parameter_Changed(CSG_Parameters * pParameters,CSG_Parameter * pParameter)812 int CPolygons2Grid::On_Parameter_Changed(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
813 {
814 if( pParameter->Cmp_Identifier("POLYGONS") )
815 {
816 m_Grid_Target.Set_User_Defined(pParameters, pParameter->asShapes());
817 }
818
819 m_Grid_Target.On_Parameter_Changed(pParameters, pParameter);
820
821 return( CSG_Tool::On_Parameter_Changed(pParameters, pParameter) );
822 }
823
824 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)825 int CPolygons2Grid::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
826 {
827 if( pParameter->Cmp_Identifier("OUTPUT") )
828 {
829 pParameters->Set_Enabled("FIELD" , pParameter->asInt() == 1);
830 pParameters->Set_Enabled("GRID_TYPE", pParameter->asInt() == 1);
831 }
832
833 m_Grid_Target.On_Parameters_Enable(pParameters, pParameter);
834
835 return( CSG_Tool::On_Parameters_Enable(pParameters, pParameter) );
836 }
837
838
839 ///////////////////////////////////////////////////////////
840 // //
841 ///////////////////////////////////////////////////////////
842
843 //---------------------------------------------------------
Get_Data_Type(int Field)844 TSG_Data_Type CPolygons2Grid::Get_Data_Type(int Field)
845 {
846 CSG_Shapes *pShapes = Parameters("POLYGONS")->asShapes();
847
848 if( Field >= 0 && (Field >= pShapes->Get_Field_Count() || !SG_Data_Type_is_Numeric(pShapes->Get_Field_Type(Field))) )
849 {
850 Field = OUTPUT_INDEX; // index number
851 }
852
853 if( Field <= OUTPUT_INDEX ) // index number
854 {
855 return( pShapes->Get_Count() < 65535 ? SG_DATATYPE_Word : SG_DATATYPE_DWord );
856 }
857
858 // if( Field >= 0 ) // attribute
859 {
860 switch( Parameters("GRID_TYPE")->asInt() )
861 {
862 case 0: return( SG_DATATYPE_Bit );
863 case 1: return( SG_DATATYPE_Byte );
864 case 2: return( SG_DATATYPE_Char );
865 case 3: return( SG_DATATYPE_Word );
866 case 4: return( SG_DATATYPE_Short );
867 case 5: return( SG_DATATYPE_DWord );
868 case 6: return( SG_DATATYPE_Int );
869 case 7: return( SG_DATATYPE_Float );
870 case 8: return( SG_DATATYPE_Double );
871 }
872
873 return( pShapes->Get_Field_Type(Field) );
874 }
875 }
876
877
878 ///////////////////////////////////////////////////////////
879 // //
880 ///////////////////////////////////////////////////////////
881
882 //---------------------------------------------------------
On_Execute(void)883 bool CPolygons2Grid::On_Execute(void)
884 {
885 //-----------------------------------------------------
886 CSG_Shapes *pPolygons = Parameters("POLYGONS")->asShapes();
887
888 m_Multiple = Parameters("MULTIPLE")->asInt();
889
890 //-----------------------------------------------------
891 int Field;
892
893 switch( Parameters("OUTPUT")->asInt() )
894 {
895 case 0: Field = OUTPUT_INDEX ; break; // index number
896 default: Field = Parameters("FIELD")->asInt(); // attribute
897 if( Field < 0 || !SG_Data_Type_is_Numeric(pPolygons->Get_Field_Type(Field)) )
898 {
899 Message_Add(_TL("WARNING: selected attribute is not numeric."));
900 }
901 break;
902 }
903
904 //-----------------------------------------------------
905 if( (m_pGrid = m_Grid_Target.Get_Grid("GRID", Get_Data_Type(Field))) == NULL )
906 {
907 return( false );
908 }
909
910 if( !pPolygons->Get_Extent().Intersects(m_pGrid->Get_Extent()) )
911 {
912 Error_Set(_TL("Polygons' and target grid's extent do not intersect."));
913
914 return( false );
915 }
916
917 if( Field < 0 )
918 {
919 m_pGrid->Set_NoData_Value(0.);
920 }
921
922 m_pGrid->Fmt_Name("%s [%s]", pPolygons->Get_Name(), Field < 0 ? _TL("ID") : pPolygons->Get_Field_Name(Field));
923 m_pGrid->Assign_NoData();
924
925 //-------------------------------------------------
926 CSG_Grid Coverage;
927
928 m_pCoverage = m_Grid_Target.Get_Grid("COVERAGE");
929
930 if( m_pCoverage == NULL )
931 {
932 Coverage.Create(m_pGrid->Get_System());
933
934 m_pCoverage = &Coverage;
935 }
936
937 m_pCoverage->Fmt_Name("%s [%s]", pPolygons->Get_Name(), _TL("Coverage"));
938 m_pCoverage->Set_NoData_Value(0.);
939 m_pCoverage->Assign(0.);
940
941 //-----------------------------------------------------
942 for(int i=0; i<pPolygons->Get_Count() && Set_Progress(i, pPolygons->Get_Count()); i++)
943 {
944 CSG_Shape_Polygon *pPolygon = (CSG_Shape_Polygon *)pPolygons->Get_Shape(i);
945
946 if( pPolygons->Get_Selection_Count() <= 0 || pPolygon->is_Selected() )
947 {
948 if( Field < 0 || !pPolygon->is_NoData(Field) )
949 {
950 if( pPolygon->Intersects(m_pGrid->Get_Extent()) )
951 {
952 Set_Polygon(pPolygon, Field < 0 ? i + 1 : pPolygon->asDouble(Field));
953 }
954 }
955 }
956 }
957
958 //-----------------------------------------------------
959 if( m_Multiple == 2 ) // average
960 {
961 #ifndef _DEBUG
962 #pragma omp parallel for
963 #endif
964 for(sLong i=0; i<m_pGrid->Get_NCells(); i++)
965 {
966 double Area = m_pCoverage->asDouble(i);
967
968 if( Area > 0. )
969 {
970 m_pGrid->Mul_Value(i, 1. / Area);
971 }
972 }
973 }
974
975 //-----------------------------------------------------
976 return( true );
977 }
978
979
980 ///////////////////////////////////////////////////////////
981 // //
982 ///////////////////////////////////////////////////////////
983
984 //---------------------------------------------------------
Set_Value(int x,int y,double Value,double Coverage)985 inline void CPolygons2Grid::Set_Value(int x, int y, double Value, double Coverage)
986 {
987 if( m_pGrid->is_InGrid(x, y, false) )
988 {
989 if( m_pCoverage->asDouble(x, y) <= 0. )
990 {
991 m_pGrid ->Set_Value(x, y, m_Multiple != 2 ? Value : Coverage * Value);
992 m_pCoverage->Set_Value(x, y, Coverage);
993 }
994 else switch( m_Multiple )
995 {
996 case 0: // minimum
997 if( m_pCoverage->asDouble(x, y) > Coverage )
998 {
999 m_pGrid ->Set_Value(x, y, Value );
1000 m_pCoverage->Set_Value(x, y, Coverage);
1001 }
1002 break;
1003
1004 default: // maximum
1005 if( m_pCoverage->asDouble(x, y) < Coverage )
1006 {
1007 m_pGrid ->Set_Value(x, y, Value );
1008 m_pCoverage->Set_Value(x, y, Coverage);
1009 }
1010 break;
1011
1012 case 2: // average
1013 {
1014 m_pGrid ->Add_Value(x, y, Coverage * Value);
1015 m_pCoverage->Add_Value(x, y, Coverage);
1016 }
1017 break;
1018 }
1019 }
1020 }
1021
1022
1023 ///////////////////////////////////////////////////////////
1024 // //
1025 ///////////////////////////////////////////////////////////
1026
1027 //---------------------------------------------------------
Set_Polygon(CSG_Shape_Polygon * pPolygon,double Value)1028 void CPolygons2Grid::Set_Polygon(CSG_Shape_Polygon *pPolygon, double Value)
1029 {
1030 //-----------------------------------------------------
1031 CSG_Grid_System s(m_pGrid->Get_System());
1032
1033 int xA = s.Get_xWorld_to_Grid(pPolygon->Get_Extent().Get_XMin()); if( xA < 0 ) xA = 0;
1034 int xB = s.Get_xWorld_to_Grid(pPolygon->Get_Extent().Get_XMax()); if( xB >= s.Get_NX() ) xB = s.Get_NX() - 1;
1035 int yA = s.Get_yWorld_to_Grid(pPolygon->Get_Extent().Get_YMin()); if( yA < 0 ) yA = 0;
1036 int yB = s.Get_yWorld_to_Grid(pPolygon->Get_Extent().Get_YMax()); if( yB >= s.Get_NY() ) yB = s.Get_NY() - 1;
1037
1038 //-----------------------------------------------------
1039 CSG_Shapes Intersect(SHAPE_TYPE_Polygon);
1040
1041 CSG_Shape_Polygon *pCell = (CSG_Shape_Polygon *)Intersect.Add_Shape();
1042
1043 TSG_Rect r;
1044
1045 r.yMax = s.Get_yGrid_to_World(yA) - 0.5 * s.Get_Cellsize();
1046
1047 for(int y=yA; y<=yB; y++)
1048 {
1049 r.yMin = r.yMax; r.yMax += s.Get_Cellsize();
1050
1051 r.xMax = s.Get_xGrid_to_World(xA) - 0.5 * s.Get_Cellsize();
1052
1053 for(int x=xA; x<=xB; x++)
1054 {
1055 r.xMin = r.xMax; r.xMax += s.Get_Cellsize();
1056
1057 pCell->Add_Point(r.xMin, r.yMin);
1058 pCell->Add_Point(r.xMin, r.yMax);
1059 pCell->Add_Point(r.xMax, r.yMax);
1060 pCell->Add_Point(r.xMax, r.yMin);
1061
1062 if( SG_Polygon_Intersection(pCell, pPolygon) )
1063 {
1064 Set_Value(x, y, Value, pCell->Get_Area());
1065 }
1066
1067 pCell->Del_Parts();
1068 }
1069 }
1070 }
1071
1072
1073 ///////////////////////////////////////////////////////////
1074 // //
1075 // //
1076 // //
1077 ///////////////////////////////////////////////////////////
1078
1079 //---------------------------------------------------------
1080