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. rast=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",
  201. rast = patch_rast)
  202. def _rasterize(self, grass_region, layer, cat, out_rast):
  203. #TODO different thread may be problem when user edits map
  204. environs = os.environ.copy()
  205. environs['GRASS_VECTOR_TEMPORARY'] = '1'
  206. ret, text, msg = RunCommand("v.category",
  207. input=self.vectMap,
  208. getErrorMsg = True,
  209. option='report',
  210. read=True,
  211. env=environs)
  212. ret, text, msg = RunCommand("v.build",
  213. map = self.vectMap,
  214. getErrorMsg = True,
  215. read = True,
  216. env = environs)
  217. if ret != 0:
  218. GException(_("v.build failed:\n%s" % msg))
  219. environs = os.environ.copy()
  220. environs["GRASS_REGION"] = grass_region["GRASS_REGION"]
  221. environs['GRASS_VECTOR_TEMPORARY'] = '1'
  222. ret, text, msg = RunCommand("v.to.rast",
  223. input = self.vectMap,
  224. use = "cat",
  225. layer = str(layer),
  226. cat = str(cat),
  227. output = out_rast,
  228. getErrorMsg = True,
  229. read = True,
  230. overwrite = True,
  231. env = environs)
  232. if ret != 0:
  233. GException(_("v.to.rast failed:\n%s" % msg))
  234. def _create_grass_region_env(self, bbox):
  235. r = self.an_data.GetRegion()
  236. new_r = {}
  237. if bbox["maxy"] <= r["s"]:
  238. return 0
  239. elif bbox["maxy"] >= r["n"]:
  240. new_r["n"] = bbox["maxy"]
  241. else:
  242. new_r["n"] = ceil((bbox["maxy"] - r["s"]) / r["nsres"]) * r["nsres"] + r["s"]
  243. if bbox["miny"] >= r["n"]:
  244. return 0
  245. elif bbox["miny"] <= r["s"]:
  246. new_r["s"] = bbox["miny"]
  247. else:
  248. new_r["s"] = floor((bbox["miny"] - r["s"]) / r["nsres"]) * r["nsres"] + r["s"]
  249. if bbox["maxx"] <= r["w"]:
  250. return 0
  251. elif bbox["maxx"] >= r["e"]:
  252. new_r["e"] = bbox["maxx"]
  253. else:
  254. new_r["e"] = ceil((bbox["maxx"] - r["w"]) / r["ewres"]) * r["ewres"] + r["w"]
  255. if bbox["minx"] >= r["e"]:
  256. return 0
  257. elif bbox["minx"] <= r["w"]:
  258. new_r["w"] = bbox["minx"]
  259. else:
  260. new_r["w"] = floor((bbox["minx"] - r["w"]) / r["ewres"]) * r["ewres"] + r["w"]
  261. #TODO check regions resolution
  262. new_r["nsres"] = r["nsres"]
  263. new_r["ewres"] = r["ewres"]
  264. return {"GRASS_REGION" : grass.region_env(**new_r)}
  265. class AnalyzedData:
  266. """!Represents analyzed data (bands, region).
  267. """
  268. def __init__(self):
  269. self.bands = []
  270. self.bands_info = {}
  271. self.region = None
  272. def GetRegion(self):
  273. return self.region
  274. def Create(self, bands):
  275. self.bands = bands[:]
  276. self.region = None
  277. self.region = GetRegion()
  278. if self.region["rows"] * self.region["cols"] > MAX_NCELLS:
  279. GException("too big region")
  280. self.bands_info = {}
  281. for b in self.bands[:]:
  282. i = GetRasterInfo(b)
  283. if i is None:
  284. GException("raster %s is not CELL type" % (b))
  285. self.bands_info[b] = i
  286. #TODO check size of raster
  287. return True
  288. def GetBands(self):
  289. return self.bands
  290. def GetBandInfo(self, band_id):
  291. band = self.bands[band_id]
  292. return self.bands_info[band]
  293. class ScattPlotsCondsData:
  294. """!Data structure for selected areas in scatter plot(condtions).
  295. """
  296. def __init__(self, an_data):
  297. self.an_data = an_data
  298. #TODO
  299. self.max_n_cats = 10
  300. self.dtype = 'uint8'
  301. self.type = 1;
  302. self.CleanUp()
  303. def CleanUp(self):
  304. self.cats = {}
  305. self.n_scatts = -1
  306. self.n_bands = -1
  307. for cat_id in self.cats.keys():
  308. self.DeleteCategory(cat_id)
  309. def Create(self, n_bands):
  310. self.CleanUp()
  311. self.n_scatts = (n_bands - 1) * n_bands / 2;
  312. self.n_bands = n_bands
  313. self.AddCategory(cat_id = 0)
  314. def AddCategory(self, cat_id):
  315. if cat_id not in self.cats.keys():
  316. self.cats[cat_id] = {}
  317. return cat_id
  318. return -1
  319. def DeleteCategory(self, cat_id):
  320. if cat_id not in self.cats.keys():
  321. return False
  322. for scatt in self.cats[cat_id].itervalues():
  323. grass.try_remove(scatt['np_vals'])
  324. del scatt['np_vals']
  325. del self.cats[cat_id]
  326. return True
  327. def GetCategories(self):
  328. return self.cats.keys()
  329. def GetCatScatts(self, cat_id):
  330. if not self.cats.has_key(cat_id):
  331. return False
  332. return self.cats[cat_id].keys()
  333. def AddScattPlot(self, cat_id, scatt_id):
  334. if not self.cats.has_key(cat_id):
  335. return -1
  336. if self.cats[cat_id].has_key(scatt_id):
  337. return 0
  338. b_i = self.GetBandsInfo(scatt_id)
  339. shape = (b_i['b2']['max'] - b_i['b2']['min'] + 1, b_i['b1']['max'] - b_i['b1']['min'] + 1)
  340. np_vals = np.memmap(grass.tempfile(), dtype=self.dtype, mode='w+', shape = shape)
  341. self.cats[cat_id][scatt_id] = {'np_vals' : np_vals}
  342. return 1
  343. def GetBandsInfo(self, scatt_id):
  344. b1, b2 = idScattToidBands(scatt_id, len(self.an_data.GetBands()))
  345. b1_info = self.an_data.GetBandInfo(b1)
  346. b2_info = self.an_data.GetBandInfo(b2)
  347. bands_info = {'b1' : b1_info,
  348. 'b2' : b2_info}
  349. return bands_info
  350. def DeleScattPlot(self, cat_id, scatt_id):
  351. if not self.cats.has_key(cat_id):
  352. return False
  353. if not self.cats[cat_id].has_key(scatt_id):
  354. return False
  355. del self.cats[cat_id][scatt_id]
  356. return True
  357. def GetValuesArr(self, cat_id, scatt_id):
  358. if not self.cats.has_key(cat_id):
  359. return None
  360. if not self.cats[cat_id].has_key(scatt_id):
  361. return None
  362. return self.cats[cat_id][scatt_id]['np_vals']
  363. def GetData(self, requested_dt):
  364. cats = {}
  365. for cat_id, scatt_ids in requested_dt.iteritems():
  366. if not cats.has_key(cat_id):
  367. cats[cat_id] = {}
  368. for scatt_id in scatt_ids:
  369. # if key is missing condition is always True (full scatter plor is computed)
  370. if self.cats[cat_id].has_key(scatt_id):
  371. cats[cat_id][scatt_id] = {'np_vals' : self.cats[cat_id][scatt_id]['np_vals'],
  372. 'bands_info' : self.GetBandsInfo(scatt_id)}
  373. return cats
  374. def SetData(self, cats):
  375. for cat_id, scatt_ids in cats.iteritems():
  376. for scatt_id in scatt_ids:
  377. # if key is missing condition is always True (full scatter plor is computed)
  378. if self.cats[cat_id].has_key(scatt_id):
  379. self.cats[cat_id][scatt_id]['np_vals'] = cats[cat_id][scatt_id]['np_vals']
  380. def GetScatt(self, scatt_id, cats_ids = None):
  381. scatts = {}
  382. for cat_id in self.cats.iterkeys():
  383. if cats_ids and cat_id not in cats_ids:
  384. continue
  385. if not self.cats[cat_id].has_key(scatt_id):
  386. continue
  387. scatts[cat_id] = {'np_vals' : self.cats[cat_id][scatt_id]['np_vals'],
  388. 'bands_info' : self.GetBandsInfo(scatt_id)}
  389. return scatts
  390. class ScattPlotsData(ScattPlotsCondsData):
  391. """!Data structure for computed points (classes) in scatter plots.\
  392. """
  393. def __init__(self, an_data):
  394. self.cats_rasts = {}
  395. self.cats_rasts_conds = {}
  396. self.scatts_ids = []
  397. ScattPlotsCondsData.__init__(self, an_data)
  398. self.dtype = 'uint32'
  399. #TODO
  400. self.type = 0
  401. def AddCategory(self, cat_id):
  402. cat_id = ScattPlotsCondsData.AddCategory(self, cat_id)
  403. if cat_id < 0:
  404. return cat_id
  405. for scatt_id in self.scatts_ids:
  406. ScattPlotsCondsData.AddScattPlot(self, cat_id, scatt_id)
  407. if cat_id == 0:
  408. self.cats_rasts_conds[cat_id] = None
  409. self.cats_rasts[cat_id] = None
  410. else:
  411. self.cats_rasts_conds[cat_id] = grass.tempfile()
  412. self.cats_rasts[cat_id] = "temp_cat_rast_%d_%d" % (cat_id, os.getpid())
  413. region = self.an_data.GetRegion()
  414. CreateCatRast(region, self.cats_rasts_conds[cat_id])
  415. return cat_id
  416. def DeleteCategory(self, cat_id):
  417. ScattPlotsCondsData.DeleteCategory(self, cat_id)
  418. grass.try_remove(self.cats_rasts_conds[cat_id])
  419. del self.cats_rasts_conds[cat_id]
  420. RunCommand("g.remove",
  421. rast=self.cats_rasts[cat_id])
  422. del self.cats_rasts[cat_id]
  423. return True
  424. def AddScattPlot(self, scatt_id):
  425. if scatt_id in self.scatts_ids:
  426. return False
  427. self.scatts_ids.append(scatt_id)
  428. for cat_id in self.cats.iterkeys():
  429. ScattPlotsCondsData.AddScattPlot(self, cat_id, scatt_id)
  430. self.cats[cat_id][scatt_id]['ellipse'] = None
  431. return True
  432. def DeleteScatterPlot(self, scatt_id):
  433. if scatt_id not in self.scatts_ids:
  434. return False
  435. self.scatts_ids.remove(scatt_id)
  436. for cat_id in self.cats.iterkeys():
  437. ScattPlotsCondsData.DeleteScattPlot(self, cat_id, scatt_id)
  438. return True
  439. def GetEllipses(self, scatt_id, styles):
  440. if scatt_id not in self.scatts_ids:
  441. return False
  442. scatts = {}
  443. for cat_id in self.cats.iterkeys():
  444. if cat_id == 0:
  445. continue
  446. nstd = styles[cat_id]['nstd']
  447. scatts[cat_id] = self._getEllipse(cat_id, scatt_id, nstd)
  448. return scatts
  449. def _getEllipse(self, cat_id, scatt_id, nstd):
  450. # Joe Kington
  451. # http://stackoverflow.com/questions/12301071/multidimensional-confidence-intervals
  452. data = np.copy(self.cats[cat_id][scatt_id]['np_vals'])
  453. b = self.GetBandsInfo(scatt_id)
  454. sel_pts = np.where(data > 0)
  455. x = sel_pts[1]
  456. y = sel_pts[0]
  457. flatten_data = data.reshape([-1])
  458. flatten_sel_pts = np.nonzero(flatten_data)
  459. weights = flatten_data[flatten_sel_pts]
  460. if len(weights) == 0:
  461. return None
  462. x_avg = np.average(x, weights=weights)
  463. y_avg = np.average(y, weights=weights)
  464. pos = np.array([x_avg + b['b1']['min'], y_avg + b['b2']['min']])
  465. x_diff = (x - x_avg)
  466. y_diff = (y - y_avg)
  467. x_diff = (x - x_avg)
  468. y_diff = (y - y_avg)
  469. diffs = x_diff * y_diff.T
  470. cov = np.dot(diffs, weights) / (np.sum(weights) - 1)
  471. diffs = x_diff * x_diff.T
  472. var_x = np.dot(diffs, weights) / (np.sum(weights) - 1)
  473. diffs = y_diff * y_diff.T
  474. var_y = np.dot(diffs, weights) / (np.sum(weights) - 1)
  475. cov = np.array([[var_x, cov],[cov, var_y]])
  476. def eigsorted(cov):
  477. vals, vecs = np.linalg.eigh(cov)
  478. order = vals.argsort()[::-1]
  479. return vals[order], vecs[:,order]
  480. vals, vecs = eigsorted(cov)
  481. theta = np.degrees(np.arctan2(*vecs[:,0][::-1]))
  482. # Width and height are "full" widths, not radius
  483. width, height = 2 * nstd * np.sqrt(vals)
  484. ellipse = {'pos' : pos,
  485. 'width' : width,
  486. 'height' : height,
  487. 'theta' : theta}
  488. del data
  489. del flatten_data
  490. del flatten_sel_pts
  491. del weights
  492. del sel_pts
  493. return ellipse
  494. def CleanUp(self):
  495. ScattPlotsCondsData.CleanUp(self)
  496. for tmp in self.cats_rasts_conds.itervalues():
  497. grass.try_remove(tmp)
  498. for tmp in self.cats_rasts.itervalues():
  499. RunCommand("g.remove",
  500. rast=tmp,
  501. getErrorMsg=True)
  502. self.cats_rasts = {}
  503. self.cats_rasts_conds = {}
  504. def GetCatRast(self, cat_id):
  505. if self.cats_rasts.has_key(cat_id):
  506. return self.cats_rasts[cat_id]
  507. return None
  508. def GetCatRastCond(self, cat_id):
  509. return self.cats_rasts_conds[cat_id]
  510. def GetCatsRastsConds(self):
  511. max_cat_id = max(self.cats_rasts_conds.keys())
  512. cats_rasts_conds = [''] * (max_cat_id + 1)
  513. for i_cat_id, i_rast in self.cats_rasts_conds.iteritems():
  514. cats_rasts_conds[i_cat_id] = i_rast
  515. return cats_rasts_conds
  516. def GetCatsRasts(self):
  517. max_cat_id = max(self.cats_rasts.keys())
  518. cats_rasts = [''] * (max_cat_id + 1)
  519. for i_cat_id, i_rast in self.cats_rasts.iteritems():
  520. cats_rasts[i_cat_id] = i_rast
  521. return cats_rasts
  522. # not used, using iclass_perimeter algorithm instead of scipy convolve2d
  523. """
  524. def RasterizePolygon(pol, height, min_h, width, min_w):
  525. # Joe Kington
  526. # http://stackoverflow.com/questions/3654289/scipy-create-2d-polygon-mask
  527. #poly_verts = [(1,1), (1,4), (4,4),(4,1), (1,1)]
  528. nx = width
  529. ny = height
  530. x, y = np.meshgrid(np.arange(-0.5 + min_w, nx + 0.5 + min_w, dtype=float),
  531. np.arange(-0.5 + min_h, ny + 0.5 + min_h, dtype=float))
  532. x, y = x.flatten(), y.flatten()
  533. points = np.vstack((x,y)).T
  534. p = Path(pol)
  535. grid = p.contains_points(points)
  536. grid = grid.reshape((ny + 1, nx + 1))
  537. raster = np.zeros((height, width), dtype=np.uint8)#TODO bool
  538. #TODO shift by 0.5
  539. B = np.ones((2,2))/4
  540. raster = convolve2d(grid, B, 'valid')
  541. return raster
  542. """
  543. def idScattToidBands(scatt_id, n_bands):
  544. """!Get bands ids from scatter plot id."""
  545. n_b1 = n_bands - 1
  546. band_1 = (int) ((2 * n_b1 + 1 - sqrt(((2 * n_b1 + 1) * (2 * n_b1 + 1) - 8 * scatt_id))) / 2)
  547. band_2 = scatt_id - (band_1 * (2 * n_b1 + 1) - band_1 * band_1) / 2 + band_1 + 1
  548. return band_1, band_2
  549. def idBandsToidScatt(band_1_id, band_2_id, n_bands):
  550. """!Get scatter plot id from band ids."""
  551. if band_2_id < band_1_id:
  552. tmp = band_1_id
  553. band_1_id = band_2_id
  554. band_2_id = tmp
  555. n_b1 = n_bands - 1
  556. scatt_id = (band_1_id * (2 * n_b1 + 1) - band_1_id * band_1_id) / 2 + band_2_id - band_1_id - 1
  557. return scatt_id
  558. def GetRegion():
  559. ret, region, msg = RunCommand("g.region",
  560. flags = "gp",
  561. getErrorMsg = True,
  562. read = True)
  563. if ret != 0:
  564. raise GException("g.region failed:\n%s" % msg)
  565. return _parseRegion(region)
  566. def _parseRegion(region_str):
  567. region = {}
  568. region_str = region_str.splitlines()
  569. for param in region_str:
  570. k, v = param.split("=")
  571. if k in ["rows", "cols", "cells"]:
  572. v = int(v)
  573. else:
  574. v = float(v)
  575. region[k] = v
  576. return region
  577. def GetRasterInfo(rast):
  578. ret, out, msg = RunCommand("r.info",
  579. map = rast,
  580. flags = "rg",
  581. getErrorMsg = True,
  582. read = True)
  583. if ret != 0:
  584. raise GException("r.info failed:\n%s" % msg)
  585. out = out.split("\n")
  586. raster_info = {}
  587. for b in out:
  588. if not b.strip():
  589. continue
  590. k, v = b.split("=")
  591. if k == "datatype":
  592. if v != "CELL":
  593. return None
  594. pass
  595. elif k in ['rows', 'cols', 'cells', 'min', 'max']:
  596. v = int(v)
  597. else:
  598. v = float(v)
  599. raster_info[k] = v
  600. raster_info['range'] = raster_info['max'] - raster_info['min'] + 1
  601. return raster_info