geometry.py 51 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570
  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. if isinstance(pnt, Point):
  314. return pnt.coords() == self.coords()
  315. return Point(*pnt).coords() == self.coords()
  316. def coords(self):
  317. """Return a tuple with the point coordinates. ::
  318. >>> pnt = Point(10, 100)
  319. >>> pnt.coords()
  320. (10.0, 100.0)
  321. If the point is 2D return a x, y tuple. But if we change the ``z``
  322. the Point object become a 3D point, therefore the method return a
  323. x, y, z tuple. ::
  324. >>> pnt.z = 1000.
  325. >>> pnt.coords()
  326. (10.0, 100.0, 1000.0)
  327. ..
  328. """
  329. if self.is2D:
  330. return self.x, self.y
  331. else:
  332. return self.x, self.y, self.z
  333. def get_wkt(self):
  334. """Return a "well know text" (WKT) geometry string. ::
  335. >>> pnt = Point(10, 100)
  336. >>> pnt.get_wkt()
  337. 'POINT(10.000000 100.000000)'
  338. .. warning::
  339. Only ``POINT`` (2/3D) are supported, ``POINTM`` and ``POINT`` with:
  340. ``XYZM`` are not supported yet.
  341. """
  342. return "POINT(%s)" % ' '.join(['%f' % coord
  343. for coord in self.coords()])
  344. def get_wkb(self):
  345. """Return a "well know binary" (WKB) geometry buffer
  346. .. warning::
  347. Not implemented yet.
  348. """
  349. pass
  350. def distance(self, pnt):
  351. """Calculate distance of 2 points, using the Vect_points_distance
  352. C function, If one of the point have z == None, return the 2D distance.
  353. :param pnt: the point for calculate the distance
  354. :type pnt: a Point object or a tuple with the coordinates
  355. >>> pnt0 = Point(0, 0, 0)
  356. >>> pnt1 = Point(1, 0)
  357. >>> pnt0.distance(pnt1)
  358. 1.0
  359. >>> pnt1.z = 1
  360. >>> pnt1
  361. Point(1.000000, 0.000000, 1.000000)
  362. >>> pnt0.distance(pnt1)
  363. 1.4142135623730951
  364. """
  365. if self.is2D or pnt.is2D:
  366. return libvect.Vect_points_distance(self.x, self.y, 0,
  367. pnt.x, pnt.y, 0, 0)
  368. else:
  369. return libvect.Vect_points_distance(self.x, self.y, self.z,
  370. pnt.x, pnt.y, pnt.z, 1)
  371. def buffer(self, dist=None, dist_x=None, dist_y=None, angle=0,
  372. round_=True, tol=0.1):
  373. """Return the buffer area around the point, using the
  374. ``Vect_point_buffer2`` C function.
  375. :param dist: the distance around the point
  376. :type dist: num
  377. :param dist_x: the distance along x
  378. :type dist_x: num
  379. :param dist_y: the distance along y
  380. :type dist_y: num
  381. :param angle: the angle between 0x and major axis
  382. :type angle: num
  383. :param round_: to make corners round
  384. :type round_: bool
  385. :param tol: fix the maximum distance between theoretical arc and
  386. output segments
  387. :type tol: float
  388. :returns: the buffer as Area object
  389. >>> pnt = Point(0, 0)
  390. >>> area = pnt.buffer(10)
  391. >>> area.boundary #doctest: +ELLIPSIS
  392. Line([Point(10.000000, 0.000000),...Point(10.000000, 0.000000)])
  393. >>> area.centroid
  394. Point(0.000000, 0.000000)
  395. >>> area.isles
  396. []
  397. """
  398. if dist is not None:
  399. dist_x = dist
  400. dist_y = dist
  401. elif not dist_x or not dist_y:
  402. raise TypeError('TypeError: buffer expected 1 arguments, got 0')
  403. bound = Line()
  404. p_points = ctypes.pointer(bound.c_points)
  405. libvect.Vect_point_buffer2(self.x, self.y,
  406. dist_x, dist_y,
  407. angle, int(round_), tol,
  408. p_points)
  409. return Area(boundary=bound, centroid=self)
  410. class Line(Geo):
  411. """Instantiate a new Line with a list of tuple, or with a list of Point. ::
  412. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  413. >>> line #doctest: +NORMALIZE_WHITESPACE
  414. Line([Point(0.000000, 0.000000),
  415. Point(1.000000, 1.000000),
  416. Point(2.000000, 0.000000),
  417. Point(1.000000, -1.000000)])
  418. ..
  419. """
  420. # geometry type
  421. gtype = libvect.GV_LINE
  422. def __init__(self, points=None, **kargs):
  423. super(Line, self).__init__(**kargs)
  424. if points is not None:
  425. for pnt in points:
  426. self.append(pnt)
  427. def __getitem__(self, key):
  428. """Get line point of given index, slice allowed. ::
  429. >>> line = Line([(0, 0), (1, 1), (2, 2), (3, 3)])
  430. >>> line[1]
  431. Point(1.000000, 1.000000)
  432. >>> line[-1]
  433. Point(3.000000, 3.000000)
  434. >>> line[:2]
  435. [Point(0.000000, 0.000000), Point(1.000000, 1.000000)]
  436. ..
  437. """
  438. #TODO:
  439. # line[0].x = 10 is not working
  440. #pnt.c_px = ctypes.pointer(self.c_points.contents.x[indx])
  441. # pnt.c_px = ctypes.cast(id(self.c_points.contents.x[indx]),
  442. # ctypes.POINTER(ctypes.c_double))
  443. if isinstance(key, slice):
  444. #import pdb; pdb.set_trace()
  445. #Get the start, stop, and step from the slice
  446. return [Point(self.c_points.contents.x[indx],
  447. self.c_points.contents.y[indx],
  448. None if self.is2D else self.c_points.contents.z[indx])
  449. for indx in range(*key.indices(len(self)))]
  450. elif isinstance(key, int):
  451. if key < 0: # Handle negative indices
  452. key += self.c_points.contents.n_points
  453. if key >= self.c_points.contents.n_points:
  454. raise IndexError('Index out of range')
  455. return Point(self.c_points.contents.x[key],
  456. self.c_points.contents.y[key],
  457. None if self.is2D else self.c_points.contents.z[key])
  458. else:
  459. raise ValueError("Invalid argument type: %r." % key)
  460. def __setitem__(self, indx, pnt):
  461. """Change the coordinate of point. ::
  462. >>> line = Line([(0, 0), (1, 1)])
  463. >>> line[0] = (2, 2)
  464. >>> line
  465. Line([Point(2.000000, 2.000000), Point(1.000000, 1.000000)])
  466. ..
  467. """
  468. x, y, z = get_xyz(pnt)
  469. self.c_points.contents.x[indx] = x
  470. self.c_points.contents.y[indx] = y
  471. self.c_points.contents.z[indx] = z
  472. def __iter__(self):
  473. """Return a Point generator of the Line"""
  474. return (self.__getitem__(i) for i in range(self.__len__()))
  475. def __len__(self):
  476. """Return the number of points of the line."""
  477. return self.c_points.contents.n_points
  478. def __str__(self):
  479. return self.get_wkt()
  480. def __repr__(self):
  481. return "Line([%s])" % ', '.join([repr(pnt) for pnt in self.__iter__()])
  482. def get_pnt(self, distance, angle=0, slope=0):
  483. """Return a Point object on line in the specified distance, using the
  484. `Vect_point_on_line` C function.
  485. Raise a ValueError If the distance exceed the Line length. ::
  486. >>> line = Line([(0, 0), (1, 1)])
  487. >>> line.get_pnt(5) #doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
  488. Traceback (most recent call last):
  489. ...
  490. ValueError: The distance exceed the lenght of the line,
  491. that is: 1.414214
  492. >>> line.get_pnt(1)
  493. Point(0.707107, 0.707107)
  494. ..
  495. """
  496. # instantiate an empty Point object
  497. maxdist = self.length()
  498. if distance > maxdist:
  499. str_err = "The distance exceed the lenght of the line, that is: %f"
  500. raise ValueError(str_err % maxdist)
  501. pnt = Point(0, 0, -9999)
  502. libvect.Vect_point_on_line(self.c_points, distance,
  503. pnt.c_points.contents.x,
  504. pnt.c_points.contents.y,
  505. pnt.c_points.contents.z,
  506. angle, slope)
  507. pnt.is2D = self.is2D
  508. return pnt
  509. def append(self, pnt):
  510. """Appends one point to the end of a line, using the
  511. ``Vect_append_point`` C function.
  512. :param pnt: the point to add to line
  513. :type pnt: a Point object or a tuple with the coordinates
  514. >>> line = Line()
  515. >>> line.append((10, 100))
  516. >>> line
  517. Line([Point(10.000000, 100.000000)])
  518. >>> line.append((20, 200))
  519. >>> line
  520. Line([Point(10.000000, 100.000000), Point(20.000000, 200.000000)])
  521. Like python list.
  522. """
  523. x, y, z = get_xyz(pnt)
  524. libvect.Vect_append_point(self.c_points, x, y, z)
  525. def bbox(self, bbox=None):
  526. """Return the bounding box of the line, using ``Vect_line_box``
  527. C function. ::
  528. >>> line = Line([(0, 0), (0, 1), (2, 1), (2, 0)])
  529. >>> bbox = line.bbox()
  530. >>> bbox
  531. Bbox(1.0, 0.0, 2.0, 0.0)
  532. ..
  533. """
  534. bbox = bbox if bbox else Bbox()
  535. libvect.Vect_line_box(self.c_points, bbox.c_bbox)
  536. return bbox
  537. def extend(self, line, forward=True):
  538. """Appends points to the end of a line.
  539. :param line: it is possible to extend a line, give a list of points,
  540. or directly with a line_pnts struct.
  541. :type line: Line object ot list of points
  542. :param forward: if forward is True the line is extend forward otherwise
  543. is extend backward. The method use the
  544. `Vect_append_points` C function.
  545. :type forward: bool
  546. >>> line = Line([(0, 0), (1, 1)])
  547. >>> line.extend( Line([(2, 2), (3, 3)]) )
  548. >>> line #doctest: +NORMALIZE_WHITESPACE
  549. Line([Point(0.000000, 0.000000),
  550. Point(1.000000, 1.000000),
  551. Point(2.000000, 2.000000),
  552. Point(3.000000, 3.000000)])
  553. """
  554. # set direction
  555. if forward:
  556. direction = libvect.GV_FORWARD
  557. else:
  558. direction = libvect.GV_BACKWARD
  559. # check if is a Line object
  560. if isinstance(line, Line):
  561. c_points = line.c_points
  562. else:
  563. # instantiate a Line object
  564. lin = Line()
  565. for pnt in line:
  566. # add the points to the line
  567. lin.append(pnt)
  568. c_points = lin.c_points
  569. libvect.Vect_append_points(self.c_points, c_points, direction)
  570. def insert(self, indx, pnt):
  571. """Insert new point at index position and move all old points at
  572. that position and above up, using ``Vect_line_insert_point``
  573. C function.
  574. :param indx: the index where add new point
  575. :type indx: int
  576. :param pnt: the point to add
  577. :type pnt: a Point object
  578. >>> line = Line([(0, 0), (1, 1)])
  579. >>> line.insert(0, Point(1.000000, -1.000000) )
  580. >>> line #doctest: +NORMALIZE_WHITESPACE
  581. Line([Point(1.000000, -1.000000),
  582. Point(0.000000, 0.000000),
  583. Point(1.000000, 1.000000)])
  584. """
  585. if indx < 0: # Handle negative indices
  586. indx += self.c_points.contents.n_points
  587. if indx >= self.c_points.contents.n_points:
  588. raise IndexError('Index out of range')
  589. x, y, z = get_xyz(pnt)
  590. libvect.Vect_line_insert_point(self.c_points, indx, x, y, z)
  591. def length(self):
  592. """Calculate line length, 3D-length in case of 3D vector line, using
  593. `Vect_line_length` C function. ::
  594. >>> line = Line([(0, 0), (1, 1), (0, 1)])
  595. >>> line.length()
  596. 2.414213562373095
  597. ..
  598. """
  599. return libvect.Vect_line_length(self.c_points)
  600. def length_geodesic(self):
  601. """Calculate line length, usig `Vect_line_geodesic_length` C function.
  602. ::
  603. >>> line = Line([(0, 0), (1, 1), (0, 1)])
  604. >>> line.length_geodesic()
  605. 2.414213562373095
  606. ..
  607. """
  608. return libvect.Vect_line_geodesic_length(self.c_points)
  609. def distance(self, pnt):
  610. """Calculate the distance between line and a point.
  611. :param pnt: the point to calculate distance
  612. :type pnt: a Point object or a tuple with the coordinates
  613. Return a tuple with:
  614. * the closest point on the line,
  615. * the distance between these two points,
  616. * distance of point from segment beginning
  617. * distance of point from line
  618. The distance is compute using the ``Vect_line_distance`` C function.
  619. >>> line = Line([(0, 0), (0, 2)])
  620. >>> line.distance(Point(1, 1))
  621. (Point(0.000000, 1.000000), 1.0, 1.0, 1.0)
  622. """
  623. # instantite outputs
  624. cx = ctypes.c_double(0)
  625. cy = ctypes.c_double(0)
  626. cz = ctypes.c_double(0)
  627. dist = ctypes.c_double(0)
  628. sp_dist = ctypes.c_double(0)
  629. lp_dist = ctypes.c_double(0)
  630. libvect.Vect_line_distance(self.c_points,
  631. pnt.x, pnt.y, 0 if pnt.is2D else pnt.z,
  632. 0 if self.is2D else 1,
  633. ctypes.byref(cx), ctypes.byref(cy),
  634. ctypes.byref(cz), ctypes.byref(dist),
  635. ctypes.byref(sp_dist),
  636. ctypes.byref(lp_dist))
  637. # instantiate the Point class
  638. point = Point(cx.value, cy.value, cz.value)
  639. point.is2D = self.is2D
  640. return point, dist.value, sp_dist.value, lp_dist.value
  641. def get_first_cat(self):
  642. """Fetches FIRST category number for given vector line and field, using
  643. the ``Vect_get_line_cat`` C function.
  644. .. warning::
  645. Not implemented yet.
  646. """
  647. # TODO: add this method.
  648. libvect.Vect_get_line_cat(self.map, self.id, self.field)
  649. pass
  650. def pop(self, indx):
  651. """Return the point in the index position and remove from the Line.
  652. :param indx: the index where add new point
  653. :type indx: int
  654. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  655. >>> midle_pnt = line.pop(1)
  656. >>> midle_pnt
  657. Point(1.000000, 1.000000)
  658. >>> line
  659. Line([Point(0.000000, 0.000000), Point(2.000000, 2.000000)])
  660. """
  661. if indx < 0: # Handle negative indices
  662. indx += self.c_points.contents.n_points
  663. if indx >= self.c_points.contents.n_points:
  664. raise IndexError('Index out of range')
  665. pnt = self.__getitem__(indx)
  666. libvect.Vect_line_delete_point(self.c_points, indx)
  667. return pnt
  668. def delete(self, indx):
  669. """Remove the point in the index position.
  670. :param indx: the index where add new point
  671. :type indx: int
  672. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  673. >>> line.delete(-1)
  674. >>> line
  675. Line([Point(0.000000, 0.000000), Point(1.000000, 1.000000)])
  676. """
  677. if indx < 0: # Handle negative indices
  678. indx += self.c_points.contents.n_points
  679. if indx >= self.c_points.contents.n_points:
  680. raise IndexError('Index out of range')
  681. libvect.Vect_line_delete_point(self.c_points, indx)
  682. def prune(self):
  683. """Remove duplicate points, i.e. zero length segments, using
  684. `Vect_line_prune` C function. ::
  685. >>> line = Line([(0, 0), (1, 1), (1, 1), (2, 2)])
  686. >>> line.prune()
  687. >>> line #doctest: +NORMALIZE_WHITESPACE
  688. Line([Point(0.000000, 0.000000),
  689. Point(1.000000, 1.000000),
  690. Point(2.000000, 2.000000)])
  691. ..
  692. """
  693. libvect.Vect_line_prune(self.c_points)
  694. def prune_thresh(self, threshold):
  695. """Remove points in threshold, using the ``Vect_line_prune_thresh``
  696. C funtion.
  697. :param threshold: the threshold value where prune points
  698. :type threshold: num
  699. >>> line = Line([(0, 0), (1.0, 1.0), (1.2, 0.9), (2, 2)])
  700. >>> line.prune_thresh(0.5)
  701. >>> line #doctest: +SKIP +NORMALIZE_WHITESPACE
  702. Line([Point(0.000000, 0.000000),
  703. Point(1.000000, 1.000000),
  704. Point(2.000000, 2.000000)])
  705. .. warning ::
  706. prune_thresh is not working yet.
  707. """
  708. libvect.Vect_line_prune(self.c_points, ctypes.c_double(threshold))
  709. def remove(self, pnt):
  710. """Delete point at given index and move all points above down, using
  711. `Vect_line_delete_point` C function.
  712. :param pnt: the point to remove
  713. :type pnt: a Point object or a tuple with the coordinates
  714. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  715. >>> line.remove((2, 2))
  716. >>> line[-1]
  717. Point(1.000000, 1.000000)
  718. ..
  719. """
  720. for indx, point in enumerate(self.__iter__()):
  721. if pnt == point:
  722. libvect.Vect_line_delete_point(self.c_points, indx)
  723. return
  724. raise ValueError('list.remove(x): x not in list')
  725. def reverse(self):
  726. """Reverse the order of vertices, using `Vect_line_reverse`
  727. C function. ::
  728. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  729. >>> line.reverse()
  730. >>> line #doctest: +NORMALIZE_WHITESPACE
  731. Line([Point(2.000000, 2.000000),
  732. Point(1.000000, 1.000000),
  733. Point(0.000000, 0.000000)])
  734. ..
  735. """
  736. libvect.Vect_line_reverse(self.c_points)
  737. def segment(self, start, end):
  738. # TODO improve documentation
  739. """Create line segment. using the ``Vect_line_segment`` C function.
  740. """
  741. line = Line()
  742. libvect.Vect_line_segment(self.c_points, start, end, line.c_points)
  743. return line
  744. def tolist(self):
  745. """Return a list of tuple. ::
  746. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  747. >>> line.tolist()
  748. [(0.0, 0.0), (1.0, 1.0), (2.0, 0.0), (1.0, -1.0)]
  749. ..
  750. """
  751. return [pnt.coords() for pnt in self.__iter__()]
  752. def toarray(self):
  753. """Return an array of coordinates. ::
  754. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  755. >>> line.toarray() #doctest: +NORMALIZE_WHITESPACE
  756. array([[ 0., 0.],
  757. [ 1., 1.],
  758. [ 2., 0.],
  759. [ 1., -1.]])
  760. ..
  761. """
  762. return np.array(self.tolist())
  763. def get_wkt(self):
  764. """Return a Well Known Text string of the line. ::
  765. >>> line = Line([(0, 0), (1, 1), (1, 2)])
  766. >>> line.get_wkt() #doctest: +ELLIPSIS
  767. 'LINESTRING(0.000000 0.000000, ..., 1.000000 2.000000)'
  768. ..
  769. """
  770. return "LINESTRING(%s)" % ', '.join([
  771. ' '.join(['%f' % coord for coord in pnt.coords()])
  772. for pnt in self.__iter__()])
  773. def from_wkt(self, wkt):
  774. """Create a line reading a WKT string.
  775. :param wkt: the WKT string containing the LINESTRING
  776. :type wkt: str
  777. >>> line = Line()
  778. >>> line.from_wkt("LINESTRING(0 0,1 1,1 2)")
  779. >>> line #doctest: +NORMALIZE_WHITESPACE
  780. Line([Point(0.000000, 0.000000),
  781. Point(1.000000, 1.000000),
  782. Point(1.000000, 2.000000)])
  783. ..
  784. """
  785. match = re.match('LINESTRING\((.*)\)', wkt)
  786. if match:
  787. self.reset()
  788. for coord in match.groups()[0].strip().split(','):
  789. self.append(tuple([float(e) for e in coord.split(' ')]))
  790. else:
  791. return None
  792. def get_wkb(self):
  793. """Return a WKB buffer.
  794. .. warning::
  795. Not implemented yet.
  796. """
  797. pass
  798. def buffer(self, dist=None, dist_x=None, dist_y=None,
  799. angle=0, round_=True, caps=True, tol=0.1):
  800. """Return the buffer area around the line, using the
  801. ``Vect_line_buffer2`` C function.
  802. :param dist: the distance around the line
  803. :type dist: num
  804. :param dist_x: the distance along x
  805. :type dist_x: num
  806. :param dist_y: the distance along y
  807. :type dist_y: num
  808. :param angle: the angle between 0x and major axis
  809. :type angle: num
  810. :param round_: to make corners round
  811. :type round_: bool
  812. :param tol: fix the maximum distance between theoretical arc and
  813. output segments
  814. :type tol: float
  815. :returns: the buffer as Area object
  816. >>> line = Line([(0, 0), (0, 2)])
  817. >>> area = line.buffer(10)
  818. >>> area.boundary #doctest: +ELLIPSIS
  819. Line([Point(-10.000000, 0.000000),...Point(-10.000000, 0.000000)])
  820. >>> area.centroid
  821. Point(0.000000, 0.000000)
  822. >>> area.isles
  823. []
  824. ..
  825. """
  826. if dist is not None:
  827. dist_x = dist
  828. dist_y = dist
  829. elif not dist_x or not dist_y:
  830. raise TypeError('TypeError: buffer expected 1 arguments, got 0')
  831. p_bound = ctypes.pointer(ctypes.pointer(libvect.line_pnts()))
  832. pp_isle = ctypes.pointer(ctypes.pointer(
  833. ctypes.pointer(libvect.line_pnts())))
  834. n_isles = ctypes.pointer(ctypes.c_int())
  835. libvect.Vect_line_buffer2(self.c_points,
  836. dist_x, dist_y, angle,
  837. int(round_), int(caps), tol,
  838. p_bound, pp_isle, n_isles)
  839. return Area(boundary=Line(c_points=p_bound.contents),
  840. centroid=self[0],
  841. isles=[Line(c_points=pp_isle[i].contents)
  842. for i in range(n_isles.contents.value)])
  843. def reset(self):
  844. """Reset line, using `Vect_reset_line` C function. ::
  845. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  846. >>> len(line)
  847. 4
  848. >>> line.reset()
  849. >>> len(line)
  850. 0
  851. >>> line
  852. Line([])
  853. ..
  854. """
  855. libvect.Vect_reset_line(self.c_points)
  856. class Node(object):
  857. pass
  858. class Boundary(Line):
  859. """
  860. """
  861. # geometry type
  862. gtype = libvect.GV_BOUNDARY
  863. def __init__(self, lines=None, left=None, right=None,
  864. **kargs):
  865. v_id = kargs.get('v_id', 0)
  866. self.dir = libvect.GV_FORWARD if v_id > 0 else libvect.GV_BACKWARD
  867. super(Boundary, self).__init__(**kargs)
  868. self.c_left = ctypes.pointer(ctypes.c_int())
  869. self.c_right = ctypes.pointer(ctypes.c_int())
  870. #self.get_left_right()
  871. @property
  872. def left_id(self):
  873. return self.c_left.contents.value
  874. @property
  875. def right_id(self):
  876. return self.c_right.contents.value
  877. def __repr__(self):
  878. return "Boundary(v_id=%r)" % self.id
  879. def _get_centroid(self, side, idonly=False):
  880. if side > 0:
  881. v_id = libvect.Vect_get_area_centroid(self.c_mapinfo, side)
  882. v_id = v_id if v_id else None
  883. if idonly:
  884. return v_id
  885. else:
  886. cntr = Centroid(v_id=v_id, c_mapinfo=self.c_mapinfo)
  887. return cntr
  888. def get_left_centroid(self, idonly=False):
  889. """Return left value
  890. :param idonly: True to return only the cat of feature
  891. :type idonly: bool
  892. """
  893. return self._get_centroid(self.left_id, idonly)
  894. def get_right_centroid(self, idonly=False):
  895. """Return right value
  896. :param idonly: True to return only the cat of feature
  897. :type idonly: bool
  898. """
  899. return self._get_centroid(self.left_id, idonly)
  900. def get_left_right(self):
  901. """Return left and right value"""
  902. libvect.Vect_get_line_areas(self.c_mapinfo, self.id,
  903. self.c_left, self.c_right)
  904. return self.c_left.contents.value, self.c_right.contents.value
  905. def area(self):
  906. """Return the area of the polygon.
  907. >>> bound = Boundary(points=[(0, 0), (0, 2), (2, 2), (2, 0),
  908. ... (0, 0)])
  909. >>> bound.area()
  910. 4.0
  911. """
  912. libgis.G_begin_polygon_area_calculations()
  913. return libgis.G_area_of_polygon(self.c_points.contents.x,
  914. self.c_points.contents.y,
  915. self.c_points.contents.n_points)
  916. class Centroid(Point):
  917. """The Centroid class inherit from the Point class.
  918. Centroid contains an attribute with the C Map_info struct, and attributes
  919. with the id of the Area. ::
  920. >>> centroid = Centroid(x=0, y=10)
  921. >>> centroid
  922. Centoid(0.000000, 10.000000)
  923. >>> from grass.pygrass.vector import VectorTopo
  924. >>> geo = VectorTopo('geology')
  925. >>> geo.open()
  926. >>> centroid = Centroid(v_id=1, c_mapinfo=geo.c_mapinfo)
  927. >>> centroid
  928. Centoid(893202.874416, 297339.312795)
  929. ..
  930. """
  931. # geometry type
  932. gtype = libvect.GV_CENTROID
  933. def __init__(self, area_id=None, **kargs):
  934. super(Centroid, self).__init__(**kargs)
  935. self.area_id = area_id
  936. if self.id and self.c_mapinfo and self.area_id is None:
  937. self.area_id = self.get_area_id()
  938. elif self.c_mapinfo and self.area_id and self.id is None:
  939. self.id = self.get_centroid_id()
  940. if self.area_id is not None:
  941. self.read()
  942. #self.c_pline = ctypes.pointer(libvect.P_line()) if topology else None
  943. def __repr__(self):
  944. return "Centoid(%s)" % ', '.join(['%f' % co for co in self.coords()])
  945. def get_centroid_id(self):
  946. """Return the centroid_id, using the c_mapinfo and an area_id
  947. attributes of the class, and calling the Vect_get_area_centroid
  948. C function, if no centroid_id were found return None"""
  949. centroid_id = libvect.Vect_get_area_centroid(self.c_mapinfo,
  950. self.area_id)
  951. return centroid_id if centroid_id != 0 else None
  952. def get_area_id(self):
  953. """Return the area_id, using the c_mapinfo and an centroid_id
  954. attributes of the class, and calling the Vect_get_centroid_area
  955. C function, if no area_id were found return None"""
  956. area_id = libvect.Vect_get_centroid_area(self.c_mapinfo,
  957. self.id)
  958. return area_id if area_id != 0 else None
  959. class Isle(Geo):
  960. """An Isle is an area contained by another area.
  961. """
  962. def __init__(self, **kargs):
  963. super(Isle, self).__init__(**kargs)
  964. #self.area_id = area_id
  965. def __repr__(self):
  966. return "Isle(%d)" % (self.id)
  967. def boundaries(self):
  968. """Return a list of boundaries"""
  969. ilist = Ilist()
  970. libvect.Vect_get_isle_boundaries(self.c_mapinfo, self.id,
  971. ilist.c_ilist)
  972. return ilist
  973. def bbox(self, bbox=None):
  974. """Return bounding box of Isle"""
  975. bbox = bbox if bbox else Bbox()
  976. libvect.Vect_get_isle_box(self.c_mapinfo, self.id, bbox.c_bbox)
  977. return bbox
  978. def points(self):
  979. """Return a Line object with the outer ring points"""
  980. line = Line()
  981. libvect.Vect_get_isle_points(self.c_mapinfo, self.id, line.c_points)
  982. return line
  983. def points_geos(self):
  984. """Return a Line object with the outer ring points
  985. """
  986. return libvect.Vect_get_isle_points_geos(self.c_mapinfo, self.id)
  987. def area_id(self):
  988. """Returns area id for isle."""
  989. return libvect.Vect_get_isle_area(self.c_mapinfo, self.id)
  990. def alive(self):
  991. """Check if isle is alive or dead (topology required)"""
  992. return bool(libvect.Vect_isle_alive(self.c_mapinfo, self.id))
  993. def contain_pnt(self, pnt):
  994. """Check if point is in area.
  995. :param pnt: the point to remove
  996. :type pnt: a Point object or a tuple with the coordinates
  997. """
  998. bbox = self.bbox()
  999. return bool(libvect.Vect_point_in_island(pnt.x, pnt.y,
  1000. self.c_mapinfo, self.id,
  1001. bbox.c_bbox.contents))
  1002. def area(self):
  1003. """Return the area value of an Isle"""
  1004. border = self.points()
  1005. return libgis.G_area_of_polygon(border.c_points.contents.x,
  1006. border.c_points.contents.y,
  1007. border.c_points.contents.n_points)
  1008. def perimeter(self):
  1009. """Return the perimeter value of an Isle.
  1010. ::
  1011. double Vect_area_perimeter()
  1012. """
  1013. border = self.points()
  1014. return libvect.Vect_area_perimeter(border.c_points)
  1015. class Isles(object):
  1016. def __init__(self, c_mapinfo, area_id=None):
  1017. self.c_mapinfo = c_mapinfo
  1018. self.area_id = area_id
  1019. self._isles_id = None
  1020. self._isles = None
  1021. if area_id:
  1022. self._isles_id = self.get_isles_id()
  1023. self._isles = self.get_isles()
  1024. def __len__(self):
  1025. return libvect.Vect_get_area_num_isles(self.c_mapinfo, self.area_id)
  1026. def __repr__(self):
  1027. return "Isles(%r)" % self.area_id
  1028. def __getitem__(self, key):
  1029. if self._isles is None:
  1030. self.get_isles()
  1031. return self._isles[key]
  1032. def get_isles_id(self):
  1033. """Return the id of isles"""
  1034. return [libvect.Vect_get_area_isle(self.c_mapinfo, self.area_id, i)
  1035. for i in range(self.__len__())]
  1036. def get_isles(self):
  1037. """Return isles"""
  1038. return [Isle(v_id=isle_id, c_mapinfo=self.c_mapinfo)
  1039. for isle_id in self._isles_id]
  1040. def select_by_bbox(self, bbox):
  1041. """Vect_select_isles_by_box
  1042. .. warning::
  1043. Not implemented yet.
  1044. """
  1045. pass
  1046. class Area(Geo):
  1047. """
  1048. 'Vect_build_line_area',
  1049. 'Vect_find_area',
  1050. 'Vect_get_area_box',
  1051. 'Vect_get_area_points_geos',
  1052. 'Vect_get_centroid_area',
  1053. 'Vect_get_isle_area',
  1054. 'Vect_get_line_areas',
  1055. 'Vect_get_num_areas',
  1056. 'Vect_get_point_in_area',
  1057. 'Vect_isle_find_area',
  1058. 'Vect_point_in_area',
  1059. 'Vect_point_in_area_outer_ring',
  1060. 'Vect_read_area_geos',
  1061. 'Vect_remove_small_areas',
  1062. 'Vect_select_areas_by_box',
  1063. 'Vect_select_areas_by_polygon']
  1064. """
  1065. # geometry type
  1066. gtype = libvect.GV_AREA
  1067. def __init__(self, boundary=None, centroid=None, isles=None, **kargs):
  1068. super(Area, self).__init__(**kargs)
  1069. self.boundary = None
  1070. self.centroid = None
  1071. self.isles = None
  1072. if boundary and centroid:
  1073. self.boundary = boundary
  1074. self.centroid = centroid
  1075. self.isles = isles if isles else []
  1076. # set the attributes
  1077. if self.attrs and self.cat:
  1078. self.attrs.cat = self.cat
  1079. def __repr__(self):
  1080. return "Area(%d)" % self.id if self.id else "Area( )"
  1081. def init_from_id(self, area_id=None):
  1082. """Return an Area object"""
  1083. if area_id is None and self.id is None:
  1084. raise ValueError("You need to give or set the area_id")
  1085. self.id = area_id if area_id is not None else self.id
  1086. # get boundary
  1087. self.get_points()
  1088. # get isles
  1089. self.get_isles()
  1090. def get_points(self, line=None):
  1091. """Return a Line object with the outer ring
  1092. :param line: a Line object to fill with info from points of area
  1093. :type line: a Line object
  1094. """
  1095. line = Line() if line is None else line
  1096. libvect.Vect_get_area_points(self.c_mapinfo, self.id, line.c_points)
  1097. return line
  1098. def get_centroid(self, centroid=None):
  1099. """Return the centroid
  1100. :param centroid: a Centroid object to fill with info from centroid of area
  1101. :type centroid: a Centroid object
  1102. """
  1103. centroid_id = libvect.Vect_get_area_centroid(self.c_mapinfo, self.id)
  1104. if centroid:
  1105. centroid.id = centroid_id
  1106. centroid.read()
  1107. return centroid
  1108. return Centroid(v_id=centroid_id, c_mapinfo=self.c_mapinfo,
  1109. area_id=self.id)
  1110. def num_isles(self):
  1111. return libvect.Vect_get_area_num_isles(self.c_mapinfo, self.id)
  1112. def get_isles(self, isles=None):
  1113. """Instantiate the boundary attribute reading area_id"""
  1114. if isles is not None:
  1115. isles.area_id = self.id
  1116. return isles
  1117. return Isles(self.c_mapinfo, self.id)
  1118. def area(self):
  1119. """Returns area of area without areas of isles.
  1120. double Vect_get_area_area (const struct Map_info *Map, int area)
  1121. """
  1122. return libvect.Vect_get_area_area(self.c_mapinfo, self.id)
  1123. def alive(self):
  1124. """Check if area is alive or dead (topology required)
  1125. """
  1126. return bool(libvect.Vect_area_alive(self.c_mapinfo, self.id))
  1127. def bbox(self, bbox=None):
  1128. """Return the Bbox of area
  1129. :param bbox: a Bbox object to fill with info from bounding box of area
  1130. :type bbox: a Bbox object
  1131. """
  1132. bbox = bbox if bbox else Bbox()
  1133. libvect.Vect_get_area_box(self.c_mapinfo, self.id, bbox.c_bbox)
  1134. return bbox
  1135. def buffer(self, dist=None, dist_x=None, dist_y=None,
  1136. angle=0, round_=True, caps=True, tol=0.1):
  1137. """Return the buffer area around the area, using the
  1138. ``Vect_area_buffer2`` C function.
  1139. :param dist: the distance around the area
  1140. :type dist: num
  1141. :param dist_x: the distance along x
  1142. :type dist_x: num
  1143. :param dist_y: the distance along y
  1144. :type dist_y: num
  1145. :param angle: the angle between 0x and major axis
  1146. :type angle: num
  1147. :param round_: to make corners round
  1148. :type round_: bool
  1149. :param tol: fix the maximum distance between theoretical arc and
  1150. output segments
  1151. :type tol: float
  1152. :returns: the buffer as Area object
  1153. """
  1154. if dist is not None:
  1155. dist_x = dist
  1156. dist_y = dist
  1157. elif not dist_x or not dist_y:
  1158. raise TypeError('TypeError: buffer expected 1 arguments, got 0')
  1159. p_bound = ctypes.pointer(ctypes.pointer(libvect.line_pnts()))
  1160. pp_isle = ctypes.pointer(ctypes.pointer(
  1161. ctypes.pointer(libvect.line_pnts())))
  1162. n_isles = ctypes.pointer(ctypes.c_int())
  1163. libvect.Vect_area_buffer2(self.c_mapinfo, self.id,
  1164. dist_x, dist_y, angle,
  1165. int(round_), int(caps), tol,
  1166. p_bound, pp_isle, n_isles)
  1167. return Area(boundary=Line(c_points=p_bound.contents),
  1168. centroid=self.centroid,
  1169. isles=[Line(c_points=pp_isle[i].contents)
  1170. for i in range(n_isles.contents.value)])
  1171. def boundaries(self, ilist=False):
  1172. """Creates list of boundaries for given area.
  1173. int Vect_get_area_boundaries(const struct Map_info *Map,
  1174. int area, struct ilist *List)
  1175. """
  1176. ilst = Ilist()
  1177. libvect.Vect_get_area_boundaries(self.c_mapinfo, self.id,
  1178. ilst.c_ilist)
  1179. if ilist:
  1180. return ilist
  1181. return [Boundary(v_id, c_mapinfo=self.c_mapinfo) for v_id in ilst]
  1182. def cats(self, cats=None):
  1183. """Get area categories.
  1184. :param cats: a Cats object to fill with info with area categories
  1185. :type cats: a Cats object
  1186. """
  1187. cats = cats if cats else Cats()
  1188. libvect.Vect_get_area_cats(self.c_mapinfo, self.id, cats.c_cats)
  1189. return cats
  1190. def get_first_cat(self):
  1191. """Find FIRST category of given field and area.
  1192. int Vect_get_area_cat(const struct Map_info *Map, int area, int field)
  1193. """
  1194. pass
  1195. def contain_pnt(self, pnt, bbox=None):
  1196. """Check if point is in area.
  1197. :param pnt: the point to analyze
  1198. :type pnt: a Point object or a tuple with the coordinates
  1199. :param bbox: the bounding box where run the analysis
  1200. :type bbox: a Bbox object
  1201. """
  1202. bbox = bbox if bbox else self.bbox()
  1203. return bool(libvect.Vect_point_in_area(pnt.x, pnt.y,
  1204. self.c_mapinfo, self.id,
  1205. bbox.c_bbox))
  1206. def perimeter(self):
  1207. """Calculate area perimeter.
  1208. double Vect_area_perimeter (const struct line_pnts *Points)
  1209. """
  1210. border = self.get_points()
  1211. return libvect.Vect_area_perimeter(border.c_points)
  1212. def read(self, line=None, centroid=None, isles=None):
  1213. self.boundary = self.get_points(line)
  1214. self.centroid = self.get_centroid(centroid)
  1215. #self.isles = self.get_isles(isles)
  1216. libvect.Vect_read_line(self.c_mapinfo, None, self.c_cats,
  1217. self.centroid.id)
  1218. #
  1219. # Define a dictionary to convert the feature type to name and or object
  1220. #
  1221. GV_TYPE = {libvect.GV_POINT: {'label': 'point', 'obj': Point},
  1222. libvect.GV_LINE: {'label': 'line', 'obj': Line},
  1223. libvect.GV_BOUNDARY: {'label': 'boundary', 'obj': Boundary},
  1224. libvect.GV_CENTROID: {'label': 'centroid', 'obj': Centroid},
  1225. libvect.GV_FACE: {'label': 'face', 'obj': None},
  1226. libvect.GV_KERNEL: {'label': 'kernel', 'obj': None},
  1227. libvect.GV_AREA: {'label': 'area', 'obj': Area},
  1228. libvect.GV_VOLUME: {'label': 'volume', 'obj': None}, }
  1229. GEOOBJ = {"areas": Area,
  1230. "dblinks": None,
  1231. "faces": None,
  1232. "holes": None,
  1233. "islands": Isle,
  1234. "kernels": None,
  1235. "line_points": None,
  1236. "points": Point,
  1237. "lines": Line,
  1238. "nodes": Node,
  1239. "volumes": None}
  1240. def c_read_next_line(c_mapinfo, c_points, c_cats):
  1241. v_id = c_mapinfo.contents.next_line
  1242. v_id = v_id if v_id != 0 else None
  1243. ftype = libvect.Vect_read_next_line(c_mapinfo, c_points, c_cats)
  1244. if ftype == -2:
  1245. raise StopIteration()
  1246. if ftype == -1:
  1247. raise
  1248. return ftype, v_id, c_points, c_cats
  1249. def read_next_line(c_mapinfo, table=None, writable=False,
  1250. c_points=None, c_cats=None, is2D=True):
  1251. """Return the next geometry feature of a vector map."""
  1252. c_points = c_points if c_points else ctypes.pointer(libvect.line_pnts())
  1253. c_cats = c_cats if c_cats else ctypes.pointer(libvect.line_cats())
  1254. ftype, v_id, c_points, c_cats = c_read_next_line(c_mapinfo, c_points,
  1255. c_cats)
  1256. return GV_TYPE[ftype]['obj'](v_id=v_id, c_mapinfo=c_mapinfo,
  1257. c_points=c_points, c_cats=c_cats,
  1258. table=table, writable=writable, is2D=is2D)
  1259. def c_read_line(feature_id, c_mapinfo, c_points, c_cats):
  1260. nmax = libvect.Vect_get_num_lines(c_mapinfo)
  1261. if feature_id < 0: # Handle negative indices
  1262. feature_id += nmax + 1
  1263. if feature_id > nmax:
  1264. raise IndexError('Index out of range')
  1265. if feature_id > 0:
  1266. ftype = libvect.Vect_read_line(c_mapinfo, c_points, c_cats, feature_id)
  1267. return ftype, c_points, c_cats
  1268. else:
  1269. raise ValueError('The index must be >0, %r given.' % feature_id)
  1270. def read_line(feature_id, c_mapinfo, table=None, writable=False,
  1271. c_points=None, c_cats=None, is2D=True):
  1272. """Return a geometry object given the feature id and the c_mapinfo.
  1273. """
  1274. c_points = c_points if c_points else ctypes.pointer(libvect.line_pnts())
  1275. c_cats = c_cats if c_cats else ctypes.pointer(libvect.line_cats())
  1276. ftype, c_points, c_cats = c_read_line(feature_id, c_mapinfo,
  1277. c_points, c_cats)
  1278. if GV_TYPE[ftype]['obj'] is not None:
  1279. return GV_TYPE[ftype]['obj'](v_id=feature_id, c_mapinfo=c_mapinfo,
  1280. c_points=c_points, c_cats=c_cats,
  1281. table=table, writable=writable, is2D=is2D)