geometry.py 52 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586
  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Wed Jul 18 10:46:25 2012
  4. @author: pietro
  5. """
  6. import ctypes
  7. import re
  8. import numpy as np
  9. import grass.lib.gis as libgis
  10. import grass.lib.vector as libvect
  11. from grass.pygrass.errors import GrassError
  12. from grass.pygrass.vector.basic import Ilist, Bbox, Cats
  13. from grass.pygrass.vector import sql
  14. WKT = {'POINT\((.*)\)': 'point', # 'POINT\(\s*([+-]*\d+\.*\d*)+\s*\)'
  15. 'LINESTRING\((.*)\)': 'line'}
  16. def read_WKT(string):
  17. """Read the string and return a geometry object
  18. WKT:
  19. POINT(0 0)
  20. LINESTRING(0 0,1 1,1 2)
  21. POLYGON((0 0,4 0,4 4,0 4,0 0),(1 1, 2 1, 2 2, 1 2,1 1))
  22. MULTIPOINT(0 0,1 2)
  23. MULTILINESTRING((0 0,1 1,1 2),(2 3,3 2,5 4))
  24. MULTIPOLYGON(((0 0,4 0,4 4,0 4,0 0),(1 1,2 1,2 2,1 2,1 1)),
  25. ((-1 -1,-1 -2,-2 -2,-2 -1,-1 -1)))
  26. GEOMETRYCOLLECTION(POINT(2 3),LINESTRING(2 3,3 4))
  27. EWKT:
  28. POINT(0 0 0) -- XYZ
  29. SRID=32632;POINT(0 0) -- XY with SRID
  30. POINTM(0 0 0) -- XYM
  31. POINT(0 0 0 0) -- XYZM
  32. SRID=4326;MULTIPOINTM(0 0 0,1 2 1) -- XYM with SRID
  33. MULTILINESTRING((0 0 0,1 1 0,1 2 1),(2 3 1,3 2 1,5 4 1))
  34. POLYGON((0 0 0,4 0 0,4 4 0,0 4 0,0 0 0),(1 1 0,2 1 0,2 2 0,1 2 0,1 1 0))
  35. MULTIPOLYGON(((0 0 0,4 0 0,4 4 0,0 4 0,0 0 0),
  36. (1 1 0,2 1 0,2 2 0,1 2 0,1 1 0)),
  37. ((-1 -1 0,-1 -2 0,-2 -2 0,-2 -1 0,-1 -1 0)))
  38. GEOMETRYCOLLECTIONM( POINTM(2 3 9), LINESTRINGM(2 3 4, 3 4 5) )
  39. MULTICURVE( (0 0, 5 5), CIRCULARSTRING(4 0, 4 4, 8 4) )
  40. POLYHEDRALSURFACE( ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)),
  41. ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)),
  42. ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)),
  43. ((1 1 0, 1 1 1, 1 0 1, 1 0 0, 1 1 0)),
  44. ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0)),
  45. ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)) )
  46. TRIANGLE ((0 0, 0 9, 9 0, 0 0))
  47. TIN( ((0 0 0, 0 0 1, 0 1 0, 0 0 0)), ((0 0 0, 0 1 0, 1 1 0, 0 0 0)) )
  48. """
  49. for regexp, obj in WKT.items():
  50. if re.match(regexp, string):
  51. geo = 10
  52. return obj(geo)
  53. def read_WKB(buff):
  54. """Read the binary buffer and return a geometry object"""
  55. pass
  56. def intersects(lineA, lineB, with_z=False):
  57. """Return a list of points ::
  58. >>> lineA = Line([(0, 0), (4, 0)])
  59. >>> lineB = Line([(2, 2), (2, -2)])
  60. >>> intersects(lineA, lineB)
  61. Line([Point(2.000000, 0.000000)])
  62. """
  63. line = Line()
  64. if libvect.Vect_line_get_intersections(lineA.c_points, lineB.c_points,
  65. line.c_points, int(with_z)):
  66. return line
  67. else:
  68. return []
  69. #=============================================
  70. # GEOMETRY
  71. #=============================================
  72. def get_xyz(pnt):
  73. """Return a tuple with: x, y, z. ::
  74. >>> pnt = Point(0, 0)
  75. >>> get_xyz(pnt)
  76. (0.0, 0.0, 0.0)
  77. >>> get_xyz((1, 1))
  78. (1, 1, 0.0)
  79. >>> get_xyz((1, 1, 2))
  80. (1, 1, 2)
  81. >>> get_xyz((1, 1, 2, 2)) #doctest: +ELLIPSIS
  82. Traceback (most recent call last):
  83. ...
  84. ValueError: The the format of the point is not supported: (1, 1, 2, 2)
  85. ..
  86. """
  87. if isinstance(pnt, Point):
  88. if pnt.is2D:
  89. x, y = pnt.x, pnt.y
  90. z = 0.
  91. else:
  92. x, y, z = pnt.x, pnt.y, pnt.z
  93. else:
  94. if len(pnt) == 2:
  95. x, y = pnt
  96. z = 0.
  97. elif len(pnt) == 3:
  98. x, y, z = pnt
  99. else:
  100. str_error = "The the format of the point is not supported: {0!r}"
  101. raise ValueError(str_error.format(pnt))
  102. return x, y, z
  103. class Attrs(object):
  104. def __init__(self, cat, table, writable=False):
  105. self._cat = None
  106. self.cond = ''
  107. self.table = table
  108. self.cat = cat
  109. self.writable = writable
  110. def _get_cat(self):
  111. return self._cat
  112. def _set_cat(self, value):
  113. self._cat = value
  114. if value:
  115. # update condition
  116. self.cond = "%s=%d" % (self.table.key, value)
  117. cat = property(fget=_get_cat, fset=_set_cat,
  118. doc="Set and obtain cat value")
  119. def __getitem__(self, key):
  120. """Return the value stored in the attribute table. ::
  121. >>> from grass.pygrass.vector import VectorTopo
  122. >>> schools = VectorTopo('schools')
  123. >>> schools.open('r')
  124. >>> school = schools[1]
  125. >>> attrs = Attrs(school.cat, schools.table)
  126. >>> attrs['TAG']
  127. u'568'
  128. """
  129. #SELECT {cols} FROM {tname} WHERE {condition};
  130. try:
  131. cur = self.table.execute(sql.SELECT_WHERE.format(cols=key,
  132. tname=self.table.name,
  133. condition=self.cond))
  134. except:
  135. import ipdb; ipdb.set_trace()
  136. results = cur.fetchone()
  137. if results is not None:
  138. return results[0] if len(results) == 1 else results
  139. def __setitem__(self, key, value):
  140. """Set value of a given column of a table attribute. ::
  141. >>> from grass.pygrass.vector import VectorTopo
  142. >>> from grass.pygrass.functions import copy, remove
  143. >>> copy('schools', 'myschools', 'vect')
  144. >>> schools = VectorTopo('myschools')
  145. >>> schools.open('r')
  146. >>> school = schools[1]
  147. >>> attrs = Attrs(school.cat, schools.table, True)
  148. >>> attrs['TAG'] = 'New Label'
  149. >>> attrs['TAG']
  150. u'New Label'
  151. >>> attrs.table.conn.close()
  152. >>> remove('myschools','vect')
  153. """
  154. if self.writable:
  155. #UPDATE {tname} SET {new_col} = {old_col} WHERE {condition}
  156. values = '%s=%r' % (key, value)
  157. self.table.execute(sql.UPDATE_WHERE.format(tname=self.table.name,
  158. values=values,
  159. condition=self.cond))
  160. #self.table.conn.commit()
  161. else:
  162. str_err = "You can only read the attributes if the map is \
  163. in another mapset"
  164. raise GrassError(str_err)
  165. def __dict__(self):
  166. """Return a dict of the attribute table row."""
  167. dic = {}
  168. for key, val in zip(self.keys(), self.values()):
  169. dic[key] = val
  170. return dic
  171. def values(self):
  172. """Return the values of the attribute table row.
  173. >>> from grass.pygrass.vector import VectorTopo
  174. >>> schools = VectorTopo('schools')
  175. >>> schools.open('r')
  176. >>> school = schools[1]
  177. >>> attrs = Attrs(school.cat, schools.table)
  178. >>> attrs.values() # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
  179. (1,
  180. ...
  181. None)
  182. """
  183. #SELECT {cols} FROM {tname} WHERE {condition}
  184. cur = self.table.execute(sql.SELECT_WHERE.format(cols='*',
  185. tname=self.table.name,
  186. condition=self.cond))
  187. return cur.fetchone()
  188. def keys(self):
  189. """Return the column name of the attribute table.
  190. >>> from grass.pygrass.vector import VectorTopo
  191. >>> schools = VectorTopo('schools')
  192. >>> schools.open('r')
  193. >>> school = schools[1]
  194. >>> attrs = Attrs(school.cat, schools.table)
  195. >>> attrs.keys() # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
  196. [u'cat',
  197. ...
  198. u'NOTES']
  199. """
  200. return self.table.columns.names()
  201. def commit(self):
  202. """Save the changes"""
  203. self.table.conn.commit()
  204. class Geo(object):
  205. """
  206. >>> geo0 = Geo()
  207. >>> points = ctypes.pointer(libvect.line_pnts())
  208. >>> cats = ctypes.pointer(libvect.line_cats())
  209. >>> geo1 = Geo(c_points=points, c_cats=cats)
  210. """
  211. gtype = None
  212. def __init__(self, v_id=None, c_mapinfo=None, c_points=None, c_cats=None,
  213. table=None, writable=False, is2D=True):
  214. self.id = v_id # vector id
  215. self.c_mapinfo = c_mapinfo
  216. self.is2D = is2D
  217. read = False
  218. # set c_points
  219. if c_points is None:
  220. self.c_points = ctypes.pointer(libvect.line_pnts())
  221. read = True
  222. else:
  223. self.c_points = c_points
  224. # set c_cats
  225. if c_cats is None:
  226. self.c_cats = ctypes.pointer(libvect.line_cats())
  227. read = True
  228. else:
  229. self.c_cats = c_cats
  230. # set the attributes
  231. self.attrs = None
  232. if table is not None:
  233. self.attrs = Attrs(self.cat, table, writable)
  234. if self.id is not None and self.c_mapinfo is not None and read:
  235. self.read()
  236. @property
  237. def cat(self):
  238. if self.c_cats.contents.cat:
  239. return self.c_cats.contents.cat.contents.value
  240. def is_with_topology(self):
  241. if self.c_mapinfo is not None:
  242. return self.c_mapinfo.contents.level == 2
  243. else:
  244. return False
  245. def read(self):
  246. """Read and set the coordinates of the centroid from the vector map,
  247. using the centroid_id and calling the Vect_read_line C function"""
  248. ftype, c_points, c_cats = c_read_line(self.id, self.c_mapinfo,
  249. self.c_points, self.c_cats)
  250. class Point(Geo):
  251. """Instantiate a Point object that could be 2 or 3D, default
  252. parameters are 0.
  253. ::
  254. >>> pnt = Point()
  255. >>> pnt.x
  256. 0.0
  257. >>> pnt.y
  258. 0.0
  259. >>> pnt.z
  260. >>> pnt.is2D
  261. True
  262. >>> pnt
  263. Point(0.000000, 0.000000)
  264. >>> pnt.z = 0
  265. >>> pnt.is2D
  266. False
  267. >>> pnt
  268. Point(0.000000, 0.000000, 0.000000)
  269. >>> print(pnt)
  270. POINT(0.000000 0.000000 0.000000)
  271. ..
  272. """
  273. # geometry type
  274. gtype = libvect.GV_POINT
  275. def __init__(self, x=0, y=0, z=None, **kargs):
  276. super(Point, self).__init__(**kargs)
  277. if self.id is not None:
  278. self.read()
  279. else:
  280. self.is2D = True if z is None else False
  281. z = z if z is not None else 0
  282. libvect.Vect_append_point(self.c_points, x, y, z)
  283. def _get_x(self):
  284. return self.c_points.contents.x[0]
  285. def _set_x(self, value):
  286. self.c_points.contents.x[0] = value
  287. x = property(fget=_get_x, fset=_set_x,
  288. doc="Set and obtain x coordinate")
  289. def _get_y(self):
  290. return self.c_points.contents.y[0]
  291. def _set_y(self, value):
  292. self.c_points.contents.y[0] = value
  293. y = property(fget=_get_y, fset=_set_y,
  294. doc="Set and obtain y coordinate")
  295. def _get_z(self):
  296. if self.is2D:
  297. return None
  298. return self.c_points.contents.z[0]
  299. def _set_z(self, value):
  300. if value is None:
  301. self.is2D = True
  302. self.c_points.contents.z[0] = 0
  303. else:
  304. self.c_points.contents.z[0] = value
  305. self.is2D = False
  306. z = property(fget=_get_z, fset=_set_z,
  307. doc="Set and obtain z coordinate")
  308. def __str__(self):
  309. return self.get_wkt()
  310. def __repr__(self):
  311. return "Point(%s)" % ', '.join(['%f' % coor for coor in self.coords()])
  312. def __eq__(self, pnt):
  313. """Return True if the coordinates are the same.
  314. >>> p0 = Point()
  315. >>> p1 = Point()
  316. >>> p2 = Point(1, 1)
  317. >>> p0 == p1
  318. True
  319. >>> p1 == p2
  320. False
  321. """
  322. if isinstance(pnt, Point):
  323. return pnt.coords() == self.coords()
  324. return Point(*pnt).coords() == self.coords()
  325. def __ne__(self, other):
  326. return not self == other
  327. # Restore Python 2 hashing beaviour on Python 3
  328. __hash__ = object.__hash__
  329. def coords(self):
  330. """Return a tuple with the point coordinates. ::
  331. >>> pnt = Point(10, 100)
  332. >>> pnt.coords()
  333. (10.0, 100.0)
  334. If the point is 2D return a x, y tuple. But if we change the ``z``
  335. the Point object become a 3D point, therefore the method return a
  336. x, y, z tuple. ::
  337. >>> pnt.z = 1000.
  338. >>> pnt.coords()
  339. (10.0, 100.0, 1000.0)
  340. ..
  341. """
  342. if self.is2D:
  343. return self.x, self.y
  344. else:
  345. return self.x, self.y, self.z
  346. def get_wkt(self):
  347. """Return a "well know text" (WKT) geometry string. ::
  348. >>> pnt = Point(10, 100)
  349. >>> pnt.get_wkt()
  350. 'POINT(10.000000 100.000000)'
  351. .. warning::
  352. Only ``POINT`` (2/3D) are supported, ``POINTM`` and ``POINT`` with:
  353. ``XYZM`` are not supported yet.
  354. """
  355. return "POINT(%s)" % ' '.join(['%f' % coord
  356. for coord in self.coords()])
  357. def get_wkb(self):
  358. """Return a "well know binary" (WKB) geometry buffer
  359. .. warning::
  360. Not implemented yet.
  361. """
  362. pass
  363. def distance(self, pnt):
  364. """Calculate distance of 2 points, using the Vect_points_distance
  365. C function, If one of the point have z == None, return the 2D distance.
  366. :param pnt: the point for calculate the distance
  367. :type pnt: a Point object or a tuple with the coordinates
  368. >>> pnt0 = Point(0, 0, 0)
  369. >>> pnt1 = Point(1, 0)
  370. >>> pnt0.distance(pnt1)
  371. 1.0
  372. >>> pnt1.z = 1
  373. >>> pnt1
  374. Point(1.000000, 0.000000, 1.000000)
  375. >>> pnt0.distance(pnt1)
  376. 1.4142135623730951
  377. """
  378. if self.is2D or pnt.is2D:
  379. return libvect.Vect_points_distance(self.x, self.y, 0,
  380. pnt.x, pnt.y, 0, 0)
  381. else:
  382. return libvect.Vect_points_distance(self.x, self.y, self.z,
  383. pnt.x, pnt.y, pnt.z, 1)
  384. def buffer(self, dist=None, dist_x=None, dist_y=None, angle=0,
  385. round_=True, tol=0.1):
  386. """Return the buffer area around the point, using the
  387. ``Vect_point_buffer2`` C function.
  388. :param dist: the distance around the point
  389. :type dist: num
  390. :param dist_x: the distance along x
  391. :type dist_x: num
  392. :param dist_y: the distance along y
  393. :type dist_y: num
  394. :param angle: the angle between 0x and major axis
  395. :type angle: num
  396. :param round_: to make corners round
  397. :type round_: bool
  398. :param tol: fix the maximum distance between theoretical arc and
  399. output segments
  400. :type tol: float
  401. :returns: the buffer as Area object
  402. >>> pnt = Point(0, 0)
  403. >>> area = pnt.buffer(10)
  404. >>> area.boundary #doctest: +ELLIPSIS
  405. Line([Point(10.000000, 0.000000),...Point(10.000000, 0.000000)])
  406. >>> area.centroid
  407. Point(0.000000, 0.000000)
  408. >>> area.isles
  409. []
  410. """
  411. if dist is not None:
  412. dist_x = dist
  413. dist_y = dist
  414. elif not dist_x or not dist_y:
  415. raise TypeError('TypeError: buffer expected 1 arguments, got 0')
  416. bound = Line()
  417. p_points = ctypes.pointer(bound.c_points)
  418. libvect.Vect_point_buffer2(self.x, self.y,
  419. dist_x, dist_y,
  420. angle, int(round_), tol,
  421. p_points)
  422. return Area(boundary=bound, centroid=self)
  423. class Line(Geo):
  424. """Instantiate a new Line with a list of tuple, or with a list of Point. ::
  425. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  426. >>> line #doctest: +NORMALIZE_WHITESPACE
  427. Line([Point(0.000000, 0.000000),
  428. Point(1.000000, 1.000000),
  429. Point(2.000000, 0.000000),
  430. Point(1.000000, -1.000000)])
  431. ..
  432. """
  433. # geometry type
  434. gtype = libvect.GV_LINE
  435. def __init__(self, points=None, **kargs):
  436. super(Line, self).__init__(**kargs)
  437. if points is not None:
  438. for pnt in points:
  439. self.append(pnt)
  440. def __getitem__(self, key):
  441. """Get line point of given index, slice allowed. ::
  442. >>> line = Line([(0, 0), (1, 1), (2, 2), (3, 3)])
  443. >>> line[1]
  444. Point(1.000000, 1.000000)
  445. >>> line[-1]
  446. Point(3.000000, 3.000000)
  447. >>> line[:2]
  448. [Point(0.000000, 0.000000), Point(1.000000, 1.000000)]
  449. ..
  450. """
  451. #TODO:
  452. # line[0].x = 10 is not working
  453. #pnt.c_px = ctypes.pointer(self.c_points.contents.x[indx])
  454. # pnt.c_px = ctypes.cast(id(self.c_points.contents.x[indx]),
  455. # ctypes.POINTER(ctypes.c_double))
  456. if isinstance(key, slice):
  457. #import pdb; pdb.set_trace()
  458. #Get the start, stop, and step from the slice
  459. return [Point(self.c_points.contents.x[indx],
  460. self.c_points.contents.y[indx],
  461. None if self.is2D else self.c_points.contents.z[indx])
  462. for indx in range(*key.indices(len(self)))]
  463. elif isinstance(key, int):
  464. if key < 0: # Handle negative indices
  465. key += self.c_points.contents.n_points
  466. if key >= self.c_points.contents.n_points:
  467. raise IndexError('Index out of range')
  468. return Point(self.c_points.contents.x[key],
  469. self.c_points.contents.y[key],
  470. None if self.is2D else self.c_points.contents.z[key])
  471. else:
  472. raise ValueError("Invalid argument type: %r." % key)
  473. def __setitem__(self, indx, pnt):
  474. """Change the coordinate of point. ::
  475. >>> line = Line([(0, 0), (1, 1)])
  476. >>> line[0] = (2, 2)
  477. >>> line
  478. Line([Point(2.000000, 2.000000), Point(1.000000, 1.000000)])
  479. ..
  480. """
  481. x, y, z = get_xyz(pnt)
  482. self.c_points.contents.x[indx] = x
  483. self.c_points.contents.y[indx] = y
  484. self.c_points.contents.z[indx] = z
  485. def __iter__(self):
  486. """Return a Point generator of the Line"""
  487. return (self.__getitem__(i) for i in range(self.__len__()))
  488. def __len__(self):
  489. """Return the number of points of the line."""
  490. return self.c_points.contents.n_points
  491. def __str__(self):
  492. return self.get_wkt()
  493. def __repr__(self):
  494. return "Line([%s])" % ', '.join([repr(pnt) for pnt in self.__iter__()])
  495. def get_pnt(self, distance, angle=0, slope=0):
  496. """Return a Point object on line in the specified distance, using the
  497. `Vect_point_on_line` C function.
  498. Raise a ValueError If the distance exceed the Line length. ::
  499. >>> line = Line([(0, 0), (1, 1)])
  500. >>> line.get_pnt(5) #doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
  501. Traceback (most recent call last):
  502. ...
  503. ValueError: The distance exceed the lenght of the line,
  504. that is: 1.414214
  505. >>> line.get_pnt(1)
  506. Point(0.707107, 0.707107)
  507. ..
  508. """
  509. # instantiate an empty Point object
  510. maxdist = self.length()
  511. if distance > maxdist:
  512. str_err = "The distance exceed the lenght of the line, that is: %f"
  513. raise ValueError(str_err % maxdist)
  514. pnt = Point(0, 0, -9999)
  515. libvect.Vect_point_on_line(self.c_points, distance,
  516. pnt.c_points.contents.x,
  517. pnt.c_points.contents.y,
  518. pnt.c_points.contents.z,
  519. angle, slope)
  520. pnt.is2D = self.is2D
  521. return pnt
  522. def append(self, pnt):
  523. """Appends one point to the end of a line, using the
  524. ``Vect_append_point`` C function.
  525. :param pnt: the point to add to line
  526. :type pnt: a Point object or a tuple with the coordinates
  527. >>> line = Line()
  528. >>> line.append((10, 100))
  529. >>> line
  530. Line([Point(10.000000, 100.000000)])
  531. >>> line.append((20, 200))
  532. >>> line
  533. Line([Point(10.000000, 100.000000), Point(20.000000, 200.000000)])
  534. Like python list.
  535. """
  536. x, y, z = get_xyz(pnt)
  537. libvect.Vect_append_point(self.c_points, x, y, z)
  538. def bbox(self, bbox=None):
  539. """Return the bounding box of the line, using ``Vect_line_box``
  540. C function. ::
  541. >>> line = Line([(0, 0), (0, 1), (2, 1), (2, 0)])
  542. >>> bbox = line.bbox()
  543. >>> bbox
  544. Bbox(1.0, 0.0, 2.0, 0.0)
  545. ..
  546. """
  547. bbox = bbox if bbox else Bbox()
  548. libvect.Vect_line_box(self.c_points, bbox.c_bbox)
  549. return bbox
  550. def extend(self, line, forward=True):
  551. """Appends points to the end of a line.
  552. :param line: it is possible to extend a line, give a list of points,
  553. or directly with a line_pnts struct.
  554. :type line: Line object ot list of points
  555. :param forward: if forward is True the line is extend forward otherwise
  556. is extend backward. The method use the
  557. `Vect_append_points` C function.
  558. :type forward: bool
  559. >>> line = Line([(0, 0), (1, 1)])
  560. >>> line.extend( Line([(2, 2), (3, 3)]) )
  561. >>> line #doctest: +NORMALIZE_WHITESPACE
  562. Line([Point(0.000000, 0.000000),
  563. Point(1.000000, 1.000000),
  564. Point(2.000000, 2.000000),
  565. Point(3.000000, 3.000000)])
  566. """
  567. # set direction
  568. if forward:
  569. direction = libvect.GV_FORWARD
  570. else:
  571. direction = libvect.GV_BACKWARD
  572. # check if is a Line object
  573. if isinstance(line, Line):
  574. c_points = line.c_points
  575. else:
  576. # instantiate a Line object
  577. lin = Line()
  578. for pnt in line:
  579. # add the points to the line
  580. lin.append(pnt)
  581. c_points = lin.c_points
  582. libvect.Vect_append_points(self.c_points, c_points, direction)
  583. def insert(self, indx, pnt):
  584. """Insert new point at index position and move all old points at
  585. that position and above up, using ``Vect_line_insert_point``
  586. C function.
  587. :param indx: the index where add new point
  588. :type indx: int
  589. :param pnt: the point to add
  590. :type pnt: a Point object
  591. >>> line = Line([(0, 0), (1, 1)])
  592. >>> line.insert(0, Point(1.000000, -1.000000) )
  593. >>> line #doctest: +NORMALIZE_WHITESPACE
  594. Line([Point(1.000000, -1.000000),
  595. Point(0.000000, 0.000000),
  596. Point(1.000000, 1.000000)])
  597. """
  598. if indx < 0: # Handle negative indices
  599. indx += self.c_points.contents.n_points
  600. if indx >= self.c_points.contents.n_points:
  601. raise IndexError('Index out of range')
  602. x, y, z = get_xyz(pnt)
  603. libvect.Vect_line_insert_point(self.c_points, indx, x, y, z)
  604. def length(self):
  605. """Calculate line length, 3D-length in case of 3D vector line, using
  606. `Vect_line_length` C function. ::
  607. >>> line = Line([(0, 0), (1, 1), (0, 1)])
  608. >>> line.length()
  609. 2.414213562373095
  610. ..
  611. """
  612. return libvect.Vect_line_length(self.c_points)
  613. def length_geodesic(self):
  614. """Calculate line length, usig `Vect_line_geodesic_length` C function.
  615. ::
  616. >>> line = Line([(0, 0), (1, 1), (0, 1)])
  617. >>> line.length_geodesic()
  618. 2.414213562373095
  619. ..
  620. """
  621. return libvect.Vect_line_geodesic_length(self.c_points)
  622. def distance(self, pnt):
  623. """Calculate the distance between line and a point.
  624. :param pnt: the point to calculate distance
  625. :type pnt: a Point object or a tuple with the coordinates
  626. Return a tuple with:
  627. * the closest point on the line,
  628. * the distance between these two points,
  629. * distance of point from segment beginning
  630. * distance of point from line
  631. The distance is compute using the ``Vect_line_distance`` C function.
  632. >>> line = Line([(0, 0), (0, 2)])
  633. >>> line.distance(Point(1, 1))
  634. (Point(0.000000, 1.000000), 1.0, 1.0, 1.0)
  635. """
  636. # instantite outputs
  637. cx = ctypes.c_double(0)
  638. cy = ctypes.c_double(0)
  639. cz = ctypes.c_double(0)
  640. dist = ctypes.c_double(0)
  641. sp_dist = ctypes.c_double(0)
  642. lp_dist = ctypes.c_double(0)
  643. libvect.Vect_line_distance(self.c_points,
  644. pnt.x, pnt.y, 0 if pnt.is2D else pnt.z,
  645. 0 if self.is2D else 1,
  646. ctypes.byref(cx), ctypes.byref(cy),
  647. ctypes.byref(cz), ctypes.byref(dist),
  648. ctypes.byref(sp_dist),
  649. ctypes.byref(lp_dist))
  650. # instantiate the Point class
  651. point = Point(cx.value, cy.value, cz.value)
  652. point.is2D = self.is2D
  653. return point, dist.value, sp_dist.value, lp_dist.value
  654. def get_first_cat(self):
  655. """Fetches FIRST category number for given vector line and field, using
  656. the ``Vect_get_line_cat`` C function.
  657. .. warning::
  658. Not implemented yet.
  659. """
  660. # TODO: add this method.
  661. libvect.Vect_get_line_cat(self.map, self.id, self.field)
  662. pass
  663. def pop(self, indx):
  664. """Return the point in the index position and remove from the Line.
  665. :param indx: the index where add new point
  666. :type indx: int
  667. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  668. >>> midle_pnt = line.pop(1)
  669. >>> midle_pnt
  670. Point(1.000000, 1.000000)
  671. >>> line
  672. Line([Point(0.000000, 0.000000), Point(2.000000, 2.000000)])
  673. """
  674. if indx < 0: # Handle negative indices
  675. indx += self.c_points.contents.n_points
  676. if indx >= self.c_points.contents.n_points:
  677. raise IndexError('Index out of range')
  678. pnt = self.__getitem__(indx)
  679. libvect.Vect_line_delete_point(self.c_points, indx)
  680. return pnt
  681. def delete(self, indx):
  682. """Remove the point in the index position.
  683. :param indx: the index where add new point
  684. :type indx: int
  685. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  686. >>> line.delete(-1)
  687. >>> line
  688. Line([Point(0.000000, 0.000000), Point(1.000000, 1.000000)])
  689. """
  690. if indx < 0: # Handle negative indices
  691. indx += self.c_points.contents.n_points
  692. if indx >= self.c_points.contents.n_points:
  693. raise IndexError('Index out of range')
  694. libvect.Vect_line_delete_point(self.c_points, indx)
  695. def prune(self):
  696. """Remove duplicate points, i.e. zero length segments, using
  697. `Vect_line_prune` C function. ::
  698. >>> line = Line([(0, 0), (1, 1), (1, 1), (2, 2)])
  699. >>> line.prune()
  700. >>> line #doctest: +NORMALIZE_WHITESPACE
  701. Line([Point(0.000000, 0.000000),
  702. Point(1.000000, 1.000000),
  703. Point(2.000000, 2.000000)])
  704. ..
  705. """
  706. libvect.Vect_line_prune(self.c_points)
  707. def prune_thresh(self, threshold):
  708. """Remove points in threshold, using the ``Vect_line_prune_thresh``
  709. C funtion.
  710. :param threshold: the threshold value where prune points
  711. :type threshold: num
  712. >>> line = Line([(0, 0), (1.0, 1.0), (1.2, 0.9), (2, 2)])
  713. >>> line.prune_thresh(0.5)
  714. >>> line #doctest: +SKIP +NORMALIZE_WHITESPACE
  715. Line([Point(0.000000, 0.000000),
  716. Point(1.000000, 1.000000),
  717. Point(2.000000, 2.000000)])
  718. .. warning ::
  719. prune_thresh is not working yet.
  720. """
  721. libvect.Vect_line_prune(self.c_points, ctypes.c_double(threshold))
  722. def remove(self, pnt):
  723. """Delete point at given index and move all points above down, using
  724. `Vect_line_delete_point` C function.
  725. :param pnt: the point to remove
  726. :type pnt: a Point object or a tuple with the coordinates
  727. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  728. >>> line.remove((2, 2))
  729. >>> line[-1]
  730. Point(1.000000, 1.000000)
  731. ..
  732. """
  733. for indx, point in enumerate(self.__iter__()):
  734. if pnt == point:
  735. libvect.Vect_line_delete_point(self.c_points, indx)
  736. return
  737. raise ValueError('list.remove(x): x not in list')
  738. def reverse(self):
  739. """Reverse the order of vertices, using `Vect_line_reverse`
  740. C function. ::
  741. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  742. >>> line.reverse()
  743. >>> line #doctest: +NORMALIZE_WHITESPACE
  744. Line([Point(2.000000, 2.000000),
  745. Point(1.000000, 1.000000),
  746. Point(0.000000, 0.000000)])
  747. ..
  748. """
  749. libvect.Vect_line_reverse(self.c_points)
  750. def segment(self, start, end):
  751. # TODO improve documentation
  752. """Create line segment. using the ``Vect_line_segment`` C function.
  753. """
  754. line = Line()
  755. libvect.Vect_line_segment(self.c_points, start, end, line.c_points)
  756. return line
  757. def tolist(self):
  758. """Return a list of tuple. ::
  759. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  760. >>> line.tolist()
  761. [(0.0, 0.0), (1.0, 1.0), (2.0, 0.0), (1.0, -1.0)]
  762. ..
  763. """
  764. return [pnt.coords() for pnt in self.__iter__()]
  765. def toarray(self):
  766. """Return an array of coordinates. ::
  767. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  768. >>> line.toarray() #doctest: +NORMALIZE_WHITESPACE
  769. array([[ 0., 0.],
  770. [ 1., 1.],
  771. [ 2., 0.],
  772. [ 1., -1.]])
  773. ..
  774. """
  775. return np.array(self.tolist())
  776. def get_wkt(self):
  777. """Return a Well Known Text string of the line. ::
  778. >>> line = Line([(0, 0), (1, 1), (1, 2)])
  779. >>> line.get_wkt() #doctest: +ELLIPSIS
  780. 'LINESTRING(0.000000 0.000000, ..., 1.000000 2.000000)'
  781. ..
  782. """
  783. return "LINESTRING(%s)" % ', '.join([
  784. ' '.join(['%f' % coord for coord in pnt.coords()])
  785. for pnt in self.__iter__()])
  786. def from_wkt(self, wkt):
  787. """Create a line reading a WKT string.
  788. :param wkt: the WKT string containing the LINESTRING
  789. :type wkt: str
  790. >>> line = Line()
  791. >>> line.from_wkt("LINESTRING(0 0,1 1,1 2)")
  792. >>> line #doctest: +NORMALIZE_WHITESPACE
  793. Line([Point(0.000000, 0.000000),
  794. Point(1.000000, 1.000000),
  795. Point(1.000000, 2.000000)])
  796. ..
  797. """
  798. match = re.match('LINESTRING\((.*)\)', wkt)
  799. if match:
  800. self.reset()
  801. for coord in match.groups()[0].strip().split(','):
  802. self.append(tuple([float(e) for e in coord.split(' ')]))
  803. else:
  804. return None
  805. def get_wkb(self):
  806. """Return a WKB buffer.
  807. .. warning::
  808. Not implemented yet.
  809. """
  810. pass
  811. def buffer(self, dist=None, dist_x=None, dist_y=None,
  812. angle=0, round_=True, caps=True, tol=0.1):
  813. """Return the buffer area around the line, using the
  814. ``Vect_line_buffer2`` C function.
  815. :param dist: the distance around the line
  816. :type dist: num
  817. :param dist_x: the distance along x
  818. :type dist_x: num
  819. :param dist_y: the distance along y
  820. :type dist_y: num
  821. :param angle: the angle between 0x and major axis
  822. :type angle: num
  823. :param round_: to make corners round
  824. :type round_: bool
  825. :param tol: fix the maximum distance between theoretical arc and
  826. output segments
  827. :type tol: float
  828. :returns: the buffer as Area object
  829. >>> line = Line([(0, 0), (0, 2)])
  830. >>> area = line.buffer(10)
  831. >>> area.boundary #doctest: +ELLIPSIS
  832. Line([Point(-10.000000, 0.000000),...Point(-10.000000, 0.000000)])
  833. >>> area.centroid
  834. Point(0.000000, 0.000000)
  835. >>> area.isles
  836. []
  837. ..
  838. """
  839. if dist is not None:
  840. dist_x = dist
  841. dist_y = dist
  842. elif not dist_x or not dist_y:
  843. raise TypeError('TypeError: buffer expected 1 arguments, got 0')
  844. p_bound = ctypes.pointer(ctypes.pointer(libvect.line_pnts()))
  845. pp_isle = ctypes.pointer(ctypes.pointer(
  846. ctypes.pointer(libvect.line_pnts())))
  847. n_isles = ctypes.pointer(ctypes.c_int())
  848. libvect.Vect_line_buffer2(self.c_points,
  849. dist_x, dist_y, angle,
  850. int(round_), int(caps), tol,
  851. p_bound, pp_isle, n_isles)
  852. return Area(boundary=Line(c_points=p_bound.contents),
  853. centroid=self[0],
  854. isles=[Line(c_points=pp_isle[i].contents)
  855. for i in range(n_isles.contents.value)])
  856. def reset(self):
  857. """Reset line, using `Vect_reset_line` C function. ::
  858. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  859. >>> len(line)
  860. 4
  861. >>> line.reset()
  862. >>> len(line)
  863. 0
  864. >>> line
  865. Line([])
  866. ..
  867. """
  868. libvect.Vect_reset_line(self.c_points)
  869. class Node(object):
  870. pass
  871. class Boundary(Line):
  872. """
  873. """
  874. # geometry type
  875. gtype = libvect.GV_BOUNDARY
  876. def __init__(self, lines=None, left=None, right=None,
  877. **kargs):
  878. v_id = kargs.get('v_id', 0)
  879. self.dir = libvect.GV_FORWARD if v_id > 0 else libvect.GV_BACKWARD
  880. super(Boundary, self).__init__(**kargs)
  881. self.c_left = ctypes.pointer(ctypes.c_int())
  882. self.c_right = ctypes.pointer(ctypes.c_int())
  883. #self.get_left_right()
  884. @property
  885. def left_id(self):
  886. return self.c_left.contents.value
  887. @property
  888. def right_id(self):
  889. return self.c_right.contents.value
  890. def __repr__(self):
  891. return "Boundary(v_id=%r)" % self.id
  892. def _get_centroid(self, side, idonly=False):
  893. if side > 0:
  894. v_id = libvect.Vect_get_area_centroid(self.c_mapinfo, side)
  895. v_id = v_id if v_id else None
  896. if idonly:
  897. return v_id
  898. else:
  899. cntr = Centroid(v_id=v_id, c_mapinfo=self.c_mapinfo)
  900. return cntr
  901. def get_left_centroid(self, idonly=False):
  902. """Return left value
  903. :param idonly: True to return only the cat of feature
  904. :type idonly: bool
  905. """
  906. return self._get_centroid(self.left_id, idonly)
  907. def get_right_centroid(self, idonly=False):
  908. """Return right value
  909. :param idonly: True to return only the cat of feature
  910. :type idonly: bool
  911. """
  912. return self._get_centroid(self.left_id, idonly)
  913. def get_left_right(self):
  914. """Return left and right value"""
  915. libvect.Vect_get_line_areas(self.c_mapinfo, self.id,
  916. self.c_left, self.c_right)
  917. return self.c_left.contents.value, self.c_right.contents.value
  918. def area(self):
  919. """Return the area of the polygon.
  920. >>> bound = Boundary(points=[(0, 0), (0, 2), (2, 2), (2, 0),
  921. ... (0, 0)])
  922. >>> bound.area()
  923. 4.0
  924. """
  925. libgis.G_begin_polygon_area_calculations()
  926. return libgis.G_area_of_polygon(self.c_points.contents.x,
  927. self.c_points.contents.y,
  928. self.c_points.contents.n_points)
  929. class Centroid(Point):
  930. """The Centroid class inherit from the Point class.
  931. Centroid contains an attribute with the C Map_info struct, and attributes
  932. with the id of the Area. ::
  933. >>> centroid = Centroid(x=0, y=10)
  934. >>> centroid
  935. Centoid(0.000000, 10.000000)
  936. >>> from grass.pygrass.vector import VectorTopo
  937. >>> geo = VectorTopo('geology')
  938. >>> geo.open(mode='r')
  939. >>> centroid = Centroid(v_id=1, c_mapinfo=geo.c_mapinfo)
  940. >>> centroid
  941. Centoid(893202.874416, 297339.312795)
  942. ..
  943. """
  944. # geometry type
  945. gtype = libvect.GV_CENTROID
  946. def __init__(self, area_id=None, **kargs):
  947. super(Centroid, self).__init__(**kargs)
  948. self.area_id = area_id
  949. if self.id and self.c_mapinfo and self.area_id is None:
  950. self.area_id = self.get_area_id()
  951. elif self.c_mapinfo and self.area_id and self.id is None:
  952. self.id = self.get_centroid_id()
  953. if self.area_id is not None:
  954. self.read()
  955. #self.c_pline = ctypes.pointer(libvect.P_line()) if topology else None
  956. def __repr__(self):
  957. return "Centoid(%s)" % ', '.join(['%f' % co for co in self.coords()])
  958. def get_centroid_id(self):
  959. """Return the centroid_id, using the c_mapinfo and an area_id
  960. attributes of the class, and calling the Vect_get_area_centroid
  961. C function, if no centroid_id were found return None"""
  962. centroid_id = libvect.Vect_get_area_centroid(self.c_mapinfo,
  963. self.area_id)
  964. return centroid_id if centroid_id != 0 else None
  965. def get_area_id(self):
  966. """Return the area_id, using the c_mapinfo and an centroid_id
  967. attributes of the class, and calling the Vect_get_centroid_area
  968. C function, if no area_id were found return None"""
  969. area_id = libvect.Vect_get_centroid_area(self.c_mapinfo,
  970. self.id)
  971. return area_id if area_id != 0 else None
  972. class Isle(Geo):
  973. """An Isle is an area contained by another area.
  974. """
  975. def __init__(self, **kargs):
  976. super(Isle, self).__init__(**kargs)
  977. #self.area_id = area_id
  978. def __repr__(self):
  979. return "Isle(%d)" % (self.id)
  980. def boundaries(self):
  981. """Return a list of boundaries"""
  982. ilist = Ilist()
  983. libvect.Vect_get_isle_boundaries(self.c_mapinfo, self.id,
  984. ilist.c_ilist)
  985. return ilist
  986. def bbox(self, bbox=None):
  987. """Return bounding box of Isle"""
  988. bbox = bbox if bbox else Bbox()
  989. libvect.Vect_get_isle_box(self.c_mapinfo, self.id, bbox.c_bbox)
  990. return bbox
  991. def points(self):
  992. """Return a Line object with the outer ring points"""
  993. line = Line()
  994. libvect.Vect_get_isle_points(self.c_mapinfo, self.id, line.c_points)
  995. return line
  996. def points_geos(self):
  997. """Return a Line object with the outer ring points
  998. """
  999. return libvect.Vect_get_isle_points_geos(self.c_mapinfo, self.id)
  1000. def area_id(self):
  1001. """Returns area id for isle."""
  1002. return libvect.Vect_get_isle_area(self.c_mapinfo, self.id)
  1003. def alive(self):
  1004. """Check if isle is alive or dead (topology required)"""
  1005. return bool(libvect.Vect_isle_alive(self.c_mapinfo, self.id))
  1006. def contain_pnt(self, pnt):
  1007. """Check if point is in area.
  1008. :param pnt: the point to remove
  1009. :type pnt: a Point object or a tuple with the coordinates
  1010. """
  1011. bbox = self.bbox()
  1012. return bool(libvect.Vect_point_in_island(pnt.x, pnt.y,
  1013. self.c_mapinfo, self.id,
  1014. bbox.c_bbox.contents))
  1015. def area(self):
  1016. """Return the area value of an Isle"""
  1017. border = self.points()
  1018. return libgis.G_area_of_polygon(border.c_points.contents.x,
  1019. border.c_points.contents.y,
  1020. border.c_points.contents.n_points)
  1021. def perimeter(self):
  1022. """Return the perimeter value of an Isle.
  1023. ::
  1024. double Vect_area_perimeter()
  1025. """
  1026. border = self.points()
  1027. return libvect.Vect_area_perimeter(border.c_points)
  1028. class Isles(object):
  1029. def __init__(self, c_mapinfo, area_id=None):
  1030. self.c_mapinfo = c_mapinfo
  1031. self.area_id = area_id
  1032. self._isles_id = None
  1033. self._isles = None
  1034. if area_id:
  1035. self._isles_id = self.get_isles_id()
  1036. self._isles = self.get_isles()
  1037. def __len__(self):
  1038. return libvect.Vect_get_area_num_isles(self.c_mapinfo, self.area_id)
  1039. def __repr__(self):
  1040. return "Isles(%r)" % self.area_id
  1041. def __getitem__(self, key):
  1042. if self._isles is None:
  1043. self.get_isles()
  1044. return self._isles[key]
  1045. def get_isles_id(self):
  1046. """Return the id of isles"""
  1047. return [libvect.Vect_get_area_isle(self.c_mapinfo, self.area_id, i)
  1048. for i in range(self.__len__())]
  1049. def get_isles(self):
  1050. """Return isles"""
  1051. return [Isle(v_id=isle_id, c_mapinfo=self.c_mapinfo)
  1052. for isle_id in self._isles_id]
  1053. def select_by_bbox(self, bbox):
  1054. """Vect_select_isles_by_box
  1055. .. warning::
  1056. Not implemented yet.
  1057. """
  1058. pass
  1059. class Area(Geo):
  1060. """
  1061. 'Vect_build_line_area',
  1062. 'Vect_find_area',
  1063. 'Vect_get_area_box',
  1064. 'Vect_get_area_points_geos',
  1065. 'Vect_get_centroid_area',
  1066. 'Vect_get_isle_area',
  1067. 'Vect_get_line_areas',
  1068. 'Vect_get_num_areas',
  1069. 'Vect_get_point_in_area',
  1070. 'Vect_isle_find_area',
  1071. 'Vect_point_in_area',
  1072. 'Vect_point_in_area_outer_ring',
  1073. 'Vect_read_area_geos',
  1074. 'Vect_remove_small_areas',
  1075. 'Vect_select_areas_by_box',
  1076. 'Vect_select_areas_by_polygon']
  1077. """
  1078. # geometry type
  1079. gtype = libvect.GV_AREA
  1080. def __init__(self, boundary=None, centroid=None, isles=None, **kargs):
  1081. super(Area, self).__init__(**kargs)
  1082. self.boundary = None
  1083. self.centroid = None
  1084. self.isles = None
  1085. if boundary and centroid:
  1086. self.boundary = boundary
  1087. self.centroid = centroid
  1088. self.isles = isles if isles else []
  1089. # set the attributes
  1090. if self.attrs and self.cat:
  1091. self.attrs.cat = self.cat
  1092. def __repr__(self):
  1093. return "Area(%d)" % self.id if self.id else "Area( )"
  1094. def init_from_id(self, area_id=None):
  1095. """Return an Area object"""
  1096. if area_id is None and self.id is None:
  1097. raise ValueError("You need to give or set the area_id")
  1098. self.id = area_id if area_id is not None else self.id
  1099. # get boundary
  1100. self.get_points()
  1101. # get isles
  1102. self.get_isles()
  1103. def get_points(self, line=None):
  1104. """Return a Line object with the outer ring
  1105. :param line: a Line object to fill with info from points of area
  1106. :type line: a Line object
  1107. """
  1108. line = Line() if line is None else line
  1109. libvect.Vect_get_area_points(self.c_mapinfo, self.id, line.c_points)
  1110. return line
  1111. def get_centroid(self, centroid=None):
  1112. """Return the centroid
  1113. :param centroid: a Centroid object to fill with info from centroid of area
  1114. :type centroid: a Centroid object
  1115. """
  1116. centroid_id = libvect.Vect_get_area_centroid(self.c_mapinfo, self.id)
  1117. if centroid:
  1118. centroid.id = centroid_id
  1119. centroid.read()
  1120. return centroid
  1121. return Centroid(v_id=centroid_id, c_mapinfo=self.c_mapinfo,
  1122. area_id=self.id)
  1123. def num_isles(self):
  1124. return libvect.Vect_get_area_num_isles(self.c_mapinfo, self.id)
  1125. def get_isles(self, isles=None):
  1126. """Instantiate the boundary attribute reading area_id"""
  1127. if isles is not None:
  1128. isles.area_id = self.id
  1129. return isles
  1130. return Isles(self.c_mapinfo, self.id)
  1131. def area(self):
  1132. """Returns area of area without areas of isles.
  1133. double Vect_get_area_area (const struct Map_info *Map, int area)
  1134. """
  1135. return libvect.Vect_get_area_area(self.c_mapinfo, self.id)
  1136. def alive(self):
  1137. """Check if area is alive or dead (topology required)
  1138. """
  1139. return bool(libvect.Vect_area_alive(self.c_mapinfo, self.id))
  1140. def bbox(self, bbox=None):
  1141. """Return the Bbox of area
  1142. :param bbox: a Bbox object to fill with info from bounding box of area
  1143. :type bbox: a Bbox object
  1144. """
  1145. bbox = bbox if bbox else Bbox()
  1146. libvect.Vect_get_area_box(self.c_mapinfo, self.id, bbox.c_bbox)
  1147. return bbox
  1148. def buffer(self, dist=None, dist_x=None, dist_y=None,
  1149. angle=0, round_=True, caps=True, tol=0.1):
  1150. """Return the buffer area around the area, using the
  1151. ``Vect_area_buffer2`` C function.
  1152. :param dist: the distance around the area
  1153. :type dist: num
  1154. :param dist_x: the distance along x
  1155. :type dist_x: num
  1156. :param dist_y: the distance along y
  1157. :type dist_y: num
  1158. :param angle: the angle between 0x and major axis
  1159. :type angle: num
  1160. :param round_: to make corners round
  1161. :type round_: bool
  1162. :param tol: fix the maximum distance between theoretical arc and
  1163. output segments
  1164. :type tol: float
  1165. :returns: the buffer as Area object
  1166. """
  1167. if dist is not None:
  1168. dist_x = dist
  1169. dist_y = dist
  1170. elif not dist_x or not dist_y:
  1171. raise TypeError('TypeError: buffer expected 1 arguments, got 0')
  1172. p_bound = ctypes.pointer(ctypes.pointer(libvect.line_pnts()))
  1173. pp_isle = ctypes.pointer(ctypes.pointer(
  1174. ctypes.pointer(libvect.line_pnts())))
  1175. n_isles = ctypes.pointer(ctypes.c_int())
  1176. libvect.Vect_area_buffer2(self.c_mapinfo, self.id,
  1177. dist_x, dist_y, angle,
  1178. int(round_), int(caps), tol,
  1179. p_bound, pp_isle, n_isles)
  1180. return Area(boundary=Line(c_points=p_bound.contents),
  1181. centroid=self.centroid,
  1182. isles=[Line(c_points=pp_isle[i].contents)
  1183. for i in range(n_isles.contents.value)])
  1184. def boundaries(self, ilist=False):
  1185. """Creates list of boundaries for given area.
  1186. int Vect_get_area_boundaries(const struct Map_info *Map,
  1187. int area, struct ilist *List)
  1188. """
  1189. ilst = Ilist()
  1190. libvect.Vect_get_area_boundaries(self.c_mapinfo, self.id,
  1191. ilst.c_ilist)
  1192. if ilist:
  1193. return ilist
  1194. return [Boundary(v_id, c_mapinfo=self.c_mapinfo) for v_id in ilst]
  1195. def cats(self, cats=None):
  1196. """Get area categories.
  1197. :param cats: a Cats object to fill with info with area categories
  1198. :type cats: a Cats object
  1199. """
  1200. cats = cats if cats else Cats()
  1201. libvect.Vect_get_area_cats(self.c_mapinfo, self.id, cats.c_cats)
  1202. return cats
  1203. def get_first_cat(self):
  1204. """Find FIRST category of given field and area.
  1205. int Vect_get_area_cat(const struct Map_info *Map, int area, int field)
  1206. """
  1207. pass
  1208. def contain_pnt(self, pnt, bbox=None):
  1209. """Check if point is in area.
  1210. :param pnt: the point to analyze
  1211. :type pnt: a Point object or a tuple with the coordinates
  1212. :param bbox: the bounding box where run the analysis
  1213. :type bbox: a Bbox object
  1214. """
  1215. bbox = bbox if bbox else self.bbox()
  1216. return bool(libvect.Vect_point_in_area(pnt.x, pnt.y,
  1217. self.c_mapinfo, self.id,
  1218. bbox.c_bbox))
  1219. def perimeter(self):
  1220. """Calculate area perimeter.
  1221. double Vect_area_perimeter (const struct line_pnts *Points)
  1222. """
  1223. border = self.get_points()
  1224. return libvect.Vect_area_perimeter(border.c_points)
  1225. def read(self, line=None, centroid=None, isles=None):
  1226. self.boundary = self.get_points(line)
  1227. self.centroid = self.get_centroid(centroid)
  1228. #self.isles = self.get_isles(isles)
  1229. libvect.Vect_read_line(self.c_mapinfo, None, self.c_cats,
  1230. self.centroid.id)
  1231. #
  1232. # Define a dictionary to convert the feature type to name and or object
  1233. #
  1234. GV_TYPE = {libvect.GV_POINT: {'label': 'point', 'obj': Point},
  1235. libvect.GV_LINE: {'label': 'line', 'obj': Line},
  1236. libvect.GV_BOUNDARY: {'label': 'boundary', 'obj': Boundary},
  1237. libvect.GV_CENTROID: {'label': 'centroid', 'obj': Centroid},
  1238. libvect.GV_FACE: {'label': 'face', 'obj': None},
  1239. libvect.GV_KERNEL: {'label': 'kernel', 'obj': None},
  1240. libvect.GV_AREA: {'label': 'area', 'obj': Area},
  1241. libvect.GV_VOLUME: {'label': 'volume', 'obj': None}, }
  1242. GEOOBJ = {"areas": Area,
  1243. "dblinks": None,
  1244. "faces": None,
  1245. "holes": None,
  1246. "islands": Isle,
  1247. "kernels": None,
  1248. "line_points": None,
  1249. "points": Point,
  1250. "lines": Line,
  1251. "nodes": Node,
  1252. "volumes": None}
  1253. def c_read_next_line(c_mapinfo, c_points, c_cats):
  1254. v_id = c_mapinfo.contents.next_line
  1255. v_id = v_id if v_id != 0 else None
  1256. ftype = libvect.Vect_read_next_line(c_mapinfo, c_points, c_cats)
  1257. if ftype == -2:
  1258. raise StopIteration()
  1259. if ftype == -1:
  1260. raise
  1261. return ftype, v_id, c_points, c_cats
  1262. def read_next_line(c_mapinfo, table=None, writable=False,
  1263. c_points=None, c_cats=None, is2D=True):
  1264. """Return the next geometry feature of a vector map."""
  1265. c_points = c_points if c_points else ctypes.pointer(libvect.line_pnts())
  1266. c_cats = c_cats if c_cats else ctypes.pointer(libvect.line_cats())
  1267. ftype, v_id, c_points, c_cats = c_read_next_line(c_mapinfo, c_points,
  1268. c_cats)
  1269. return GV_TYPE[ftype]['obj'](v_id=v_id, c_mapinfo=c_mapinfo,
  1270. c_points=c_points, c_cats=c_cats,
  1271. table=table, writable=writable, is2D=is2D)
  1272. def c_read_line(feature_id, c_mapinfo, c_points, c_cats):
  1273. nmax = libvect.Vect_get_num_lines(c_mapinfo)
  1274. if feature_id < 0: # Handle negative indices
  1275. feature_id += nmax + 1
  1276. if feature_id > nmax:
  1277. raise IndexError('Index out of range')
  1278. if feature_id > 0:
  1279. ftype = libvect.Vect_read_line(c_mapinfo, c_points, c_cats, feature_id)
  1280. return ftype, c_points, c_cats
  1281. else:
  1282. raise ValueError('The index must be >0, %r given.' % feature_id)
  1283. def read_line(feature_id, c_mapinfo, table=None, writable=False,
  1284. c_points=None, c_cats=None, is2D=True):
  1285. """Return a geometry object given the feature id and the c_mapinfo.
  1286. """
  1287. c_points = c_points if c_points else ctypes.pointer(libvect.line_pnts())
  1288. c_cats = c_cats if c_cats else ctypes.pointer(libvect.line_cats())
  1289. ftype, c_points, c_cats = c_read_line(feature_id, c_mapinfo,
  1290. c_points, c_cats)
  1291. if GV_TYPE[ftype]['obj'] is not None:
  1292. return GV_TYPE[ftype]['obj'](v_id=feature_id, c_mapinfo=c_mapinfo,
  1293. c_points=c_points, c_cats=c_cats,
  1294. table=table, writable=writable, is2D=is2D)