__init__.py 30 KB

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