__init__.py 30 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869
  1. # -*- coding: utf-8 -*-
  2. from os.path import join, exists
  3. import grass.lib.gis as libgis
  4. libgis.G_gisinit('')
  5. import grass.lib.vector as libvect
  6. import ctypes
  7. #
  8. # import pygrass modules
  9. #
  10. from grass.pygrass.vector.vector_type import VTYPE
  11. from grass.pygrass.errors import GrassError, must_be_open
  12. from grass.pygrass.gis import Location
  13. from grass.pygrass.vector.geometry import GEOOBJ as _GEOOBJ
  14. from grass.pygrass.vector.geometry import read_line, read_next_line
  15. from grass.pygrass.vector.geometry import Area as _Area
  16. from grass.pygrass.vector.abstract import Info
  17. from grass.pygrass.vector.basic import Bbox, Cats, Ilist
  18. _NUMOF = {"areas": libvect.Vect_get_num_areas,
  19. "dblinks": libvect.Vect_get_num_dblinks,
  20. "faces": libvect.Vect_get_num_faces,
  21. "holes": libvect.Vect_get_num_holes,
  22. "islands": libvect.Vect_get_num_islands,
  23. "kernels": libvect.Vect_get_num_kernels,
  24. "points": libvect.Vect_get_num_lines,
  25. "lines": libvect.Vect_get_num_lines,
  26. "nodes": libvect.Vect_get_num_nodes,
  27. "updated_lines": libvect.Vect_get_num_updated_lines,
  28. "updated_nodes": libvect.Vect_get_num_updated_nodes,
  29. "volumes": libvect.Vect_get_num_volumes}
  30. # For test purposes
  31. test_vector_name = "vector_doctest_map"
  32. #=============================================
  33. # VECTOR
  34. #=============================================
  35. class Vector(Info):
  36. """Vector class is the grass vector format without topology
  37. >>> from grass.pygrass.vector import Vector
  38. >>> test_vect = Vector(test_vector_name)
  39. >>> test_vect.is_open()
  40. False
  41. >>> test_vect.mapset
  42. ''
  43. >>> test_vect.exist()
  44. True
  45. >>> test_vect.overwrite
  46. False
  47. """
  48. def __init__(self, name, mapset='', *args, **kwargs):
  49. # Set map name and mapset
  50. super(Vector, self).__init__(name, mapset, *args, **kwargs)
  51. self._topo_level = 1
  52. self._class_name = 'Vector'
  53. self.overwrite = False
  54. self._cats = []
  55. def __repr__(self):
  56. if self.exist():
  57. return "%s(%r, %r)" % (self._class_name, self.name, self.mapset)
  58. else:
  59. return "%s(%r)" % (self._class_name, self.name)
  60. def __iter__(self):
  61. """::
  62. >>> test_vect = Vector(test_vector_name)
  63. >>> test_vect.open(mode='r')
  64. >>> features = [feature for feature in test_vect]
  65. >>> features[:3]
  66. [Point(10.000000, 6.000000), Point(12.000000, 6.000000), Point(14.000000, 6.000000)]
  67. >>> test_vect.close()
  68. ..
  69. """
  70. #return (self.read(f_id) for f_id in xrange(self.num_of_features()))
  71. return self
  72. @must_be_open
  73. def next(self):
  74. """::
  75. >>> test_vect = Vector(test_vector_name)
  76. >>> test_vect.open(mode='r')
  77. >>> test_vect.next()
  78. Point(10.000000, 6.000000)
  79. >>> test_vect.next()
  80. Point(12.000000, 6.000000)
  81. >>> test_vect.close()
  82. ..
  83. """
  84. return read_next_line(self.c_mapinfo, self.table, self.writeable,
  85. is2D=not self.is_3D())
  86. @must_be_open
  87. def rewind(self):
  88. """Rewind vector map to cause reads to start at beginning."""
  89. if libvect.Vect_rewind(self.c_mapinfo) == -1:
  90. raise GrassError("Vect_rewind raise an error.")
  91. @must_be_open
  92. def write(self, geo_obj, cat=None, attrs=None):
  93. """Write geometry features and attributes.
  94. :param geo_obj: a geometry grass object define in
  95. grass.pygrass.vector.geometry
  96. :type geo_obj: geometry GRASS object
  97. :param attrs: a list with the values that will be insert in the
  98. attribute table.
  99. :type attrs: list
  100. :param cat: The category of the geometry feature, otherwise the
  101. c_cats attribute of the geometry object will be used.
  102. :type cat: integer
  103. Open a new vector map ::
  104. >>> new = VectorTopo('newvect')
  105. >>> new.exist()
  106. False
  107. define the new columns of the attribute table ::
  108. >>> cols = [(u'cat', 'INTEGER PRIMARY KEY'),
  109. ... (u'name', 'TEXT')]
  110. open the vector map in write mode
  111. >>> new.open('w', tab_name='newvect', tab_cols=cols)
  112. import a geometry feature ::
  113. >>> from grass.pygrass.vector.geometry import Point
  114. create two points ::
  115. >>> point0 = Point(0, 0)
  116. >>> point1 = Point(1, 1)
  117. then write the two points on the map, with ::
  118. >>> new.write(point0, cat=1, attrs=('pub',))
  119. >>> new.write(point1, cat=2, attrs=('resturant',))
  120. commit the db changes ::
  121. >>> new.table.conn.commit()
  122. >>> new.table.execute().fetchall()
  123. [(1, u'pub'), (2, u'resturant')]
  124. close the vector map ::
  125. >>> new.close()
  126. >>> new.exist()
  127. True
  128. then play with the map ::
  129. >>> new.open(mode='r')
  130. >>> new.read(1)
  131. Point(0.000000, 0.000000)
  132. >>> new.read(2)
  133. Point(1.000000, 1.000000)
  134. >>> new.read(1).attrs['name']
  135. u'pub'
  136. >>> new.read(2).attrs['name']
  137. u'resturant'
  138. >>> new.close()
  139. >>> new.remove()
  140. """
  141. self.n_lines += 1
  142. if attrs and cat is None:
  143. cat = (self._cats[-1] if self._cats else 0) + 1
  144. if cat is not None and cat not in self._cats:
  145. self._cats.append(cat)
  146. if self.table is not None and attrs is not None:
  147. attr = [cat, ]
  148. attr.extend(attrs)
  149. cur = self.table.conn.cursor()
  150. cur.execute(self.table.columns.insert_str, attr)
  151. cur.close()
  152. if cat is not None:
  153. cats = Cats(geo_obj.c_cats)
  154. cats.reset()
  155. cats.set(cat, self.layer)
  156. if geo_obj.gtype == _Area.gtype:
  157. result = self._write_area(geo_obj)
  158. result = libvect.Vect_write_line(self.c_mapinfo, geo_obj.gtype,
  159. geo_obj.c_points, geo_obj.c_cats)
  160. if result == -1:
  161. raise GrassError("Not able to write the vector feature.")
  162. if self._topo_level == 2:
  163. # return new feature id (on level 2)
  164. geo_obj.id = result
  165. else:
  166. # return offset into file where the feature starts (on level 1)
  167. geo_obj.offset = result
  168. @must_be_open
  169. def has_color_table(self):
  170. """Return if vector has color table associated in file system;
  171. Color table stored in the vector's attribute table well be not checked
  172. >>> test_vect = Vector(test_vector_name)
  173. >>> test_vect.open(mode='r')
  174. >>> test_vect.has_color_table()
  175. False
  176. >>> test_vect.close()
  177. >>> from grass.pygrass.utils import copy, remove
  178. >>> copy(test_vector_name,'mytest_vect','vect')
  179. >>> from grass.pygrass.modules.shortcuts import vector as v
  180. >>> v.colors(map='mytest_vect', color='population', column='value')
  181. Module('v.colors')
  182. >>> mytest_vect = Vector('mytest_vect')
  183. >>> mytest_vect.open(mode='r')
  184. >>> mytest_vect.has_color_table()
  185. True
  186. >>> mytest_vect.close()
  187. >>> remove('mytest_vect', 'vect')
  188. """
  189. loc = Location()
  190. path = join(loc.path(), self.mapset, 'vector', self.name, 'colr')
  191. return True if exists(path) else False
  192. #=============================================
  193. # VECTOR WITH TOPOLOGY
  194. #=============================================
  195. class VectorTopo(Vector):
  196. """Vector class with the support of the GRASS topology.
  197. Open a vector map using the *with statement*: ::
  198. >>> with VectorTopo(test_vector_name, mode='r') as test_vect:
  199. ... for feature in test_vect[:7]:
  200. ... print feature.attrs['name']
  201. ...
  202. point
  203. point
  204. point
  205. line
  206. line
  207. line
  208. >>> test_vect.is_open()
  209. False
  210. ..
  211. """
  212. def __init__(self, name, mapset='', *args, **kwargs):
  213. super(VectorTopo, self).__init__(name, mapset, *args, **kwargs)
  214. self._topo_level = 2
  215. self._class_name = 'VectorTopo'
  216. def __len__(self):
  217. return libvect.Vect_get_num_lines(self.c_mapinfo)
  218. def __getitem__(self, key):
  219. """::
  220. >>> test_vect = VectorTopo(test_vector_name)
  221. >>> test_vect.open(mode='r')
  222. >>> test_vect[:4]
  223. [Point(10.000000, 6.000000), Point(12.000000, 6.000000), Point(14.000000, 6.000000)]
  224. >>> test_vect.close()
  225. ..
  226. """
  227. if isinstance(key, slice):
  228. return [self.read(indx)
  229. for indx in range(key.start if key.start else 1,
  230. key.stop if key.stop else len(self),
  231. key.step if key.step else 1)]
  232. elif isinstance(key, int):
  233. return self.read(key)
  234. else:
  235. raise ValueError("Invalid argument type: %r." % key)
  236. @must_be_open
  237. def num_primitive_of(self, primitive):
  238. """Return the number of primitive
  239. :param primitive: the name of primitive to query; the supported values are:
  240. * *boundary*,
  241. * *centroid*,
  242. * *face*,
  243. * *kernel*,
  244. * *line*,
  245. * *point*
  246. * *area*
  247. * *volume*
  248. :type primitive: str
  249. ::
  250. >>> test_vect = VectorTopo(test_vector_name)
  251. >>> test_vect.open(mode='r')
  252. >>> test_vect.num_primitive_of('point')
  253. 3
  254. >>> test_vect.num_primitive_of('line')
  255. 3
  256. >>> test_vect.num_primitive_of('centroid')
  257. 4
  258. >>> test_vect.num_primitive_of('boundary')
  259. 11
  260. >>> test_vect.close()
  261. ..
  262. """
  263. return libvect.Vect_get_num_primitives(self.c_mapinfo,
  264. VTYPE[primitive])
  265. @must_be_open
  266. def number_of(self, vtype):
  267. """Return the number of the choosen element type
  268. :param vtype: the name of type to query; the supported values are:
  269. *areas*, *dblinks*, *faces*, *holes*, *islands*,
  270. *kernels*, *line_points*, *lines*, *nodes*, *points*,
  271. *update_lines*, *update_nodes*, *volumes*
  272. :type vtype: str
  273. >>> test_vect = VectorTopo(test_vector_name)
  274. >>> test_vect.open(mode='r')
  275. >>> test_vect.number_of("areas")
  276. 4
  277. >>> test_vect.number_of("islands")
  278. 2
  279. >>> test_vect.number_of("holes")
  280. 0
  281. >>> test_vect.number_of("lines")
  282. 21
  283. >>> test_vect.number_of("nodes")
  284. 15
  285. >>> test_vect.number_of("pizza")
  286. ... # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
  287. Traceback (most recent call last):
  288. ...
  289. ValueError: vtype not supported, use one of: 'areas', ...
  290. >>> test_vect.close()
  291. ..
  292. """
  293. if vtype in _NUMOF.keys():
  294. return _NUMOF[vtype](self.c_mapinfo)
  295. else:
  296. keys = "', '".join(sorted(_NUMOF.keys()))
  297. raise ValueError("vtype not supported, use one of: '%s'" % keys)
  298. @must_be_open
  299. def num_primitives(self):
  300. """Return dictionary with the number of all primitives
  301. """
  302. output = {}
  303. for prim in VTYPE.keys():
  304. output[prim] = self.num_primitive_of(prim)
  305. return output
  306. @must_be_open
  307. def viter(self, vtype, idonly=False):
  308. """Return an iterator of vector features
  309. :param vtype: the name of type to query; the supported values are:
  310. *areas*, *dblinks*, *faces*, *holes*, *islands*,
  311. *kernels*, *line_points*, *lines*, *nodes*, *points*,
  312. *update_lines*, *update_nodes*, *volumes*
  313. :type vtype: str
  314. :param idonly: variable to return only the id of features instead of
  315. full features
  316. :type idonly: bool
  317. >>> test_vect = VectorTopo(test_vector_name, mode='r')
  318. >>> test_vect.open(mode='r')
  319. >>> areas = [area for area in test_vect.viter('areas')]
  320. >>> areas[:3]
  321. [Area(1), Area(2), Area(3)]
  322. to sort the result in a efficient way, use: ::
  323. >>> from operator import methodcaller as method
  324. >>> areas.sort(key=method('area'), reverse=True) # sort the list
  325. >>> for area in areas[:3]:
  326. ... print area, area.area()
  327. Area(1) 12.0
  328. Area(2) 8.0
  329. Area(4) 8.0
  330. >>> areas = [area for area in test_vect.viter('areas')]
  331. >>> for area in areas:
  332. ... print(area.centroid().cat)
  333. 3
  334. 3
  335. 3
  336. 3
  337. >>> test_vect.close()
  338. """
  339. if vtype in _GEOOBJ.keys():
  340. if _GEOOBJ[vtype] is not None:
  341. ids = (indx for indx in range(1, self.number_of(vtype) + 1))
  342. if idonly:
  343. return ids
  344. return (_GEOOBJ[vtype](v_id=indx, c_mapinfo=self.c_mapinfo,
  345. table=self.table,
  346. writeable=self.writeable)
  347. for indx in ids)
  348. else:
  349. keys = "', '".join(sorted(_GEOOBJ.keys()))
  350. raise ValueError("vtype not supported, use one of: '%s'" % keys)
  351. @must_be_open
  352. def rewind(self):
  353. """Rewind vector map to cause reads to start at beginning. ::
  354. >>> test_vect = VectorTopo(test_vector_name)
  355. >>> test_vect.open(mode='r')
  356. >>> test_vect.next()
  357. Point(10.000000, 6.000000)
  358. >>> test_vect.next()
  359. Point(12.000000, 6.000000)
  360. >>> test_vect.next()
  361. Point(14.000000, 6.000000)
  362. >>> test_vect.rewind()
  363. >>> test_vect.next()
  364. Point(10.000000, 6.000000)
  365. >>> test_vect.close()
  366. ..
  367. """
  368. libvect.Vect_rewind(self.c_mapinfo)
  369. @must_be_open
  370. def cat(self, cat_id, vtype, layer=None, generator=False, geo=None):
  371. """Return the geometry features with category == cat_id.
  372. :param cat_id: the category number
  373. :type cat_id: int
  374. :param vtype: the type of geometry feature that we are looking for
  375. :type vtype: str
  376. :param layer: the layer number that will be used
  377. :type layer: int
  378. :param generator: if True return a generator otherwise it return a
  379. list of features
  380. :type generator: bool
  381. """
  382. if geo is None and vtype not in _GEOOBJ:
  383. keys = "', '".join(sorted(_GEOOBJ.keys()))
  384. raise ValueError("vtype not supported, use one of: '%s'" % keys)
  385. Obj = _GEOOBJ[vtype] if geo is None else geo
  386. ilist = Ilist()
  387. libvect.Vect_cidx_find_all(self.c_mapinfo,
  388. layer if layer else self.layer,
  389. Obj.gtype, cat_id, ilist.c_ilist)
  390. is2D = not self.is_3D()
  391. if generator:
  392. return (read_line(feature_id=v_id, c_mapinfo=self.c_mapinfo,
  393. table=self.table, writeable=self.writeable,
  394. is2D=is2D)
  395. for v_id in ilist)
  396. else:
  397. return [read_line(feature_id=v_id, c_mapinfo=self.c_mapinfo,
  398. table=self.table, writeable=self.writeable,
  399. is2D=is2D)
  400. for v_id in ilist]
  401. @must_be_open
  402. def read(self, feature_id):
  403. """Return a geometry object given the feature id.
  404. :param int feature_id: the id of feature to obtain
  405. >>> test_vect = VectorTopo(test_vector_name)
  406. >>> test_vect.open(mode='r')
  407. >>> feature1 = test_vect.read(0) #doctest: +ELLIPSIS
  408. Traceback (most recent call last):
  409. ...
  410. ValueError: The index must be >0, 0 given.
  411. >>> feature1 = test_vect.read(5)
  412. >>> feature1
  413. Line([Point(12.000000, 4.000000), Point(12.000000, 2.000000), Point(12.000000, 0.000000)])
  414. >>> feature1.length()
  415. 4.0
  416. >>> test_vect.read(-1)
  417. Centoid(7.500000, 3.500000)
  418. >>> len(test_vect)
  419. 21
  420. >>> test_vect.read(21)
  421. Centoid(7.500000, 3.500000)
  422. >>> test_vect.read(22) #doctest: +ELLIPSIS
  423. Traceback (most recent call last):
  424. ...
  425. IndexError: Index out of range
  426. >>> test_vect.close()
  427. """
  428. return read_line(feature_id, self.c_mapinfo, self.table, self.writeable,
  429. is2D=not self.is_3D())
  430. @must_be_open
  431. def is_empty(self):
  432. """Return if a vector map is empty or not
  433. """
  434. primitives = self.num_primitives()
  435. output = True
  436. for v in primitives.values():
  437. if v != 0:
  438. output = False
  439. break
  440. return output
  441. @must_be_open
  442. def rewrite(self, line, geo_obj, attrs=None, **kargs):
  443. """Rewrite a geometry features
  444. """
  445. if self.table is not None and attrs:
  446. attr = [line, ]
  447. attr.extend(attrs)
  448. self.table.update(key=line, values=attr)
  449. elif self.table is None and attrs:
  450. print "Table for vector {name} does not exist, attributes not" \
  451. " loaded".format(name=self.name)
  452. libvect.Vect_cat_set(geo_obj.c_cats, self.layer, line)
  453. result = libvect.Vect_rewrite_line(self.c_mapinfo,
  454. line, geo_obj.gtype,
  455. geo_obj.c_points,
  456. geo_obj.c_cats)
  457. if result == -1:
  458. raise GrassError("Not able to write the vector feature.")
  459. # return offset into file where the feature starts
  460. geo_obj.offset = result
  461. @must_be_open
  462. def delete(self, feature_id):
  463. """Remove a feature by its id
  464. :param feature_id: the id of the feature
  465. :type feature_id: int
  466. """
  467. if libvect.Vect_rewrite_line(self.c_mapinfo, feature_id) == -1:
  468. raise GrassError("C funtion: Vect_rewrite_line.")
  469. @must_be_open
  470. def restore(self, geo_obj):
  471. if hasattr(geo_obj, 'offset'):
  472. if libvect.Vect_restore_line(self.c_mapinfo, geo_obj.offset,
  473. geo_obj.id) == -1:
  474. raise GrassError("C funtion: Vect_restore_line.")
  475. else:
  476. raise ValueError("The value have not an offset attribute.")
  477. @must_be_open
  478. def bbox(self):
  479. """Return the BBox of the vecor map
  480. """
  481. bbox = Bbox()
  482. if libvect.Vect_get_map_box(self.c_mapinfo, bbox.c_bbox) == 0:
  483. raise GrassError("I can not find the Bbox.")
  484. return bbox
  485. def close(self, build=True, release=True):
  486. """Close the VectorTopo map, if release is True, the memory
  487. occupied by spatial index is released"""
  488. if release:
  489. libvect.Vect_set_release_support(self.c_mapinfo)
  490. super(VectorTopo, self).close(build=build)
  491. @must_be_open
  492. def table_to_dict(self, where=None):
  493. """Return the attribute table as a dictionary with the category as keys
  494. The columns have the order of the self.table.columns.names() list.
  495. Examples
  496. >>> from grass.pygrass.vector import VectorTopo
  497. >>> from grass.pygrass.vector.basic import Bbox
  498. >>> test_vect = VectorTopo(test_vector_name)
  499. >>> test_vect.open('r')
  500. >>> test_vect.table_to_dict()
  501. {1: [1, u'point', 1.0], 2: [2, u'line', 2.0], 3: [3, u'centroid', 3.0]}
  502. >>> test_vect.table_to_dict(where="value > 2")
  503. {3: [3, u'centroid', 3.0]}
  504. >>> test_vect.table_to_dict(where="value > 0")
  505. {1: [1, u'point', 1.0], 2: [2, u'line', 2.0], 3: [3, u'centroid', 3.0]}
  506. >>> test_vect.table.filters.get_sql()
  507. u'SELECT cat,name,value FROM vector_doctest_map WHERE value > 0 ORDER BY cat;'
  508. """
  509. if self.table is not None:
  510. table_dict = {}
  511. # Get the category index
  512. cat_index = self.table.columns.names().index("cat")
  513. # Prepare a filter
  514. if where is not None:
  515. self.table.filters.where(where)
  516. self.table.filters.order_by("cat")
  517. self.table.filters.select(",".join(self.table.columns.names()))
  518. # Execute the query and fetch the result
  519. cur = self.table.execute()
  520. l = cur.fetchall()
  521. # Generate the dictionary
  522. for entry in l:
  523. table_dict[entry[cat_index]] = list(entry)
  524. return(table_dict)
  525. return None
  526. @must_be_open
  527. def features_to_wkb_list(self, bbox=None, feature_type="point", field=1):
  528. """Return all features of type point, line, boundary or centroid
  529. as a list of Well Known Binary representations (WKB)
  530. (id, cat, wkb) triplets located in a specific
  531. bounding box.
  532. :param bbox: The boundingbox to search for features,
  533. if bbox=None the boundingbox of the whole
  534. vector map layer is used
  535. :type bbox: grass.pygrass.vector.basic.Bbox
  536. :param feature_type: The type of feature that should be converted to
  537. the Well Known Binary (WKB) format. Supported are:
  538. 'point' -> libvect.GV_POINT 1
  539. 'line' -> libvect.GV_LINE 2
  540. 'boundary' -> libvect.GV_BOUNDARY 3
  541. 'centroid' -> libvect.GV_CENTROID 4
  542. :type type: string
  543. :param field: The category field
  544. :type field: integer
  545. :return: A list of triplets, or None if nothing was found
  546. The well known binary are stored in byte arrays.
  547. Examples:
  548. >>> from grass.pygrass.vector import VectorTopo
  549. >>> from grass.pygrass.vector.basic import Bbox
  550. >>> test_vect = VectorTopo(test_vector_name)
  551. >>> test_vect.open('r')
  552. >>> bbox = Bbox(north=20, south=-1, east=20, west=-1)
  553. >>> result = test_vect.features_to_wkb_list(bbox=bbox,
  554. ... feature_type="point")
  555. >>> len(result)
  556. 3
  557. >>> for entry in result:
  558. ... f_id, cat, wkb = entry
  559. ... print(f_id, cat, len(wkb))
  560. (1, 1, 21)
  561. (2, 1, 21)
  562. (3, 1, 21)
  563. >>> result = test_vect.features_to_wkb_list(bbox=None,
  564. ... feature_type="line")
  565. >>> len(result)
  566. 3
  567. >>> for entry in result:
  568. ... f_id, cat, wkb = entry
  569. ... print(f_id, cat, len(wkb))
  570. (4, 2, 57)
  571. (5, 2, 57)
  572. (6, 2, 57)
  573. >>> result = test_vect.features_to_wkb_list(bbox=bbox,
  574. ... feature_type="boundary")
  575. >>> len(result)
  576. 11
  577. >>> result = test_vect.features_to_wkb_list(bbox=None,
  578. ... feature_type="centroid")
  579. >>> len(result)
  580. 4
  581. >>> for entry in result:
  582. ... f_id, cat, wkb = entry
  583. ... print(f_id, cat, len(wkb))
  584. (19, 3, 21)
  585. (18, 3, 21)
  586. (20, 3, 21)
  587. (21, 3, 21)
  588. >>> result = test_vect.features_to_wkb_list(bbox=bbox,
  589. ... feature_type="blub")
  590. Traceback (most recent call last):
  591. ...
  592. GrassError: Unsupported feature type <blub>, supported are <point,line,boundary,centroid>
  593. >>> test_vect.close()
  594. """
  595. supported = ['point', 'line', 'boundary', 'centroid']
  596. if feature_type.lower() not in supported:
  597. raise GrassError("Unsupported feature type <%s>, "\
  598. "supported are <%s>"%(feature_type,
  599. ",".join(supported)))
  600. if bbox is None:
  601. bbox = self.bbox()
  602. bboxlist = self.find_by_bbox.geos(bbox, type=feature_type.lower(),
  603. bboxlist_only = True)
  604. if bboxlist is not None and len(bboxlist) > 0:
  605. l = []
  606. line_p = libvect.line_pnts()
  607. line_c = libvect.line_cats()
  608. size = ctypes.c_size_t()
  609. cat = ctypes.c_int()
  610. error = ctypes.c_int()
  611. for f_id in bboxlist.ids:
  612. barray = libvect.Vect_read_line_to_wkb(self.c_mapinfo,
  613. ctypes.byref(line_p),
  614. ctypes.byref(line_c),
  615. f_id,
  616. ctypes.byref(size),
  617. ctypes.byref(error))
  618. if not barray:
  619. if error == -1:
  620. raise GrassError(_("Unable to read line of feature %i"%(f_id)))
  621. if error == -2:
  622. print("Empty feature %i"%(f_id))
  623. continue
  624. ok = libvect.Vect_cat_get(ctypes.byref(line_c), field,
  625. ctypes.byref(cat))
  626. if ok < 1:
  627. pcat = None
  628. else:
  629. pcat = cat.value
  630. l.append((f_id, pcat, ctypes.string_at(barray, size.value)))
  631. libgis.G_free(barray)
  632. return l
  633. return None
  634. @must_be_open
  635. def areas_to_wkb_list(self, bbox=None, field=1):
  636. """Return all features of type point, line, boundary or centroid
  637. as a list of Well Known Binary representations (WKB)
  638. (id, cat, wkb) triplets located in a specific
  639. bounding box.
  640. :param bbox: The boundingbox to search for features,
  641. if bbox=None the boundingbox of the whole
  642. vector map layer is used
  643. :type bbox: grass.pygrass.vector.basic.Bbox
  644. :param field: The centroid category field
  645. :type field: integer
  646. :return: A list of triplets, or None if nothing was found
  647. The well known binary are stored in byte arrays.
  648. Examples:
  649. >>> from grass.pygrass.vector import VectorTopo
  650. >>> from grass.pygrass.vector.basic import Bbox
  651. >>> test_vect = VectorTopo(test_vector_name)
  652. >>> test_vect.open('r')
  653. >>> bbox = Bbox(north=20, south=-1, east=20, west=-1)
  654. >>> result = test_vect.areas_to_wkb_list(bbox=bbox)
  655. >>> len(result)
  656. 4
  657. >>> for entry in result:
  658. ... a_id, cat, wkb = entry
  659. ... print(a_id, cat, len(wkb))
  660. (1, 3, 225)
  661. (2, 3, 141)
  662. (3, 3, 93)
  663. (4, 3, 141)
  664. >>> result = test_vect.areas_to_wkb_list()
  665. >>> len(result)
  666. 4
  667. >>> for entry in result:
  668. ... a_id, cat, wkb = entry
  669. ... print(a_id, cat, len(wkb))
  670. (1, 3, 225)
  671. (2, 3, 141)
  672. (3, 3, 93)
  673. (4, 3, 141)
  674. >>> test_vect.close()
  675. """
  676. if bbox is None:
  677. bbox = self.bbox()
  678. bboxlist = self.find_by_bbox.areas(bbox, bboxlist_only = True)
  679. if bboxlist is not None and len(bboxlist) > 0:
  680. l = []
  681. line_c = libvect.line_cats()
  682. size = ctypes.c_size_t()
  683. cat = ctypes.c_int()
  684. for a_id in bboxlist.ids:
  685. barray = libvect.Vect_read_area_to_wkb(self.c_mapinfo,
  686. a_id,
  687. ctypes.byref(size))
  688. if not barray:
  689. raise GrassError(_("Unable to read area with id %i"%(a_id)))
  690. pcat = None
  691. c_ok = libvect.Vect_get_area_cats(self.c_mapinfo, a_id,
  692. ctypes.byref(line_c))
  693. if c_ok == 0: # Centroid found
  694. ok = libvect.Vect_cat_get(ctypes.byref(line_c), field,
  695. ctypes.byref(cat))
  696. if ok > 0:
  697. pcat = cat.value
  698. l.append((a_id, pcat, ctypes.string_at(barray, size.value)))
  699. libgis.G_free(barray)
  700. return l
  701. return None
  702. if __name__ == "__main__":
  703. import doctest
  704. from grass.pygrass import utils
  705. utils.create_test_vector_map(test_vector_name)
  706. doctest.testmod()
  707. """Remove the generated vector map, if exist"""
  708. from grass.pygrass.utils import get_mapset_vector
  709. from grass.script.core import run_command
  710. mset = get_mapset_vector(test_vector_name, mapset='')
  711. if mset:
  712. run_command("g.remove", flags='f', type='vector', name=test_vector_name)