iscatt_core.py 24 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854
  1. """
  2. @package iscatt.iscatt_core
  3. @brief Non GUI functions.
  4. Classes:
  5. - iscatt_core::Core
  6. - iscatt_core::CatRastUpdater
  7. - iscatt_core::AnalyzedData
  8. - iscatt_core::ScattPlotsCondsData
  9. - iscatt_core::ScattPlotsData
  10. (C) 2013 by the GRASS Development Team
  11. This program is free software under the GNU General Public License
  12. (>=v2). Read the file COPYING that comes with GRASS for details.
  13. @author Stepan Turek <stepan.turek seznam.cz> (mentor: Martin Landa)
  14. """
  15. import os
  16. import sys
  17. import time
  18. import numpy as np
  19. # used iclass perimeters algorithm instead of convolve2d
  20. #from matplotlib.path import Path
  21. #from scipy.signal import convolve2d
  22. from math import sqrt, ceil, floor
  23. from copy import deepcopy
  24. from core.gcmd import GException, GError, RunCommand
  25. import grass.script as grass
  26. from core_c import CreateCatRast, ComputeScatts, UpdateCatRast, \
  27. Rasterize, SC_SCATT_DATA, SC_SCATT_CONDITIONS
  28. MAX_SCATT_SIZE = 4100 * 4100
  29. WARN_SCATT_SIZE = 2000 * 2000
  30. MAX_NCELLS = 65536 * 65536
  31. WARN_NCELLS = 12000 * 12000
  32. class Core:
  33. """Represents scatter plot backend.
  34. """
  35. def __init__(self):
  36. self.an_data = AnalyzedData()
  37. self.scatts_dt = ScattPlotsData(self.an_data)
  38. self.scatt_conds_dt = ScattPlotsCondsData(self.an_data)
  39. self.cat_rast_updater = CatRastUpdater(self.scatts_dt, self.an_data, self)
  40. def SetData(self, bands):
  41. """Set bands for analysis.
  42. """
  43. ret = self.an_data.Create(bands)
  44. if not ret:
  45. return False
  46. n_bands = len(self.GetBands())
  47. self.scatts_dt.Create(n_bands)
  48. self.scatt_conds_dt.Create(n_bands)
  49. return True
  50. def AddCategory(self, cat_id):
  51. self.scatts_dt.AddCategory(cat_id)
  52. return self.scatt_conds_dt.AddCategory(cat_id)
  53. def DeleteCategory(self, cat_id):
  54. self.scatts_dt.DeleteCategory(cat_id)
  55. self.scatt_conds_dt.DeleteCategory(cat_id)
  56. def CleanUp(self):
  57. self.scatts_dt.CleanUp()
  58. self.scatt_conds_dt.CleanUp()
  59. def GetBands(self):
  60. return self.an_data.GetBands()
  61. def GetScattsData(self):
  62. return self.scatts_dt, self.scatt_conds_dt;
  63. def GetRegion(self):
  64. return self.an_data.GetRegion()
  65. def GetCatRast(self, cat_id):
  66. return self.scatts_dt.GetCatRast(cat_id)
  67. def AddScattPlots(self, scatt_ids):
  68. for s_id in scatt_ids:
  69. self.scatts_dt.AddScattPlot(scatt_id = s_id)
  70. cats_ids = self.scatts_dt.GetCategories()
  71. self.ComputeCatsScatts(cats_ids)
  72. def SetEditCatData(self, cat_id, scatt_id, bbox, value):
  73. if cat_id not in self.scatts_dt.GetCategories():
  74. raise GException(_("Select category for editing."))
  75. if self.scatt_conds_dt.AddScattPlot(cat_id, scatt_id) < 0:
  76. return None
  77. arr = self.scatt_conds_dt.GetValuesArr(cat_id, scatt_id)
  78. for k, v in bbox.iteritems():
  79. bbox[k] = self._validExtend(v)
  80. arr[bbox['btm_y'] : bbox['up_y'], bbox['btm_x'] : bbox['up_x']] = value
  81. self.ComputeCatsScatts([cat_id])
  82. return cat_id
  83. def ComputeCatsScatts(self, cats_ids):
  84. requested_dt = {}
  85. requested_dt_conds = {}
  86. for c in cats_ids:
  87. requested_dt_conds[c] = self.scatt_conds_dt.GetCatScatts(c)
  88. requested_dt[c] = self.scatts_dt.GetCatScatts(c)
  89. scatt_conds = self.scatt_conds_dt.GetData(requested_dt_conds)
  90. scatts = self.scatts_dt.GetData(requested_dt)
  91. bands = self.an_data.GetBands()
  92. cats_rasts = self.scatts_dt.GetCatsRasts()
  93. cats_rasts_conds = self.scatts_dt.GetCatsRastsConds()
  94. returncode, scatts = ComputeScatts(self.an_data.GetRegion(),
  95. scatt_conds,
  96. bands,
  97. len(self.GetBands()),
  98. scatts,
  99. cats_rasts_conds,
  100. cats_rasts)
  101. if returncode < 0:
  102. GException(_("Computing of scatter plots failed."))
  103. def CatRastUpdater(self):
  104. return self.cat_rast_updater
  105. def UpdateCategoryWithPolygons(self, cat_id, scatts_pols, value):
  106. start_time = time.clock()
  107. if cat_id not in self.scatts_dt.GetCategories():
  108. raise GException(_("Select category for editing."))
  109. for scatt_id, coords in scatts_pols.iteritems():
  110. if self.scatt_conds_dt.AddScattPlot(cat_id, scatt_id) < 0:
  111. return False
  112. b1, b2 = idScattToidBands(scatt_id, len(self.an_data.GetBands()))
  113. b = self.scatts_dt.GetBandsInfo(scatt_id)
  114. region = {}
  115. region['s'] = b['b2']['min'] - 0.5
  116. region['n'] = b['b2']['max'] + 0.5
  117. region['w'] = b['b1']['min'] - 0.5
  118. region['e'] = b['b1']['max'] + 0.5
  119. arr = self.scatt_conds_dt.GetValuesArr(cat_id, scatt_id)
  120. arr = Rasterize(polygon=coords,
  121. rast=arr,
  122. region=region,
  123. value=value)
  124. # previous way of rasterization / used scipy
  125. #raster_pol = RasterizePolygon(coords, b['b1']['range'], b['b1']['min'],
  126. # b['b2']['range'], b['b2']['min'])
  127. #raster_ind = np.where(raster_pol > 0)
  128. #arr = self.scatt_conds_dt.GetValuesArr(cat_id, scatt_id)
  129. #arr[raster_ind] = value
  130. #arr.flush()
  131. self.ComputeCatsScatts([cat_id])
  132. return cat_id
  133. def ExportCatRast(self, cat_id, rast_name):
  134. cat_rast = self.scatts_dt.GetCatRast(cat_id);
  135. if not cat_rast:
  136. return 1
  137. return RunCommand("g.copy",
  138. raster=cat_rast + ',' + rast_name,
  139. getErrorMsg=True,
  140. overwrite=True)
  141. def _validExtend(self, val):
  142. #TODO do it general
  143. if val > 255:
  144. val = 255
  145. elif val < 0:
  146. val = 0
  147. return val
  148. class CatRastUpdater:
  149. """Update backend data structures according to selected areas in mapwindow.
  150. """
  151. def __init__(self, scatts_dt, an_data, core):
  152. self.scatts_dt = scatts_dt
  153. self.an_data = an_data # TODO may be confusing
  154. self.core = core
  155. self.vectMap = None
  156. def SetVectMap(self, vectMap):
  157. self.vectMap = vectMap
  158. def SyncWithMap(self):
  159. #TODO possible optimization - bbox only of vertex and its two neighbours
  160. region = self.an_data.GetRegion()
  161. bbox = {}
  162. bbox['maxx'] = region['e']
  163. bbox['minx'] = region['w']
  164. bbox['maxy'] = region['n']
  165. bbox['miny'] = region['s']
  166. updated_cats = []
  167. for cat_id in self.scatts_dt.GetCategories():
  168. if cat_id == 0:
  169. continue
  170. cat = [{1 : [cat_id]}]
  171. self._updateCatRast(bbox, cat, updated_cats)
  172. return updated_cats
  173. def EditedFeature(self, new_bboxs, new_areas_cats, old_bboxs, old_areas_cats):
  174. #TODO possible optimization - bbox only of vertex and its two neighbours
  175. bboxs = old_bboxs + new_bboxs
  176. areas_cats = old_areas_cats + new_areas_cats
  177. updated_cats = []
  178. for i in range(len(areas_cats)):
  179. self._updateCatRast(bboxs[i], areas_cats[i], updated_cats)
  180. return updated_cats
  181. def _updateCatRast(self, bbox, areas_cats, updated_cats):
  182. rasterized_cats = []
  183. for c in range(len(areas_cats)):
  184. if not areas_cats[c]:
  185. continue
  186. layer = areas_cats[c].keys()[0]
  187. cat = areas_cats[c][layer][0]
  188. if cat in rasterized_cats:
  189. continue
  190. rasterized_cats.append(cat)
  191. updated_cats.append(cat)
  192. grass_region = self._create_grass_region_env(bbox)
  193. #TODO hack check if raster exists?
  194. patch_rast = "temp_scatt_patch_%d" % (os.getpid())
  195. self._rasterize(grass_region, layer, cat, patch_rast)
  196. region = self.an_data.GetRegion()
  197. ret = UpdateCatRast(patch_rast, region, self.scatts_dt.GetCatRastCond(cat))
  198. if ret < 0:
  199. GException(_("Patching category raster conditions file failed."))
  200. RunCommand("g.remove", flags='f', type='raster', name=patch_rast)
  201. def _rasterize(self, grass_region, layer, cat, out_rast):
  202. #TODO different thread may be problem when user edits map
  203. environs = os.environ.copy()
  204. environs['GRASS_VECTOR_TEMPORARY'] = '1'
  205. ret, text, msg = RunCommand("v.category",
  206. input=self.vectMap,
  207. getErrorMsg = True,
  208. option='report',
  209. read=True,
  210. env=environs)
  211. ret, text, msg = RunCommand("v.build",
  212. map = self.vectMap,
  213. getErrorMsg = True,
  214. read = True,
  215. env = environs)
  216. if ret != 0:
  217. GException(_("v.build failed:\n%s" % msg))
  218. environs = os.environ.copy()
  219. environs["GRASS_REGION"] = grass_region["GRASS_REGION"]
  220. environs['GRASS_VECTOR_TEMPORARY'] = '1'
  221. ret, text, msg = RunCommand("v.to.rast",
  222. input = self.vectMap,
  223. use = "cat",
  224. layer = str(layer),
  225. cat = str(cat),
  226. output = out_rast,
  227. getErrorMsg = True,
  228. read = True,
  229. overwrite = True,
  230. env = environs)
  231. if ret != 0:
  232. GException(_("v.to.rast failed:\n%s" % msg))
  233. def _create_grass_region_env(self, bbox):
  234. r = self.an_data.GetRegion()
  235. new_r = {}
  236. if bbox["maxy"] <= r["s"]:
  237. return 0
  238. elif bbox["maxy"] >= r["n"]:
  239. new_r["n"] = bbox["maxy"]
  240. else:
  241. new_r["n"] = ceil((bbox["maxy"] - r["s"]) / r["nsres"]) * r["nsres"] + r["s"]
  242. if bbox["miny"] >= r["n"]:
  243. return 0
  244. elif bbox["miny"] <= r["s"]:
  245. new_r["s"] = bbox["miny"]
  246. else:
  247. new_r["s"] = floor((bbox["miny"] - r["s"]) / r["nsres"]) * r["nsres"] + r["s"]
  248. if bbox["maxx"] <= r["w"]:
  249. return 0
  250. elif bbox["maxx"] >= r["e"]:
  251. new_r["e"] = bbox["maxx"]
  252. else:
  253. new_r["e"] = ceil((bbox["maxx"] - r["w"]) / r["ewres"]) * r["ewres"] + r["w"]
  254. if bbox["minx"] >= r["e"]:
  255. return 0
  256. elif bbox["minx"] <= r["w"]:
  257. new_r["w"] = bbox["minx"]
  258. else:
  259. new_r["w"] = floor((bbox["minx"] - r["w"]) / r["ewres"]) * r["ewres"] + r["w"]
  260. #TODO check regions resolution
  261. new_r["nsres"] = r["nsres"]
  262. new_r["ewres"] = r["ewres"]
  263. return {"GRASS_REGION" : grass.region_env(**new_r)}
  264. class AnalyzedData:
  265. """Represents analyzed data (bands, region).
  266. """
  267. def __init__(self):
  268. self.bands = []
  269. self.bands_info = {}
  270. self.region = None
  271. def GetRegion(self):
  272. return self.region
  273. def Create(self, bands):
  274. self.bands = bands[:]
  275. self.region = None
  276. self.region = GetRegion()
  277. if self.region["rows"] * self.region["cols"] > MAX_NCELLS:
  278. GException("too big region")
  279. self.bands_info = {}
  280. for b in self.bands[:]:
  281. i = GetRasterInfo(b)
  282. if i is None:
  283. GException("raster %s is not CELL type" % (b))
  284. self.bands_info[b] = i
  285. #TODO check size of raster
  286. return True
  287. def GetBands(self):
  288. return self.bands
  289. def GetBandInfo(self, band_id):
  290. band = self.bands[band_id]
  291. return self.bands_info[band]
  292. class ScattPlotsCondsData:
  293. """Data structure for selected areas in scatter plot(condtions).
  294. """
  295. def __init__(self, an_data):
  296. self.an_data = an_data
  297. #TODO
  298. self.max_n_cats = 10
  299. self.dtype = 'uint8'
  300. self.type = 1;
  301. self.CleanUp()
  302. def CleanUp(self):
  303. self.cats = {}
  304. self.n_scatts = -1
  305. self.n_bands = -1
  306. for cat_id in self.cats.keys():
  307. self.DeleteCategory(cat_id)
  308. def Create(self, n_bands):
  309. self.CleanUp()
  310. self.n_scatts = (n_bands - 1) * n_bands / 2;
  311. self.n_bands = n_bands
  312. self.AddCategory(cat_id = 0)
  313. def AddCategory(self, cat_id):
  314. if cat_id not in self.cats.keys():
  315. self.cats[cat_id] = {}
  316. return cat_id
  317. return -1
  318. def DeleteCategory(self, cat_id):
  319. if cat_id not in self.cats.keys():
  320. return False
  321. for scatt in self.cats[cat_id].itervalues():
  322. grass.try_remove(scatt['np_vals'])
  323. del scatt['np_vals']
  324. del self.cats[cat_id]
  325. return True
  326. def GetCategories(self):
  327. return self.cats.keys()
  328. def GetCatScatts(self, cat_id):
  329. if not self.cats.has_key(cat_id):
  330. return False
  331. return self.cats[cat_id].keys()
  332. def AddScattPlot(self, cat_id, scatt_id):
  333. if not self.cats.has_key(cat_id):
  334. return -1
  335. if self.cats[cat_id].has_key(scatt_id):
  336. return 0
  337. b_i = self.GetBandsInfo(scatt_id)
  338. shape = (b_i['b2']['max'] - b_i['b2']['min'] + 1, b_i['b1']['max'] - b_i['b1']['min'] + 1)
  339. np_vals = np.memmap(grass.tempfile(), dtype=self.dtype, mode='w+', shape=shape)
  340. self.cats[cat_id][scatt_id] = {'np_vals' : np_vals}
  341. return 1
  342. def GetBandsInfo(self, scatt_id):
  343. b1, b2 = idScattToidBands(scatt_id, len(self.an_data.GetBands()))
  344. b1_info = self.an_data.GetBandInfo(b1)
  345. b2_info = self.an_data.GetBandInfo(b2)
  346. bands_info = {'b1' : b1_info,
  347. 'b2' : b2_info}
  348. return bands_info
  349. def DeleScattPlot(self, cat_id, scatt_id):
  350. if not self.cats.has_key(cat_id):
  351. return False
  352. if not self.cats[cat_id].has_key(scatt_id):
  353. return False
  354. del self.cats[cat_id][scatt_id]
  355. return True
  356. def GetValuesArr(self, cat_id, scatt_id):
  357. if not self.cats.has_key(cat_id):
  358. return None
  359. if not self.cats[cat_id].has_key(scatt_id):
  360. return None
  361. return self.cats[cat_id][scatt_id]['np_vals']
  362. def GetData(self, requested_dt):
  363. cats = {}
  364. for cat_id, scatt_ids in requested_dt.iteritems():
  365. if not cats.has_key(cat_id):
  366. cats[cat_id] = {}
  367. for scatt_id in scatt_ids:
  368. # if key is missing condition is always True (full scatter plor is computed)
  369. if self.cats[cat_id].has_key(scatt_id):
  370. cats[cat_id][scatt_id] = {'np_vals' : self.cats[cat_id][scatt_id]['np_vals'],
  371. 'bands_info' : self.GetBandsInfo(scatt_id)}
  372. return cats
  373. def SetData(self, cats):
  374. for cat_id, scatt_ids in cats.iteritems():
  375. for scatt_id in scatt_ids:
  376. # if key is missing condition is always True (full scatter plor is computed)
  377. if self.cats[cat_id].has_key(scatt_id):
  378. self.cats[cat_id][scatt_id]['np_vals'] = cats[cat_id][scatt_id]['np_vals']
  379. def GetScatt(self, scatt_id, cats_ids = None):
  380. scatts = {}
  381. for cat_id in self.cats.iterkeys():
  382. if cats_ids and cat_id not in cats_ids:
  383. continue
  384. if not self.cats[cat_id].has_key(scatt_id):
  385. continue
  386. scatts[cat_id] = {'np_vals' : self.cats[cat_id][scatt_id]['np_vals'],
  387. 'bands_info' : self.GetBandsInfo(scatt_id)}
  388. return scatts
  389. class ScattPlotsData(ScattPlotsCondsData):
  390. """Data structure for computed points (classes) in scatter plots.\
  391. """
  392. def __init__(self, an_data):
  393. self.cats_rasts = {}
  394. self.cats_rasts_conds = {}
  395. self.scatts_ids = []
  396. ScattPlotsCondsData.__init__(self, an_data)
  397. self.dtype = 'uint32'
  398. #TODO
  399. self.type = 0
  400. def AddCategory(self, cat_id):
  401. cat_id = ScattPlotsCondsData.AddCategory(self, cat_id)
  402. if cat_id < 0:
  403. return cat_id
  404. for scatt_id in self.scatts_ids:
  405. ScattPlotsCondsData.AddScattPlot(self, cat_id, scatt_id)
  406. if cat_id == 0:
  407. self.cats_rasts_conds[cat_id] = None
  408. self.cats_rasts[cat_id] = None
  409. else:
  410. self.cats_rasts_conds[cat_id] = grass.tempfile()
  411. self.cats_rasts[cat_id] = "temp_cat_rast_%d_%d" % (cat_id, os.getpid())
  412. region = self.an_data.GetRegion()
  413. CreateCatRast(region, self.cats_rasts_conds[cat_id])
  414. return cat_id
  415. def DeleteCategory(self, cat_id):
  416. ScattPlotsCondsData.DeleteCategory(self, cat_id)
  417. grass.try_remove(self.cats_rasts_conds[cat_id])
  418. del self.cats_rasts_conds[cat_id]
  419. RunCommand("g.remove", flags='f', type='raster',
  420. name=self.cats_rasts[cat_id])
  421. del self.cats_rasts[cat_id]
  422. return True
  423. def AddScattPlot(self, scatt_id):
  424. if scatt_id in self.scatts_ids:
  425. return False
  426. self.scatts_ids.append(scatt_id)
  427. for cat_id in self.cats.iterkeys():
  428. ScattPlotsCondsData.AddScattPlot(self, cat_id, scatt_id)
  429. self.cats[cat_id][scatt_id]['ellipse'] = None
  430. return True
  431. def DeleteScatterPlot(self, scatt_id):
  432. if scatt_id not in self.scatts_ids:
  433. return False
  434. self.scatts_ids.remove(scatt_id)
  435. for cat_id in self.cats.iterkeys():
  436. ScattPlotsCondsData.DeleteScattPlot(self, cat_id, scatt_id)
  437. return True
  438. def GetEllipses(self, scatt_id, styles):
  439. if scatt_id not in self.scatts_ids:
  440. return False
  441. scatts = {}
  442. for cat_id in self.cats.iterkeys():
  443. if cat_id == 0:
  444. continue
  445. nstd = styles[cat_id]['nstd']
  446. scatts[cat_id] = self._getEllipse(cat_id, scatt_id, nstd)
  447. return scatts
  448. def _getEllipse(self, cat_id, scatt_id, nstd):
  449. # Joe Kington
  450. # http://stackoverflow.com/questions/12301071/multidimensional-confidence-intervals
  451. data = np.copy(self.cats[cat_id][scatt_id]['np_vals'])
  452. b = self.GetBandsInfo(scatt_id)
  453. sel_pts = np.where(data > 0)
  454. x = sel_pts[1]
  455. y = sel_pts[0]
  456. flatten_data = data.reshape([-1])
  457. flatten_sel_pts = np.nonzero(flatten_data)
  458. weights = flatten_data[flatten_sel_pts]
  459. if len(weights) == 0:
  460. return None
  461. x_avg = np.average(x, weights=weights)
  462. y_avg = np.average(y, weights=weights)
  463. pos = np.array([x_avg + b['b1']['min'], y_avg + b['b2']['min']])
  464. x_diff = (x - x_avg)
  465. y_diff = (y - y_avg)
  466. x_diff = (x - x_avg)
  467. y_diff = (y - y_avg)
  468. diffs = x_diff * y_diff.T
  469. cov = np.dot(diffs, weights) / (np.sum(weights) - 1)
  470. diffs = x_diff * x_diff.T
  471. var_x = np.dot(diffs, weights) / (np.sum(weights) - 1)
  472. diffs = y_diff * y_diff.T
  473. var_y = np.dot(diffs, weights) / (np.sum(weights) - 1)
  474. cov = np.array([[var_x, cov],[cov, var_y]])
  475. def eigsorted(cov):
  476. vals, vecs = np.linalg.eigh(cov)
  477. order = vals.argsort()[::-1]
  478. return vals[order], vecs[:,order]
  479. vals, vecs = eigsorted(cov)
  480. theta = np.degrees(np.arctan2(*vecs[:,0][::-1]))
  481. # Width and height are "full" widths, not radius
  482. width, height = 2 * nstd * np.sqrt(vals)
  483. ellipse = {'pos' : pos,
  484. 'width' : width,
  485. 'height' : height,
  486. 'theta' : theta}
  487. del data
  488. del flatten_data
  489. del flatten_sel_pts
  490. del weights
  491. del sel_pts
  492. return ellipse
  493. def CleanUp(self):
  494. ScattPlotsCondsData.CleanUp(self)
  495. for tmp in self.cats_rasts_conds.itervalues():
  496. grass.try_remove(tmp)
  497. for tmp in self.cats_rasts.itervalues():
  498. RunCommand("g.remove", flags='f',
  499. type='raster', name=tmp,
  500. getErrorMsg=True)
  501. self.cats_rasts = {}
  502. self.cats_rasts_conds = {}
  503. def GetCatRast(self, cat_id):
  504. if self.cats_rasts.has_key(cat_id):
  505. return self.cats_rasts[cat_id]
  506. return None
  507. def GetCatRastCond(self, cat_id):
  508. return self.cats_rasts_conds[cat_id]
  509. def GetCatsRastsConds(self):
  510. max_cat_id = max(self.cats_rasts_conds.keys())
  511. cats_rasts_conds = [''] * (max_cat_id + 1)
  512. for i_cat_id, i_rast in self.cats_rasts_conds.iteritems():
  513. cats_rasts_conds[i_cat_id] = i_rast
  514. return cats_rasts_conds
  515. def GetCatsRasts(self):
  516. max_cat_id = max(self.cats_rasts.keys())
  517. cats_rasts = [''] * (max_cat_id + 1)
  518. for i_cat_id, i_rast in self.cats_rasts.iteritems():
  519. cats_rasts[i_cat_id] = i_rast
  520. return cats_rasts
  521. # not used, using iclass_perimeter algorithm instead of scipy convolve2d
  522. """
  523. def RasterizePolygon(pol, height, min_h, width, min_w):
  524. # Joe Kington
  525. # http://stackoverflow.com/questions/3654289/scipy-create-2d-polygon-mask
  526. #poly_verts = [(1,1), (1,4), (4,4),(4,1), (1,1)]
  527. nx = width
  528. ny = height
  529. x, y = np.meshgrid(np.arange(-0.5 + min_w, nx + 0.5 + min_w, dtype=float),
  530. np.arange(-0.5 + min_h, ny + 0.5 + min_h, dtype=float))
  531. x, y = x.flatten(), y.flatten()
  532. points = np.vstack((x,y)).T
  533. p = Path(pol)
  534. grid = p.contains_points(points)
  535. grid = grid.reshape((ny + 1, nx + 1))
  536. raster = np.zeros((height, width), dtype=np.uint8)#TODO bool
  537. #TODO shift by 0.5
  538. B = np.ones((2,2))/4
  539. raster = convolve2d(grid, B, 'valid')
  540. return raster
  541. """
  542. def idScattToidBands(scatt_id, n_bands):
  543. """Get bands ids from scatter plot id."""
  544. n_b1 = n_bands - 1
  545. band_1 = (int) ((2 * n_b1 + 1 - sqrt(((2 * n_b1 + 1) * (2 * n_b1 + 1) - 8 * scatt_id))) / 2)
  546. band_2 = scatt_id - (band_1 * (2 * n_b1 + 1) - band_1 * band_1) / 2 + band_1 + 1
  547. return band_1, band_2
  548. def idBandsToidScatt(band_1_id, band_2_id, n_bands):
  549. """Get scatter plot id from band ids."""
  550. if band_2_id < band_1_id:
  551. tmp = band_1_id
  552. band_1_id = band_2_id
  553. band_2_id = tmp
  554. n_b1 = n_bands - 1
  555. scatt_id = (band_1_id * (2 * n_b1 + 1) - band_1_id * band_1_id) / 2 + band_2_id - band_1_id - 1
  556. return scatt_id
  557. def GetRegion():
  558. ret, region, msg = RunCommand("g.region",
  559. flags = "gp",
  560. getErrorMsg = True,
  561. read = True)
  562. if ret != 0:
  563. raise GException("g.region failed:\n%s" % msg)
  564. return _parseRegion(region)
  565. def _parseRegion(region_str):
  566. region = {}
  567. region_str = region_str.splitlines()
  568. for param in region_str:
  569. k, v = param.split("=")
  570. if k in ["rows", "cols", "cells"]:
  571. v = int(v)
  572. else:
  573. v = float(v)
  574. region[k] = v
  575. return region
  576. def GetRasterInfo(rast):
  577. ret, out, msg = RunCommand("r.info",
  578. map = rast,
  579. flags = "rg",
  580. getErrorMsg = True,
  581. read = True)
  582. if ret != 0:
  583. raise GException("r.info failed:\n%s" % msg)
  584. out = out.split("\n")
  585. raster_info = {}
  586. for b in out:
  587. if not b.strip():
  588. continue
  589. k, v = b.split("=")
  590. if k == "datatype":
  591. if v != "CELL":
  592. return None
  593. pass
  594. elif k in ['rows', 'cols', 'cells', 'min', 'max']:
  595. v = int(v)
  596. else:
  597. v = float(v)
  598. raster_info[k] = v
  599. raster_info['range'] = raster_info['max'] - raster_info['min'] + 1
  600. return raster_info