1""" 2@package iscatt.iscatt_core 3 4@brief Non GUI functions. 5 6Classes: 7 - iscatt_core::Core 8 - iscatt_core::CatRastUpdater 9 - iscatt_core::AnalyzedData 10 - iscatt_core::ScattPlotsCondsData 11 - iscatt_core::ScattPlotsData 12 13(C) 2013 by the GRASS Development Team 14 15This program is free software under the GNU General Public License 16(>=v2). Read the file COPYING that comes with GRASS for details. 17 18@author Stepan Turek <stepan.turek seznam.cz> (mentor: Martin Landa) 19""" 20import os 21import sys 22import six 23 24import numpy as np 25 26# used iclass perimeters algorithm instead of convolve2d 27#from matplotlib.path import Path 28#from scipy.signal import convolve2d 29 30from math import sqrt, ceil, floor 31from copy import deepcopy 32 33from core.gcmd import GException, GError, RunCommand 34 35import grass.script as grass 36 37from iscatt.core_c import CreateCatRast, ComputeScatts, UpdateCatRast, \ 38 Rasterize, SC_SCATT_DATA, SC_SCATT_CONDITIONS 39 40MAX_SCATT_SIZE = 4100 * 4100 41WARN_SCATT_SIZE = 2000 * 2000 42MAX_NCELLS = 65536 * 65536 43WARN_NCELLS = 12000 * 12000 44 45 46class Core: 47 """Represents scatter plot backend. 48 """ 49 50 def __init__(self): 51 52 self.an_data = AnalyzedData() 53 54 self.scatts_dt = ScattPlotsData(self.an_data) 55 self.scatt_conds_dt = ScattPlotsCondsData(self.an_data) 56 57 self.cat_rast_updater = CatRastUpdater( 58 self.scatts_dt, self.an_data, self) 59 60 def SetData(self, bands): 61 """Set bands for analysis. 62 """ 63 ret = self.an_data.Create(bands) 64 if not ret: 65 return False 66 67 n_bands = len(self.GetBands()) 68 69 self.scatts_dt.Create(n_bands) 70 self.scatt_conds_dt.Create(n_bands) 71 72 return True 73 74 def AddCategory(self, cat_id): 75 self.scatts_dt.AddCategory(cat_id) 76 return self.scatt_conds_dt.AddCategory(cat_id) 77 78 def DeleteCategory(self, cat_id): 79 self.scatts_dt.DeleteCategory(cat_id) 80 self.scatt_conds_dt.DeleteCategory(cat_id) 81 82 def CleanUp(self): 83 self.scatts_dt.CleanUp() 84 self.scatt_conds_dt.CleanUp() 85 86 def GetBands(self): 87 return self.an_data.GetBands() 88 89 def GetScattsData(self): 90 return self.scatts_dt, self.scatt_conds_dt 91 92 def GetRegion(self): 93 return self.an_data.GetRegion() 94 95 def GetCatRast(self, cat_id): 96 return self.scatts_dt.GetCatRast(cat_id) 97 98 def AddScattPlots(self, scatt_ids): 99 100 for s_id in scatt_ids: 101 self.scatts_dt.AddScattPlot(scatt_id=s_id) 102 103 cats_ids = self.scatts_dt.GetCategories() 104 self.ComputeCatsScatts(cats_ids) 105 106 def SetEditCatData(self, cat_id, scatt_id, bbox, value): 107 108 if cat_id not in self.scatts_dt.GetCategories(): 109 raise GException(_("Select category for editing.")) 110 111 if self.scatt_conds_dt.AddScattPlot(cat_id, scatt_id) < 0: 112 return None 113 114 arr = self.scatt_conds_dt.GetValuesArr(cat_id, scatt_id) 115 116 for k, v in six.iteritems(bbox): 117 bbox[k] = self._validExtend(v) 118 119 arr[bbox['btm_y']: bbox['up_y'], bbox['btm_x']: bbox['up_x']] = value 120 121 self.ComputeCatsScatts([cat_id]) 122 return cat_id 123 124 def ComputeCatsScatts(self, cats_ids): 125 126 requested_dt = {} 127 requested_dt_conds = {} 128 129 for c in cats_ids: 130 requested_dt_conds[c] = self.scatt_conds_dt.GetCatScatts(c) 131 requested_dt[c] = self.scatts_dt.GetCatScatts(c) 132 133 scatt_conds = self.scatt_conds_dt.GetData(requested_dt_conds) 134 scatts = self.scatts_dt.GetData(requested_dt) 135 136 bands = self.an_data.GetBands() 137 138 cats_rasts = self.scatts_dt.GetCatsRasts() 139 cats_rasts_conds = self.scatts_dt.GetCatsRastsConds() 140 141 returncode, scatts = ComputeScatts(self.an_data.GetRegion(), 142 scatt_conds, 143 bands, 144 len(self.GetBands()), 145 scatts, 146 cats_rasts_conds, 147 cats_rasts) 148 149 if returncode < 0: 150 GException(_("Computing of scatter plots failed.")) 151 152 def CatRastUpdater(self): 153 return self.cat_rast_updater 154 155 def UpdateCategoryWithPolygons(self, cat_id, scatts_pols, value): 156 if cat_id not in self.scatts_dt.GetCategories(): 157 raise GException(_("Select category for editing.")) 158 159 for scatt_id, coords in six.iteritems(scatts_pols): 160 161 if self.scatt_conds_dt.AddScattPlot(cat_id, scatt_id) < 0: 162 return False 163 164 b1, b2 = idScattToidBands(scatt_id, len(self.an_data.GetBands())) 165 b = self.scatts_dt.GetBandsInfo(scatt_id) 166 167 region = {} 168 region['s'] = b['b2']['min'] - 0.5 169 region['n'] = b['b2']['max'] + 0.5 170 171 region['w'] = b['b1']['min'] - 0.5 172 region['e'] = b['b1']['max'] + 0.5 173 174 arr = self.scatt_conds_dt.GetValuesArr(cat_id, scatt_id) 175 arr = Rasterize(polygon=coords, 176 rast=arr, 177 region=region, 178 value=value) 179 180 # previous way of rasterization / used scipy 181 # raster_pol = RasterizePolygon(coords, b['b1']['range'], b['b1']['min'], 182 # b['b2']['range'], b['b2']['min']) 183 184 #raster_ind = np.where(raster_pol > 0) 185 #arr = self.scatt_conds_dt.GetValuesArr(cat_id, scatt_id) 186 187 #arr[raster_ind] = value 188 # arr.flush() 189 190 self.ComputeCatsScatts([cat_id]) 191 return cat_id 192 193 def ExportCatRast(self, cat_id, rast_name): 194 195 cat_rast = self.scatts_dt.GetCatRast(cat_id) 196 if not cat_rast: 197 return 1 198 199 return RunCommand("g.copy", 200 raster=cat_rast + ',' + rast_name, 201 getErrorMsg=True, 202 overwrite=True) 203 204 def _validExtend(self, val): 205 # TODO do it general 206 if val > 255: 207 val = 255 208 elif val < 0: 209 val = 0 210 211 return val 212 213 214class CatRastUpdater: 215 """Update backend data structures according to selected areas in mapwindow. 216 """ 217 218 def __init__(self, scatts_dt, an_data, core): 219 self.scatts_dt = scatts_dt 220 self.an_data = an_data # TODO may be confusing 221 self.core = core 222 self.vectMap = None 223 224 def SetVectMap(self, vectMap): 225 self.vectMap = vectMap 226 227 def SyncWithMap(self): 228 # TODO possible optimization - bbox only of vertex and its two 229 # neighbours 230 231 region = self.an_data.GetRegion() 232 233 bbox = {} 234 bbox['maxx'] = region['e'] 235 bbox['minx'] = region['w'] 236 bbox['maxy'] = region['n'] 237 bbox['miny'] = region['s'] 238 239 updated_cats = [] 240 241 for cat_id in self.scatts_dt.GetCategories(): 242 if cat_id == 0: 243 continue 244 245 cat = [{1: [cat_id]}] 246 self._updateCatRast(bbox, cat, updated_cats) 247 248 return updated_cats 249 250 def EditedFeature(self, new_bboxs, new_areas_cats, 251 old_bboxs, old_areas_cats): 252 # TODO possible optimization - bbox only of vertex and its two 253 # neighbours 254 255 bboxs = old_bboxs + new_bboxs 256 areas_cats = old_areas_cats + new_areas_cats 257 258 updated_cats = [] 259 260 for i in range(len(areas_cats)): 261 self._updateCatRast(bboxs[i], areas_cats[i], updated_cats) 262 263 return updated_cats 264 265 def _updateCatRast(self, bbox, areas_cats, updated_cats): 266 267 rasterized_cats = [] 268 for c in range(len(areas_cats)): 269 270 if not areas_cats[c]: 271 continue 272 273 layer = list(areas_cats[c])[0] 274 cat = areas_cats[c][layer][0] 275 276 if cat in rasterized_cats: 277 continue 278 279 rasterized_cats.append(cat) 280 updated_cats.append(cat) 281 282 grass_region = self._create_grass_region_env(bbox) 283 284 # TODO hack check if raster exists? 285 patch_rast = "temp_scatt_patch_%d" % (os.getpid()) 286 self._rasterize(grass_region, layer, cat, patch_rast) 287 288 region = self.an_data.GetRegion() 289 ret = UpdateCatRast( 290 patch_rast, 291 region, 292 self.scatts_dt.GetCatRastCond(cat)) 293 if ret < 0: 294 GException( 295 _("Patching category raster conditions file failed.")) 296 RunCommand("g.remove", flags='f', type='raster', name=patch_rast) 297 298 def _rasterize(self, grass_region, layer, cat, out_rast): 299 300 # TODO different thread may be problem when user edits map 301 environs = os.environ.copy() 302 environs['GRASS_VECTOR_TEMPORARY'] = '1' 303 304 ret, text, msg = RunCommand("v.category", 305 input=self.vectMap, 306 getErrorMsg=True, 307 option='report', 308 read=True, 309 env=environs) 310 311 ret, text, msg = RunCommand("v.build", 312 map=self.vectMap, 313 getErrorMsg=True, 314 read=True, 315 env=environs) 316 317 if ret != 0: 318 GException(_("v.build failed:\n%s" % msg)) 319 320 environs = os.environ.copy() 321 environs["GRASS_REGION"] = grass_region["GRASS_REGION"] 322 environs['GRASS_VECTOR_TEMPORARY'] = '1' 323 324 ret, text, msg = RunCommand("v.to.rast", 325 input=self.vectMap, 326 use="cat", 327 layer=str(layer), 328 cat=str(cat), 329 output=out_rast, 330 getErrorMsg=True, 331 read=True, 332 overwrite=True, 333 env=environs) 334 335 if ret != 0: 336 GException(_("v.to.rast failed:\n%s" % msg)) 337 338 def _create_grass_region_env(self, bbox): 339 340 r = self.an_data.GetRegion() 341 new_r = {} 342 343 if bbox["maxy"] <= r["s"]: 344 return 0 345 elif bbox["maxy"] >= r["n"]: 346 new_r["n"] = bbox["maxy"] 347 else: 348 new_r["n"] = ceil( 349 (bbox["maxy"] - r["s"]) / r["nsres"]) * r["nsres"] + r["s"] 350 351 if bbox["miny"] >= r["n"]: 352 return 0 353 elif bbox["miny"] <= r["s"]: 354 new_r["s"] = bbox["miny"] 355 else: 356 new_r["s"] = floor( 357 (bbox["miny"] - r["s"]) / r["nsres"]) * r["nsres"] + r["s"] 358 359 if bbox["maxx"] <= r["w"]: 360 return 0 361 elif bbox["maxx"] >= r["e"]: 362 new_r["e"] = bbox["maxx"] 363 else: 364 new_r["e"] = ceil( 365 (bbox["maxx"] - r["w"]) / r["ewres"]) * r["ewres"] + r["w"] 366 367 if bbox["minx"] >= r["e"]: 368 return 0 369 elif bbox["minx"] <= r["w"]: 370 new_r["w"] = bbox["minx"] 371 else: 372 new_r["w"] = floor( 373 (bbox["minx"] - r["w"]) / r["ewres"]) * r["ewres"] + r["w"] 374 375 # TODO check regions resolution 376 new_r["nsres"] = r["nsres"] 377 new_r["ewres"] = r["ewres"] 378 379 return {"GRASS_REGION": grass.region_env(**new_r)} 380 381 382class AnalyzedData: 383 """Represents analyzed data (bands, region). 384 """ 385 386 def __init__(self): 387 388 self.bands = [] 389 self.bands_info = {} 390 391 self.region = None 392 393 def GetRegion(self): 394 return self.region 395 396 def Create(self, bands): 397 self.bands = bands[:] 398 self.region = None 399 400 self.region = GetRegion() 401 if self.region["rows"] * self.region["cols"] > MAX_NCELLS: 402 GException("too big region") 403 404 self.bands_info = {} 405 406 for b in self.bands[:]: 407 i = GetRasterInfo(b) 408 409 if i is None: 410 GException("raster %s is not CELL type" % (b)) 411 412 self.bands_info[b] = i 413 # TODO check size of raster 414 415 return True 416 417 def GetBands(self): 418 return self.bands 419 420 def GetBandInfo(self, band_id): 421 band = self.bands[band_id] 422 return self.bands_info[band] 423 424 425class ScattPlotsCondsData: 426 """Data structure for selected areas in scatter plot(conditions). 427 """ 428 429 def __init__(self, an_data): 430 431 self.an_data = an_data 432 433 # TODO 434 self.max_n_cats = 10 435 436 self.dtype = 'uint8' 437 self.type = 1 438 self.CleanUp() 439 440 def CleanUp(self): 441 442 self.cats = {} 443 444 self.n_scatts = -1 445 self.n_bands = -1 446 447 for cat_id in self.cats.keys(): 448 self.DeleteCategory(cat_id) 449 450 def Create(self, n_bands): 451 452 self.CleanUp() 453 454 self.n_scatts = (n_bands - 1) * n_bands / 2 455 self.n_bands = n_bands 456 457 self.AddCategory(cat_id=0) 458 459 def AddCategory(self, cat_id): 460 461 if cat_id not in self.cats.keys(): 462 self.cats[cat_id] = {} 463 return cat_id 464 return -1 465 466 def DeleteCategory(self, cat_id): 467 468 if cat_id not in self.cats.keys(): 469 return False 470 471 for scatt in six.itervalues(self.cats[cat_id]): 472 grass.try_remove(scatt['np_vals']) 473 del scatt['np_vals'] 474 475 del self.cats[cat_id] 476 477 return True 478 479 def GetCategories(self): 480 return self.cats.keys() 481 482 def GetCatScatts(self, cat_id): 483 484 if cat_id not in self.cats: 485 return False 486 487 return self.cats[cat_id].keys() 488 489 def AddScattPlot(self, cat_id, scatt_id): 490 491 if cat_id not in self.cats: 492 return -1 493 494 if scatt_id in self.cats[cat_id]: 495 return 0 496 497 b_i = self.GetBandsInfo(scatt_id) 498 499 shape = ( 500 b_i['b2']['max'] - 501 b_i['b2']['min'] + 502 1, 503 b_i['b1']['max'] - 504 b_i['b1']['min'] + 505 1) 506 507 np_vals = np.memmap( 508 grass.tempfile(), 509 dtype=self.dtype, 510 mode='w+', 511 shape=shape) 512 513 self.cats[cat_id][scatt_id] = {'np_vals': np_vals} 514 515 return 1 516 517 def GetBandsInfo(self, scatt_id): 518 b1, b2 = idScattToidBands(scatt_id, len(self.an_data.GetBands())) 519 520 b1_info = self.an_data.GetBandInfo(b1) 521 b2_info = self.an_data.GetBandInfo(b2) 522 523 bands_info = {'b1': b1_info, 524 'b2': b2_info} 525 526 return bands_info 527 528 def DeleScattPlot(self, cat_id, scatt_id): 529 530 if cat_id not in self.cats: 531 return False 532 533 if scatt_id not in self.cats[cat_id]: 534 return False 535 536 del self.cats[cat_id][scatt_id] 537 return True 538 539 def GetValuesArr(self, cat_id, scatt_id): 540 541 if cat_id not in self.cats: 542 return None 543 544 if scatt_id not in self.cats[cat_id]: 545 return None 546 547 return self.cats[cat_id][scatt_id]['np_vals'] 548 549 def GetData(self, requested_dt): 550 551 cats = {} 552 for cat_id, scatt_ids in six.iteritems(requested_dt): 553 if cat_id not in cats: 554 cats[cat_id] = {} 555 for scatt_id in scatt_ids: 556 # if key is missing condition is always True (full scatter plor 557 # is computed) 558 if scatt_id in self.cats[cat_id]: 559 cats[cat_id][scatt_id] = { 560 'np_vals': self.cats[cat_id][scatt_id]['np_vals'], 561 'bands_info': self.GetBandsInfo(scatt_id)} 562 563 return cats 564 565 def SetData(self, cats): 566 567 for cat_id, scatt_ids in six.iteritems(cats): 568 for scatt_id in scatt_ids: 569 # if key is missing condition is always True (full scatter plor 570 # is computed) 571 if scatt_id in self.cats[cat_id]: 572 self.cats[cat_id][scatt_id]['np_vals'] = cats[ 573 cat_id][scatt_id]['np_vals'] 574 575 def GetScatt(self, scatt_id, cats_ids=None): 576 scatts = {} 577 for cat_id in six.iterkeys(self.cats): 578 if cats_ids and cat_id not in cats_ids: 579 continue 580 if scatt_id not in self.cats[cat_id]: 581 continue 582 583 scatts[cat_id] = {'np_vals': self.cats[cat_id][scatt_id][ 584 'np_vals'], 'bands_info': self.GetBandsInfo(scatt_id)} 585 return scatts 586 587 588class ScattPlotsData(ScattPlotsCondsData): 589 """Data structure for computed points (classes) in scatter plots.\ 590 """ 591 592 def __init__(self, an_data): 593 594 self.cats_rasts = {} 595 self.cats_rasts_conds = {} 596 self.scatts_ids = [] 597 598 ScattPlotsCondsData.__init__(self, an_data) 599 600 self.dtype = 'uint32' 601 602 # TODO 603 self.type = 0 604 605 def AddCategory(self, cat_id): 606 cat_id = ScattPlotsCondsData.AddCategory(self, cat_id) 607 if cat_id < 0: 608 return cat_id 609 610 for scatt_id in self.scatts_ids: 611 ScattPlotsCondsData.AddScattPlot(self, cat_id, scatt_id) 612 613 if cat_id == 0: 614 self.cats_rasts_conds[cat_id] = None 615 self.cats_rasts[cat_id] = None 616 else: 617 self.cats_rasts_conds[cat_id] = grass.tempfile() 618 self.cats_rasts[cat_id] = "temp_cat_rast_%d_%d" % ( 619 cat_id, os.getpid()) 620 region = self.an_data.GetRegion() 621 CreateCatRast(region, self.cats_rasts_conds[cat_id]) 622 623 return cat_id 624 625 def DeleteCategory(self, cat_id): 626 627 ScattPlotsCondsData.DeleteCategory(self, cat_id) 628 629 grass.try_remove(self.cats_rasts_conds[cat_id]) 630 del self.cats_rasts_conds[cat_id] 631 632 RunCommand("g.remove", flags='f', type='raster', 633 name=self.cats_rasts[cat_id]) 634 del self.cats_rasts[cat_id] 635 636 return True 637 638 def AddScattPlot(self, scatt_id): 639 640 if scatt_id in self.scatts_ids: 641 return False 642 643 self.scatts_ids.append(scatt_id) 644 for cat_id in six.iterkeys(self.cats): 645 ScattPlotsCondsData.AddScattPlot(self, cat_id, scatt_id) 646 self.cats[cat_id][scatt_id]['ellipse'] = None 647 648 return True 649 650 def DeleteScatterPlot(self, scatt_id): 651 652 if scatt_id not in self.scatts_ids: 653 return False 654 655 self.scatts_ids.remove(scatt_id) 656 657 for cat_id in six.iterkeys(self.cats): 658 ScattPlotsCondsData.DeleteScattPlot(self, cat_id, scatt_id) 659 660 return True 661 662 def GetEllipses(self, scatt_id, styles): 663 if scatt_id not in self.scatts_ids: 664 return False 665 666 scatts = {} 667 for cat_id in six.iterkeys(self.cats): 668 if cat_id == 0: 669 continue 670 nstd = styles[cat_id]['nstd'] 671 scatts[cat_id] = self._getEllipse(cat_id, scatt_id, nstd) 672 673 return scatts 674 675 def _getEllipse(self, cat_id, scatt_id, nstd): 676 # Joe Kington 677 # http://stackoverflow.com/questions/12301071/multidimensional-confidence-intervals 678 679 data = np.copy(self.cats[cat_id][scatt_id]['np_vals']) 680 681 b = self.GetBandsInfo(scatt_id) 682 sel_pts = np.where(data > 0) 683 684 x = sel_pts[1] 685 y = sel_pts[0] 686 687 flatten_data = data.reshape([-1]) 688 flatten_sel_pts = np.nonzero(flatten_data) 689 weights = flatten_data[flatten_sel_pts] 690 if len(weights) == 0: 691 return None 692 693 x_avg = np.average(x, weights=weights) 694 y_avg = np.average(y, weights=weights) 695 pos = np.array([x_avg + b['b1']['min'], y_avg + b['b2']['min']]) 696 697 x_diff = (x - x_avg) 698 y_diff = (y - y_avg) 699 700 x_diff = (x - x_avg) 701 y_diff = (y - y_avg) 702 703 diffs = x_diff * y_diff.T 704 cov = np.dot(diffs, weights) / (np.sum(weights) - 1) 705 706 diffs = x_diff * x_diff.T 707 var_x = np.dot(diffs, weights) / (np.sum(weights) - 1) 708 709 diffs = y_diff * y_diff.T 710 var_y = np.dot(diffs, weights) / (np.sum(weights) - 1) 711 712 cov = np.array([[var_x, cov], [cov, var_y]]) 713 714 def eigsorted(cov): 715 vals, vecs = np.linalg.eigh(cov) 716 order = vals.argsort()[::-1] 717 return vals[order], vecs[:, order] 718 719 vals, vecs = eigsorted(cov) 720 theta = np.degrees(np.arctan2(*vecs[:, 0][::-1])) 721 722 # Width and height are "full" widths, not radius 723 width, height = 2 * nstd * np.sqrt(vals) 724 725 ellipse = {'pos': pos, 726 'width': width, 727 'height': height, 728 'theta': theta} 729 730 del data 731 del flatten_data 732 del flatten_sel_pts 733 del weights 734 del sel_pts 735 return ellipse 736 737 def CleanUp(self): 738 739 ScattPlotsCondsData.CleanUp(self) 740 for tmp in six.itervalues(self.cats_rasts_conds): 741 grass.try_remove(tmp) 742 for tmp in six.itervalues(self.cats_rasts): 743 RunCommand("g.remove", flags='f', 744 type='raster', name=tmp, 745 getErrorMsg=True) 746 747 self.cats_rasts = {} 748 self.cats_rasts_conds = {} 749 750 def GetCatRast(self, cat_id): 751 if cat_id in self.cats_rasts: 752 return self.cats_rasts[cat_id] 753 return None 754 755 def GetCatRastCond(self, cat_id): 756 return self.cats_rasts_conds[cat_id] 757 758 def GetCatsRastsConds(self): 759 max_cat_id = max(self.cats_rasts_conds.keys()) 760 761 cats_rasts_conds = [''] * (max_cat_id + 1) 762 for i_cat_id, i_rast in six.iteritems(self.cats_rasts_conds): 763 cats_rasts_conds[i_cat_id] = i_rast 764 765 return cats_rasts_conds 766 767 def GetCatsRasts(self): 768 max_cat_id = max(self.cats_rasts.keys()) 769 770 cats_rasts = [''] * (max_cat_id + 1) 771 for i_cat_id, i_rast in six.iteritems(self.cats_rasts): 772 cats_rasts[i_cat_id] = i_rast 773 774 return cats_rasts 775 776 777# not used, using iclass_perimeter algorithm instead of scipy convolve2d 778""" 779def RasterizePolygon(pol, height, min_h, width, min_w): 780 781 # Joe Kington 782 # http://stackoverflow.com/questions/3654289/scipy-create-2d-polygon-mask 783 784 #poly_verts = [(1,1), (1,4), (4,4),(4,1), (1,1)] 785 786 nx = width 787 ny = height 788 789 x, y = np.meshgrid(np.arange(-0.5 + min_w, nx + 0.5 + min_w, dtype=float), 790 np.arange(-0.5 + min_h, ny + 0.5 + min_h, dtype=float)) 791 x, y = x.flatten(), y.flatten() 792 793 points = np.vstack((x,y)).T 794 795 p = Path(pol) 796 grid = p.contains_points(points) 797 grid = grid.reshape((ny + 1, nx + 1)) 798 raster = np.zeros((height, width), dtype=np.uint8)#TODO bool 799 800 #TODO shift by 0.5 801 B = np.ones((2,2))/4 802 raster = convolve2d(grid, B, 'valid') 803 804 return raster 805""" 806 807 808def idScattToidBands(scatt_id, n_bands): 809 """Get bands ids from scatter plot id.""" 810 n_b1 = n_bands - 1 811 812 band_1 = (int)( 813 (2 * n_b1 + 1 - sqrt(((2 * n_b1 + 1) * (2 * n_b1 + 1) - 8 * scatt_id))) / 2) 814 815 band_2 = int(scatt_id - (band_1 * (2 * n_b1 + 1) - band_1 * band_1) / 2 + band_1 + 1) 816 817 return band_1, band_2 818 819 820def idBandsToidScatt(band_1_id, band_2_id, n_bands): 821 """Get scatter plot id from band ids.""" 822 if band_2_id < band_1_id: 823 tmp = band_1_id 824 band_1_id = band_2_id 825 band_2_id = tmp 826 827 n_b1 = n_bands - 1 828 829 scatt_id = int(( 830 band_1_id * (2 * n_b1 + 1) - band_1_id * band_1_id) / 2 + band_2_id - band_1_id - 1) 831 832 return scatt_id 833 834 835def GetRegion(): 836 ret, region, msg = RunCommand("g.region", 837 flags="gp", 838 getErrorMsg=True, 839 read=True) 840 841 if ret != 0: 842 raise GException("g.region failed:\n%s" % msg) 843 844 return _parseRegion(region) 845 846 847def _parseRegion(region_str): 848 849 region = {} 850 region_str = region_str.splitlines() 851 852 for param in region_str: 853 k, v = param.split("=") 854 if k in ["rows", "cols", "cells"]: 855 v = int(v) 856 else: 857 v = float(v) 858 region[k] = v 859 860 return region 861 862 863def GetRasterInfo(rast): 864 ret, out, msg = RunCommand("r.info", 865 map=rast, 866 flags="rg", 867 getErrorMsg=True, 868 read=True) 869 870 if ret != 0: 871 raise GException("r.info failed:\n%s" % msg) 872 873 out = out.splitlines() 874 raster_info = {} 875 876 for b in out: 877 if not b.strip(): 878 continue 879 k, v = b.split("=") 880 if k == "datatype": 881 if v != "CELL": 882 return None 883 pass 884 elif k in ['rows', 'cols', 'cells', 'min', 'max']: 885 v = int(v) 886 else: 887 v = float(v) 888 889 raster_info[k] = v 890 891 raster_info['range'] = raster_info['max'] - raster_info['min'] + 1 892 return raster_info 893