1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // grid_gridding //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // grid_cell_polygon_coverage.cpp //
14 // //
15 // Copyright (C) 2016 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 "grid_cell_polygon_coverage.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
CGrid_Cell_Polygon_Coverage(void)59 CGrid_Cell_Polygon_Coverage::CGrid_Cell_Polygon_Coverage(void)
60 {
61 Set_Name (_TL("Grid Cell Area Covered by Polygons"));
62
63 Set_Author ("O.Conrad (c) 2016");
64
65 Set_Description (_TW(
66 "This tool calculates for each grid cell of the selected grid system "
67 "the area that is covered by the input polygons layer. "
68 ));
69
70 //-----------------------------------------------------
71 Parameters.Add_Shapes("",
72 "POLYGONS" , _TL("Polygons"),
73 _TL(""),
74 PARAMETER_INPUT, SHAPE_TYPE_Polygon
75 );
76
77 Parameters.Add_Choice("",
78 "METHOD" , _TL("Method"),
79 _TL("Choose cell wise, if you have not many polygons. Polygon wise will show much better performance, if you have many (small) polygons."),
80 CSG_String::Format("%s|%s",
81 _TL("cell wise"),
82 _TL("polygon wise")
83 ), 1
84 );
85
86 Parameters.Add_Choice("",
87 "OUTPUT" , _TL("Output"),
88 _TL(""),
89 CSG_String::Format("%s|%s",
90 _TL("area"),
91 _TL("percentage")
92 ), 1
93 );
94
95 Parameters.Add_Bool("",
96 "SELECTION" , _TL("Only Selected Polygons"),
97 _TL(""),
98 true
99 );
100
101 //-----------------------------------------------------
102 m_Grid_Target.Create(&Parameters, false, "", "TARGET_");
103
104 m_Grid_Target.Add_Grid("AREA", _TL("Area of Coverage"), false);
105 }
106
107
108 ///////////////////////////////////////////////////////////
109 // //
110 ///////////////////////////////////////////////////////////
111
112 //---------------------------------------------------------
On_Parameter_Changed(CSG_Parameters * pParameters,CSG_Parameter * pParameter)113 int CGrid_Cell_Polygon_Coverage::On_Parameter_Changed(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
114 {
115 if( pParameter->Cmp_Identifier("POLYGONS") )
116 {
117 m_Grid_Target.Set_User_Defined(pParameters, pParameter->asShapes());
118 }
119
120 m_Grid_Target.On_Parameter_Changed(pParameters, pParameter);
121
122 return( CSG_Tool::On_Parameter_Changed(pParameters, pParameter) );
123 }
124
125 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)126 int CGrid_Cell_Polygon_Coverage::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
127 {
128 if( pParameter->Cmp_Identifier("POLYGONS") )
129 {
130 pParameters->Set_Enabled("SELECTION", pParameter->asShapes() && pParameter->asShapes()->Get_Selection_Count() > 0);
131 }
132
133 m_Grid_Target.On_Parameters_Enable(pParameters, pParameter);
134
135 return( CSG_Tool::On_Parameters_Enable(pParameters, pParameter) );
136 }
137
138
139 ///////////////////////////////////////////////////////////
140 // //
141 ///////////////////////////////////////////////////////////
142
143 //---------------------------------------------------------
144 #define GET_NPOLYGONS (bSelection ? pPolygons->Get_Selection_Count() : pPolygons->Get_Count())
145 #define GET_POLYGON(i) ((CSG_Shape_Polygon *)(bSelection ? pPolygons->Get_Selection(i) : pPolygons->Get_Shape(i)))
146
147 //---------------------------------------------------------
On_Execute(void)148 bool CGrid_Cell_Polygon_Coverage::On_Execute(void)
149 {
150 CSG_Shapes *pPolygons = Parameters("POLYGONS")->asShapes();
151
152 CSG_Grid *pArea = m_Grid_Target.Get_Grid("AREA");
153
154 if( pPolygons->Get_Count() <= 0 || pArea == NULL || !pPolygons->Get_Extent().Intersects(pArea->Get_Extent()) )
155 {
156 Error_Set(_TL("no spatial intersection between grid system and polygon layer"));
157
158 return( false );
159 }
160
161 //-----------------------------------------------------
162 bool bSelection = pPolygons->Get_Selection_Count() > 0 ? Parameters("SELECTION")->asBool() : false;
163
164 pArea->Fmt_Name("%s [%s]", pPolygons->Get_Name(), _TL("Coverage"));
165 pArea->Set_NoData_Value(0.0);
166
167 DataObject_Add(pArea);
168 DataObject_Set_Colors(pArea, 11, SG_COLORS_RED_GREEN, true);
169
170 //-----------------------------------------------------
171 if( Parameters("METHOD")->asInt() == 0 )
172 {
173 CSG_Grid_System s(pArea->Get_System());
174
175 for(int y=0; y<s.Get_NY() && Set_Progress(y, s.Get_NY()); y++)
176 {
177 double py = s.Get_YMin() + s.Get_Cellsize() * (y - 0.5);
178
179 #ifndef _DEBUG
180 #pragma omp parallel for
181 #endif
182 for(int x=0; x<s.Get_NX(); x++)
183 {
184 double px = s.Get_XMin() + s.Get_Cellsize() * (x - 0.5);
185
186 CSG_Shapes Cell(SHAPE_TYPE_Polygon);
187 CSG_Shape_Polygon *pCell = (CSG_Shape_Polygon *)Cell.Add_Shape();
188
189 pCell->Add_Point(px , py );
190 pCell->Add_Point(px + s.Get_Cellsize(), py );
191 pCell->Add_Point(px + s.Get_Cellsize(), py + s.Get_Cellsize());
192 pCell->Add_Point(px , py + s.Get_Cellsize());
193
194 //---------------------------------------------
195 if( pPolygons->Get_Extent().Intersects(pCell->Get_Extent()) )
196 {
197 for(size_t i=0; pCell->Get_Area() > 0.0 && i<GET_NPOLYGONS; i++)
198 {
199 if( !SG_Polygon_Difference(pCell, GET_POLYGON(i)) ) // completely contained or identical > difference is empty !
200 {
201 pCell->Del_Parts();
202 }
203 }
204 }
205
206 pArea->Set_Value(x, y, s.Get_Cellarea() > pCell->Get_Area() ? s.Get_Cellarea() - pCell->Get_Area() : 0.0);
207 }
208 }
209 }
210
211 //-----------------------------------------------------
212 else
213 {
214 pArea->Assign(0.0);
215
216 for(size_t i=0; i<GET_NPOLYGONS && Set_Progress(i, GET_NPOLYGONS); i++)
217 {
218 Get_Area(GET_POLYGON(i), pArea);
219 }
220 }
221
222 //-----------------------------------------------------
223 if( Parameters("OUTPUT")->asInt() == 1 ) // Percentage
224 {
225 pArea->Multiply(100.0 / pArea->Get_Cellarea());
226 }
227
228 //-----------------------------------------------------
229 return( true );
230 }
231
232
233 ///////////////////////////////////////////////////////////
234 // //
235 ///////////////////////////////////////////////////////////
236
237 //---------------------------------------------------------
Get_Area(CSG_Shape_Polygon * pPolygon,CSG_Grid * pArea)238 bool CGrid_Cell_Polygon_Coverage::Get_Area(CSG_Shape_Polygon *pPolygon, CSG_Grid *pArea)
239 {
240 CSG_Grid_System s(pArea->Get_System());
241
242 int xMin = s.Get_xWorld_to_Grid(pPolygon->Get_Extent().Get_XMin()); if( xMin < 0 ) xMin = 0;
243 int xMax = s.Get_xWorld_to_Grid(pPolygon->Get_Extent().Get_XMax()); if( xMax >= s.Get_NX() ) xMax = s.Get_NX() - 1;
244 int yMin = s.Get_yWorld_to_Grid(pPolygon->Get_Extent().Get_YMin()); if( yMin < 0 ) yMin = 0;
245 int yMax = s.Get_yWorld_to_Grid(pPolygon->Get_Extent().Get_YMax()); if( yMax >= s.Get_NY() ) yMax = s.Get_NY() - 1;
246
247 double d = 0.5 * s.Get_Cellsize();
248
249 #pragma omp parallel for
250 for(int y=yMin; y<=yMax; y++)
251 {
252 CSG_Shapes Cell(SHAPE_TYPE_Polygon); CSG_Shape_Polygon *pCell = (CSG_Shape_Polygon *)Cell.Add_Shape();
253
254 TSG_Point p; p.y = s.Get_YMin() + s.Get_Cellsize() * y;
255
256 p.x = s.Get_XMin() + s.Get_Cellsize() * xMin;
257
258 for(int x=xMin; x<=xMax; x++, p.x+=s.Get_Cellsize())
259 {
260 pCell->Add_Point(p.x - d, p.y - d);
261 pCell->Add_Point(p.x - d, p.y + d);
262 pCell->Add_Point(p.x + d, p.y + d);
263 pCell->Add_Point(p.x + d, p.y - d);
264
265 if( SG_Polygon_Intersection(pCell, pPolygon) && pCell->Get_Area() > 0.0 )
266 {
267 pArea->Add_Value(x, y, pCell->Get_Area());
268 }
269
270 pCell->Del_Parts();
271 }
272 }
273
274 return( true );
275 }
276
277
278 ///////////////////////////////////////////////////////////
279 // //
280 // //
281 // //
282 ///////////////////////////////////////////////////////////
283
284 //---------------------------------------------------------
CPolygonCategories2Grid(void)285 CPolygonCategories2Grid::CPolygonCategories2Grid(void)
286 {
287 Set_Name (_TL("Polygon Categories to Grid"));
288
289 Set_Author ("O.Conrad (c) 2018");
290
291 Set_Description (_TW(
292 "This tool has been designed to rasterize polygons representing "
293 "categories and selects that category, which has maximum coverage "
294 "of a cell. The advantage using this tool (instead the more simple "
295 "'Shapes to Grid' or 'Polygons to Grid' tools) is that it summarizes "
296 "all polygon coverages belonging to the same category. "
297 ));
298
299 //-----------------------------------------------------
300 Parameters.Add_Shapes("",
301 "POLYGONS" , _TL("Polygons"),
302 _TL(""),
303 PARAMETER_INPUT, SHAPE_TYPE_Polygon
304 );
305
306 Parameters.Add_Table_Field("POLYGONS",
307 "FIELD" , _TL("Category"),
308 _TL("The attribute field that specifies the category a polygon belongs to.")
309 );
310
311 Parameters.Add_Choice("",
312 "METHOD" , _TL("Method"),
313 _TL("Choose cell wise, if you have not many polygons. Polygon wise will show much better performance, if you have many (small) polygons."),
314 CSG_String::Format("%s|%s",
315 _TL("cell wise"),
316 _TL("polygon wise")
317 ), 1
318 );
319
320 Parameters.Add_Choice("",
321 "MULTIPLE" , _TL("Multiple Polygons"),
322 _TL("Output value for cells that intersect wiht more than one polygon."),
323 CSG_String::Format("%s|%s",
324 _TL("minimum coverage"),
325 _TL("maximum coverage")
326 ), 1
327 );
328
329 Parameters.Add_Table("",
330 "CLASSES" , _TL("Classification"),
331 _TL("Classification look-up table for interpretation of resulting grid cell values."),
332 PARAMETER_OUTPUT_OPTIONAL
333 );
334
335 //-----------------------------------------------------
336 m_Grid_Target.Create(&Parameters, false, "", "TARGET_");
337
338 m_Grid_Target.Add_Grid("CATEGORY", _TL("Category"), false);
339 m_Grid_Target.Add_Grid("COVERAGE", _TL("Coverage"), true);
340 }
341
342
343 ///////////////////////////////////////////////////////////
344 // //
345 ///////////////////////////////////////////////////////////
346
347 //---------------------------------------------------------
On_Parameter_Changed(CSG_Parameters * pParameters,CSG_Parameter * pParameter)348 int CPolygonCategories2Grid::On_Parameter_Changed(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
349 {
350 if( pParameter->Cmp_Identifier("POLYGONS") )
351 {
352 m_Grid_Target.Set_User_Defined(pParameters, pParameter->asShapes());
353 }
354
355 m_Grid_Target.On_Parameter_Changed(pParameters, pParameter);
356
357 return( CSG_Tool::On_Parameter_Changed(pParameters, pParameter) );
358 }
359
360 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)361 int CPolygonCategories2Grid::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
362 {
363 m_Grid_Target.On_Parameters_Enable(pParameters, pParameter);
364
365 return( CSG_Tool::On_Parameters_Enable(pParameters, pParameter) );
366 }
367
368
369 ///////////////////////////////////////////////////////////
370 // //
371 ///////////////////////////////////////////////////////////
372
373 //---------------------------------------------------------
On_Execute(void)374 bool CPolygonCategories2Grid::On_Execute(void)
375 {
376 CSG_Shapes *pPolygons = Parameters("POLYGONS")->asShapes();
377
378 int Field = Parameters("FIELD")->asInt();
379
380 bool bNumber = SG_Data_Type_is_Numeric(pPolygons->Get_Field_Type(Field));
381
382 CSG_Grid *pCategory = m_Grid_Target.Get_Grid("CATEGORY", bNumber ? pPolygons->Get_Field_Type(Field) : SG_DATATYPE_Int);
383
384 if( pPolygons->Get_Count() <= 0 || pCategory == NULL || !pPolygons->Get_Extent().Intersects(pCategory->Get_Extent()) )
385 {
386 Error_Set(_TL("no spatial intersection between grid system and polygon layer"));
387
388 return( false );
389 }
390
391 //-----------------------------------------------------
392 pCategory->Fmt_Name("%s [%s]", pPolygons->Get_Name(), pPolygons->Get_Field_Name(Field));
393 pCategory->Assign_NoData();
394
395 //-----------------------------------------------------
396 if( !pPolygons->Set_Index(Field, TABLE_INDEX_Ascending) )
397 {
398 Error_Set(_TL("index creation failed"));
399
400 return( false );
401 }
402
403 //-----------------------------------------------------
404 CSG_Grid Coverage, *pCoverage = m_Grid_Target.Get_Grid("COVERAGE");
405
406 if( pCoverage == NULL )
407 {
408 Coverage.Create(pCategory->Get_System());
409
410 pCoverage = &Coverage;
411 }
412
413 pCoverage->Fmt_Name("%s [%s]", pPolygons->Get_Name(), _TL("Coverage"));
414 pCoverage->Set_NoData_Value(0.);
415 pCoverage->Assign(0.);
416
417 //-----------------------------------------------------
418 CSG_Table Classes;
419
420 Classes.Add_Field("COLOR" , SG_DATATYPE_Color );
421 Classes.Add_Field("NAME" , SG_DATATYPE_String);
422 Classes.Add_Field("DESCRIPTION", SG_DATATYPE_String);
423 Classes.Add_Field("MINIMUM" , SG_DATATYPE_Double);
424 Classes.Add_Field("MAXIMUM" , SG_DATATYPE_Double);
425
426 //-----------------------------------------------------
427 CSG_String Category;
428
429 pPolygons->Select(); // clear selection
430
431 for(int i=0; i<pPolygons->Get_Count() && Set_Progress(i, pPolygons->Get_Count()); i++)
432 {
433 CSG_Shape *pPolygon = pPolygons->Get_Shape_byIndex(i);
434
435 if( Category.Cmp(pPolygon->asString(Field)) ) // new category
436 {
437 Set_Category(pPolygons, pCategory, pCoverage, Classes, Category, bNumber); // also clears selection
438
439 Category = pPolygon->asString(Field);
440 }
441
442 pPolygons->Select(pPolygon, true); // adds polygon to selection
443 }
444
445 Set_Category(pPolygons, pCategory, pCoverage, Classes, Category, bNumber);
446
447 //-----------------------------------------------------
448 DataObject_Add (pCategory);
449 DataObject_Update(pCategory);
450
451 CSG_Parameter *pLUT = DataObject_Get_Parameter(pCategory, "LUT");
452
453 if( pLUT && pLUT->asTable() && pLUT->asTable()->Create(Classes) )
454 {
455 DataObject_Set_Parameter(pCategory, pLUT);
456 DataObject_Set_Parameter(pCategory, "COLORS_TYPE", 1); // Color Classification Type: Lookup Table
457 }
458
459 if( Parameters("CLASSES")->asTable() )
460 {
461 Classes.Del_Field(4); // MAXIMUM
462 Classes.Del_Field(2); // DESCRIPTION
463 Classes.Del_Field(0); // COLOR
464
465 Classes.Set_Field_Name(0, _TL("Category"));
466 Classes.Set_Field_Name(1, _TL("Value" ));
467
468 Parameters("CLASSES")->asTable()->Create(Classes);
469 }
470
471 //-----------------------------------------------------
472 return( true );
473 }
474
475
476 ///////////////////////////////////////////////////////////
477 // //
478 ///////////////////////////////////////////////////////////
479
480 //---------------------------------------------------------
Set_Category(CSG_Shapes * pPolygons,CSG_Grid * pCategory,CSG_Grid * pCoverage,CSG_Table & Classes,const CSG_String & Category,bool bNumber)481 bool CPolygonCategories2Grid::Set_Category(CSG_Shapes *pPolygons, CSG_Grid *pCategory, CSG_Grid *pCoverage, CSG_Table &Classes, const CSG_String &Category, bool bNumber)
482 {
483 if( pPolygons->Get_Selection_Count() <= 0 )
484 {
485 return( false );
486 }
487
488 //-----------------------------------------------------
489 CSG_Grid Coverage(pCoverage->Get_System());
490
491 CGrid_Cell_Polygon_Coverage Get_Coverage;
492
493 Get_Coverage.Set_Manager(NULL);
494
495 Get_Coverage.Set_Parameter("POLYGONS" , pPolygons);
496 Get_Coverage.Set_Parameter("METHOD" , Parameters("METHOD"));
497 Get_Coverage.Set_Parameter("OUTPUT" , 0); // area (not percentage)
498 Get_Coverage.Set_Parameter("TARGET_DEFINITION", 1); // grid or grid system
499 Get_Coverage.Set_Parameter("AREA" , &Coverage);
500
501 SG_UI_ProgressAndMsg_Lock(true);
502
503 if( !Get_Coverage.Execute() )
504 {
505 SG_UI_ProgressAndMsg_Lock(false);
506
507 pPolygons->Select(); // clear selection
508
509 return( false );
510 }
511
512 SG_UI_ProgressAndMsg_Lock(false);
513
514 pPolygons->Select(); // clear selection
515
516 //-----------------------------------------------------
517 CSG_Table_Record &Class = *Classes.Add_Record();
518
519 double ClassID = bNumber ? Category.asDouble() : Classes.Get_Count();
520
521 Class.Set_Value(0, SG_Color_Get_Random());
522 Class.Set_Value(1, Category);
523 Class.Set_Value(3, ClassID);
524 Class.Set_Value(4, ClassID);
525
526 int Multiple = Parameters("MULTIPLE")->asInt();
527
528 #pragma omp parallel for
529 for(sLong i=0; i<pCategory->Get_NCells(); i++)
530 {
531 switch( Multiple )
532 {
533 case 0: // minimum
534 if( pCoverage->asDouble(i) <= 0. || Coverage.asDouble(i) < pCoverage->asDouble(i) )
535 {
536 pCategory->Set_Value(i, ClassID);
537 pCoverage->Set_Value(i, Coverage.asDouble(i));
538 }
539 break;
540
541 default: // maximum
542 if( Coverage.asDouble(i) > pCoverage->asDouble(i) )
543 {
544 pCategory->Set_Value(i, ClassID);
545 pCoverage->Set_Value(i, Coverage.asDouble(i));
546 }
547 break;
548 }
549 }
550
551 //-----------------------------------------------------
552 return( true );
553 }
554
555
556 ///////////////////////////////////////////////////////////
557 // //
558 // //
559 // //
560 ///////////////////////////////////////////////////////////
561
562 //---------------------------------------------------------
563