geometry.py 36 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150
  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Wed Jul 18 10:46:25 2012
  4. @author: pietro
  5. """
  6. import ctypes
  7. import numpy as np
  8. import re
  9. import grass.lib.gis as libgis
  10. import grass.lib.vector as libvect
  11. from basic import Ilist, Bbox, Cats
  12. WKT = {'POINT\((.*)\)': 'point', # 'POINT\(\s*([+-]*\d+\.*\d*)+\s*\)'
  13. 'LINESTRING\((.*)\)': 'line'}
  14. def read_WKT(string):
  15. """Read the string and return a geometry object
  16. WKT:
  17. POINT(0 0)
  18. LINESTRING(0 0,1 1,1 2)
  19. POLYGON((0 0,4 0,4 4,0 4,0 0),(1 1, 2 1, 2 2, 1 2,1 1))
  20. MULTIPOINT(0 0,1 2)
  21. MULTILINESTRING((0 0,1 1,1 2),(2 3,3 2,5 4))
  22. MULTIPOLYGON(((0 0,4 0,4 4,0 4,0 0),(1 1,2 1,2 2,1 2,1 1)),
  23. ((-1 -1,-1 -2,-2 -2,-2 -1,-1 -1)))
  24. GEOMETRYCOLLECTION(POINT(2 3),LINESTRING(2 3,3 4))
  25. EWKT:
  26. POINT(0 0 0) -- XYZ
  27. SRID=32632;POINT(0 0) -- XY with SRID
  28. POINTM(0 0 0) -- XYM
  29. POINT(0 0 0 0) -- XYZM
  30. SRID=4326;MULTIPOINTM(0 0 0,1 2 1) -- XYM with SRID
  31. MULTILINESTRING((0 0 0,1 1 0,1 2 1),(2 3 1,3 2 1,5 4 1))
  32. 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))
  33. MULTIPOLYGON(((0 0 0,4 0 0,4 4 0,0 4 0,0 0 0),
  34. (1 1 0,2 1 0,2 2 0,1 2 0,1 1 0)),
  35. ((-1 -1 0,-1 -2 0,-2 -2 0,-2 -1 0,-1 -1 0)))
  36. GEOMETRYCOLLECTIONM( POINTM(2 3 9), LINESTRINGM(2 3 4, 3 4 5) )
  37. MULTICURVE( (0 0, 5 5), CIRCULARSTRING(4 0, 4 4, 8 4) )
  38. POLYHEDRALSURFACE( ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)),
  39. ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)),
  40. ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)),
  41. ((1 1 0, 1 1 1, 1 0 1, 1 0 0, 1 1 0)),
  42. ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0)),
  43. ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)) )
  44. TRIANGLE ((0 0, 0 9, 9 0, 0 0))
  45. 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)) )
  46. """
  47. for regexp, obj in WKT.items():
  48. if re.match(regexp, string):
  49. geo = 10
  50. return obj(geo)
  51. def read_WKB(buff):
  52. """Read the binary buffer and return a geometry object"""
  53. pass
  54. #=============================================
  55. # GEOMETRY
  56. #=============================================
  57. def get_xyz(pnt):
  58. """Return a tuple with: x, y, z. ::
  59. >>> pnt = Point(0, 0)
  60. >>> get_xyz(pnt)
  61. (0.0, 0.0, 0.0)
  62. >>> get_xyz((1, 1))
  63. (1, 1, 0.0)
  64. >>> get_xyz((1, 1, 2))
  65. (1, 1, 2)
  66. >>> get_xyz((1, 1, 2, 2)) #doctest: +ELLIPSIS
  67. Traceback (most recent call last):
  68. ...
  69. ValueError: The the format of the point is not supported: (1, 1, 2, 2)
  70. ..
  71. """
  72. if isinstance(pnt, Point):
  73. if pnt.is2D:
  74. x, y = pnt.x, pnt.y
  75. z = 0.
  76. else:
  77. x, y, z = pnt.x, pnt.y, pnt.z
  78. else:
  79. if len(pnt) == 2:
  80. x, y = pnt
  81. z = 0.
  82. elif len(pnt) == 3:
  83. x, y, z = pnt
  84. else:
  85. str_error = "The the format of the point is not supported: {0!r}"
  86. raise ValueError(str_error.format(pnt))
  87. return x, y, z
  88. class Geo(object):
  89. """
  90. >>> geo0 = Geo()
  91. >>> points = ctypes.pointer(libvect.line_pnts())
  92. >>> cats = ctypes.pointer(libvect.line_cats())
  93. >>> geo1 = Geo(c_points=points, c_cats=cats)
  94. """
  95. def __init__(self, v_id=None, c_mapinfo=None, c_points=None, c_cats=None):
  96. self.id = v_id # vector id
  97. self.c_mapinfo = c_mapinfo
  98. # set c_points
  99. if c_points is None:
  100. self.c_points = ctypes.pointer(libvect.line_pnts())
  101. else:
  102. self.c_points = c_points
  103. # set c_cats
  104. if c_cats is None:
  105. self.c_cats = ctypes.pointer(libvect.line_cats())
  106. else:
  107. self.c_cats = c_cats
  108. def is_with_topology(self):
  109. if self.c_mapinfo is not None:
  110. return self.c_mapinfo.contents.level == 2
  111. else:
  112. return False
  113. def read(self):
  114. """Read and set the coordinates of the centroid from the vector map,
  115. using the centroid_id and calling the Vect_read_line C function"""
  116. libvect.Vect_read_line(self.c_mapinfo, self.c_points,
  117. self.c_cats, self.id)
  118. def write(self):
  119. """Write the centroid to the Map."""
  120. libvect.Vect_write_line(self.c_mapinfo, libvect.GV_CENTROID,
  121. self.c_points, self.c_cats)
  122. class Point(Geo):
  123. """Instantiate a Point object that could be 2 or 3D, default
  124. parameters are 0.
  125. ::
  126. >>> pnt = Point()
  127. >>> pnt.x
  128. 0.0
  129. >>> pnt.y
  130. 0.0
  131. >>> pnt.z
  132. >>> pnt.is2D
  133. True
  134. >>> pnt
  135. Point(0.000000, 0.000000)
  136. >>> pnt.z = 0
  137. >>> pnt.is2D
  138. False
  139. >>> pnt
  140. Point(0.000000, 0.000000, 0.000000)
  141. >>> print pnt
  142. POINT(0.000000, 0.000000, 0.000000)
  143. ..
  144. """
  145. def __init__(self, x=0, y=0, z=None, is2D=True, **kargs):
  146. super(Point, self).__init__(**kargs)
  147. if self.id is not None:
  148. self.read()
  149. self.is2D = is2D
  150. else:
  151. self.is2D = True if z is None else False
  152. z = z if z is not None else 0
  153. libvect.Vect_append_point(self.c_points, x, y, z)
  154. # geometry type
  155. self.gtype = libvect.GV_POINT
  156. def _get_x(self):
  157. return self.c_points.contents.x[0]
  158. def _set_x(self, value):
  159. self.c_points.contents.x[0] = value
  160. x = property(fget=_get_x, fset=_set_x)
  161. def _get_y(self):
  162. return self.c_points.contents.y[0]
  163. def _set_y(self, value):
  164. self.c_points.contents.y[0] = value
  165. y = property(fget=_get_y, fset=_set_y)
  166. def _get_z(self):
  167. if self.is2D:
  168. return None
  169. return self.c_points.contents.z[0]
  170. def _set_z(self, value):
  171. if value is None:
  172. self.is2D = True
  173. self.c_points.contents.z[0] = 0
  174. else:
  175. self.c_points.contents.z[0] = value
  176. self.is2D = False
  177. z = property(fget=_get_z, fset=_set_z)
  178. def __str__(self):
  179. return self.get_wkt()
  180. def __repr__(self):
  181. return "Point(%s)" % ', '.join(['%f' % coor for coor in self.coords()])
  182. def __eq__(self, pnt):
  183. if isinstance(pnt, Point):
  184. return pnt.coords() == self.coords()
  185. return Point(*pnt).coords() == self.coords()
  186. def coords(self):
  187. """Return a tuple with the point coordinates. ::
  188. >>> pnt = Point(10, 100)
  189. >>> pnt.coords()
  190. (10.0, 100.0)
  191. If the point is 2D return a x, y tuple. But if we change the ``z``
  192. the Point object become a 3D point, therefore the method return a
  193. x, y, z tuple. ::
  194. >>> pnt.z = 1000.
  195. >>> pnt.coords()
  196. (10.0, 100.0, 1000.0)
  197. ..
  198. """
  199. if self.is2D:
  200. return self.x, self.y
  201. else:
  202. return self.x, self.y, self.z
  203. def get_wkt(self):
  204. """Return a "well know text" (WKT) geometry string. ::
  205. >>> pnt = Point(10, 100)
  206. >>> pnt.get_wkt()
  207. 'POINT(10.000000, 100.000000)'
  208. .. warning::
  209. Only ``POINT`` (2/3D) are supported, ``POINTM`` and ``POINT`` with:
  210. ``XYZM`` are not supported yet.
  211. """
  212. return "POINT(%s)" % ', '.join(['%f' % coord
  213. for coord in self.coords()])
  214. def get_wkb(self):
  215. """Return a "well know binary" (WKB) geometry buffer
  216. .. warning::
  217. Not implemented yet.
  218. """
  219. pass
  220. def distance(self, pnt):
  221. """Calculate distance of 2 points, using the Vect_points_distance
  222. C function, If one of the point have z == None, return the 2D distance.
  223. ::
  224. >>> pnt0 = Point(0, 0, 0)
  225. >>> pnt1 = Point(1, 0)
  226. >>> pnt0.distance(pnt1)
  227. 1.0
  228. >>> pnt1.z = 1
  229. >>> pnt1
  230. Point(1.000000, 0.000000, 1.000000)
  231. >>> pnt0.distance(pnt1)
  232. 1.4142135623730951
  233. The distance method require a :class:Point or a tuple with
  234. the coordinates.
  235. """
  236. if self.is2D or pnt.is2D:
  237. return libvect.Vect_points_distance(self.x, self.y, 0,
  238. pnt.x, pnt.y, 0, 0)
  239. else:
  240. return libvect.Vect_points_distance(self.x, self.y, self.z,
  241. pnt.x, pnt.y, pnt.z, 1)
  242. def buffer(self, dist=None, dist_x=None, dist_y=None, angle=0,
  243. round_=True, tol=0.1):
  244. """Return an Area object using the ``Vect_point_buffer2`` C function.
  245. Creates buffer around the point (px, py).
  246. """
  247. print "Not implemented yet"
  248. raise
  249. if dist is not None:
  250. dist_x = dist
  251. dist_y = dist
  252. area = Area()
  253. libvect.Vect_point_buffer2(self.x, self.y,
  254. dist_x, dist_y,
  255. angle, int(round_), tol,
  256. area.c_points)
  257. return area
  258. class Line(Geo):
  259. """Instantiate a new Line with a list of tuple, or with a list of Point. ::
  260. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  261. >>> line #doctest: +NORMALIZE_WHITESPACE
  262. Line([Point(0.000000, 0.000000),
  263. Point(1.000000, 1.000000),
  264. Point(2.000000, 0.000000),
  265. Point(1.000000, -1.000000)])
  266. ..
  267. """
  268. def __init__(self, points=None, is2D=True, **kargs):
  269. super(Line, self).__init__(**kargs)
  270. if points is not None:
  271. for pnt in points:
  272. self.append(pnt)
  273. self.is2D = is2D
  274. # geometry type
  275. self.gtype = libvect.GV_LINE
  276. def __getitem__(self, key):
  277. """Get line point of given index, slice allowed. ::
  278. >>> line = Line([(0, 0), (1, 1), (2, 2), (3, 3)])
  279. >>> line[1]
  280. Point(1.000000, 1.000000)
  281. >>> line[-1]
  282. Point(3.000000, 3.000000)
  283. >>> line[:2]
  284. [Point(0.000000, 0.000000), Point(1.000000, 1.000000)]
  285. ..
  286. """
  287. #TODO:
  288. # line[0].x = 10 is not working
  289. #pnt.c_px = ctypes.pointer(self.c_points.contents.x[indx])
  290. # pnt.c_px = ctypes.cast(id(self.c_points.contents.x[indx]),
  291. # ctypes.POINTER(ctypes.c_double))
  292. if isinstance(key, slice):
  293. #import pdb; pdb.set_trace()
  294. #Get the start, stop, and step from the slice
  295. return [Point(self.c_points.contents.x[indx],
  296. self.c_points.contents.y[indx],
  297. None if self.is2D else self.c_points.contents.z[indx])
  298. for indx in xrange(*key.indices(len(self)))]
  299. elif isinstance(key, int):
  300. if key < 0: # Handle negative indices
  301. key += self.c_points.contents.n_points
  302. if key >= self.c_points.contents.n_points:
  303. raise IndexError('Index out of range')
  304. return Point(self.c_points.contents.x[key],
  305. self.c_points.contents.y[key],
  306. None if self.is2D else self.c_points.contents.z[key])
  307. else:
  308. raise ValueError("Invalid argument type: %r." % key)
  309. def __setitem__(self, indx, pnt):
  310. """Change the coordinate of point. ::
  311. >>> line = Line([(0, 0), (1, 1)])
  312. >>> line[0] = (2, 2)
  313. >>> line
  314. Line([Point(2.000000, 2.000000), Point(1.000000, 1.000000)])
  315. ..
  316. """
  317. x, y, z = get_xyz(pnt)
  318. self.c_points.contents.x[indx] = x
  319. self.c_points.contents.y[indx] = y
  320. self.c_points.contents.z[indx] = z
  321. def __iter__(self):
  322. """Return a Point generator of the Line"""
  323. return (self.__getitem__(i) for i in range(self.__len__()))
  324. def __len__(self):
  325. """Return the number of points of the line."""
  326. return self.c_points.contents.n_points
  327. def __str__(self):
  328. return self.get_wkt()
  329. def __repr__(self):
  330. return "Line([%s])" % ', '.join([repr(pnt) for pnt in self.__iter__()])
  331. def get_pnt(self, distance, angle=0, slope=0):
  332. """Return a Point object on line in the specified distance, using the
  333. `Vect_point_on_line` C function.
  334. Raise a ValueError If the distance exceed the Line length. ::
  335. >>> line = Line([(0, 0), (1, 1)])
  336. >>> line.get_pnt(5) #doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
  337. Traceback (most recent call last):
  338. ...
  339. ValueError: The distance exceed the lenght of the line,
  340. that is: 1.414214
  341. >>> line.get_pnt(1)
  342. Point(0.707107, 0.707107)
  343. ..
  344. """
  345. # instantiate an empty Point object
  346. maxdist = self.length()
  347. if distance > maxdist:
  348. str_err = "The distance exceed the lenght of the line, that is: %f"
  349. raise ValueError(str_err % maxdist)
  350. pnt = Point(0, 0, -9999)
  351. libvect.Vect_point_on_line(self.c_points, distance,
  352. pnt.c_points.contents.x,
  353. pnt.c_points.contents.y,
  354. pnt.c_points.contents.z,
  355. angle, slope)
  356. pnt.is2D = self.is2D
  357. return pnt
  358. def append(self, pnt):
  359. """Appends one point to the end of a line, using the
  360. ``Vect_append_point`` C function. ::
  361. >>> line = Line()
  362. >>> line.append((10, 100))
  363. >>> line
  364. Line([Point(10.000000, 100.000000)])
  365. >>> line.append((20, 200))
  366. >>> line
  367. Line([Point(10.000000, 100.000000), Point(20.000000, 200.000000)])
  368. Like python list.
  369. """
  370. x, y, z = get_xyz(pnt)
  371. libvect.Vect_append_point(self.c_points, x, y, z)
  372. def bbox(self):
  373. """Return the bounding box of the line, using ``Vect_line_box``
  374. C function. ::
  375. >>> line = Line([(0, 0), (0, 1), (2, 1), (2, 0)])
  376. >>> bbox = line.bbox()
  377. >>> bbox
  378. Bbox(1.0, 0.0, 2.0, 0.0)
  379. ..
  380. """
  381. bbox = Bbox()
  382. libvect.Vect_line_box(self.c_points, bbox.c_bbox)
  383. return bbox
  384. def extend(self, line, forward=True):
  385. """Appends points to the end of a line.
  386. It is possible to extend a line, give a list of points, or directly
  387. with a line_pnts struct.
  388. If forward is True the line is extend forward otherwise is extend
  389. backward. The method use the `Vect_append_points` C function. ::
  390. >>> line = Line([(0, 0), (1, 1)])
  391. >>> line.extend( Line([(2, 2), (3, 3)]) )
  392. >>> line #doctest: +NORMALIZE_WHITESPACE
  393. Line([Point(0.000000, 0.000000),
  394. Point(1.000000, 1.000000),
  395. Point(2.000000, 2.000000),
  396. Point(3.000000, 3.000000)])
  397. Like python list, it is possible to extend a line, with another line
  398. or with a list of points.
  399. """
  400. # set direction
  401. if forward:
  402. direction = libvect.GV_FORWARD
  403. else:
  404. direction = libvect.GV_BACKWARD
  405. # check if is a Line object
  406. if isinstance(line, Line):
  407. c_points = line.c_points
  408. else:
  409. # instantiate a Line object
  410. lin = Line()
  411. for pnt in line:
  412. # add the points to the line
  413. lin.append(pnt)
  414. c_points = lin.c_points
  415. libvect.Vect_append_points(self.c_points, c_points, direction)
  416. def insert(self, indx, pnt):
  417. """Insert new point at index position and move all old points at
  418. that position and above up, using ``Vect_line_insert_point``
  419. C function. ::
  420. >>> line = Line([(0, 0), (1, 1)])
  421. >>> line.insert(0, Point(1.000000, -1.000000) )
  422. >>> line #doctest: +NORMALIZE_WHITESPACE
  423. Line([Point(1.000000, -1.000000),
  424. Point(0.000000, 0.000000),
  425. Point(1.000000, 1.000000)])
  426. ..
  427. """
  428. if indx < 0: # Handle negative indices
  429. indx += self.c_points.contents.n_points
  430. if indx >= self.c_points.contents.n_points:
  431. raise IndexError('Index out of range')
  432. x, y, z = get_xyz(pnt)
  433. libvect.Vect_line_insert_point(self.c_points, indx, x, y, z)
  434. def length(self):
  435. """Calculate line length, 3D-length in case of 3D vector line, using
  436. `Vect_line_length` C function. ::
  437. >>> line = Line([(0, 0), (1, 1), (0, 1)])
  438. >>> line.length()
  439. 2.414213562373095
  440. ..
  441. """
  442. return libvect.Vect_line_length(self.c_points)
  443. def length_geodesic(self):
  444. """Calculate line length, usig `Vect_line_geodesic_length` C function.
  445. ::
  446. >>> line = Line([(0, 0), (1, 1), (0, 1)])
  447. >>> line.length_geodesic()
  448. 2.414213562373095
  449. ..
  450. """
  451. return libvect.Vect_line_geodesic_length(self.c_points)
  452. def distance(self, pnt):
  453. """Return a tuple with:
  454. * the closest point on the line,
  455. * the distance between these two points,
  456. * distance of point from segment beginning
  457. * distance of point from line
  458. The distance is compute using the ``Vect_line_distance`` C function.
  459. """
  460. # instantite outputs
  461. cx = ctypes.c_double(0)
  462. cy = ctypes.c_double(0)
  463. cz = ctypes.c_double(0)
  464. dist = ctypes.c_double(0)
  465. sp_dist = ctypes.c_double(0)
  466. lp_dist = ctypes.c_double(0)
  467. libvect.Vect_line_distance(self.c_points,
  468. pnt.x, pnt.y, pnt.z, 0 if self.is2D else 1,
  469. ctypes.byref(cx), ctypes.byref(cy),
  470. ctypes.byref(cz), ctypes.byref(dist),
  471. ctypes.byref(sp_dist),
  472. ctypes.byref(lp_dist))
  473. # instantiate the Point class
  474. point = Point(cx.value, cy.value, cz.value)
  475. point.is2D = self.is2D
  476. return point, dist, sp_dist, lp_dist
  477. def get_first_cat(self):
  478. """Fetches FIRST category number for given vector line and field, using
  479. the ``Vect_get_line_cat`` C function.
  480. .. warning::
  481. Not implemented yet.
  482. """
  483. # TODO: add this method.
  484. libvect.Vect_get_line_cat(self.map, self.id, self.field)
  485. pass
  486. def pop(self, indx):
  487. """Return the point in the index position and remove from the Line. ::
  488. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  489. >>> midle_pnt = line.pop(1)
  490. >>> midle_pnt
  491. Point(1.000000, 1.000000)
  492. >>> line
  493. Line([Point(0.000000, 0.000000), Point(2.000000, 2.000000)])
  494. ..
  495. """
  496. if indx < 0: # Handle negative indices
  497. indx += self.c_points.contents.n_points
  498. if indx >= self.c_points.contents.n_points:
  499. raise IndexError('Index out of range')
  500. pnt = self.__getitem__(indx)
  501. libvect.Vect_line_delete_point(self.c_points, indx)
  502. return pnt
  503. def delete(self, indx):
  504. """Remove the point in the index position. ::
  505. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  506. >>> line.delete(-1)
  507. >>> line
  508. Line([Point(0.000000, 0.000000), Point(1.000000, 1.000000)])
  509. ..
  510. """
  511. if indx < 0: # Handle negative indices
  512. indx += self.c_points.contents.n_points
  513. if indx >= self.c_points.contents.n_points:
  514. raise IndexError('Index out of range')
  515. libvect.Vect_line_delete_point(self.c_points, indx)
  516. def prune(self):
  517. """Remove duplicate points, i.e. zero length segments, using
  518. `Vect_line_prune` C function. ::
  519. >>> line = Line([(0, 0), (1, 1), (1, 1), (2, 2)])
  520. >>> line.prune()
  521. >>> line #doctest: +NORMALIZE_WHITESPACE
  522. Line([Point(0.000000, 0.000000),
  523. Point(1.000000, 1.000000),
  524. Point(2.000000, 2.000000)])
  525. ..
  526. """
  527. libvect.Vect_line_prune(self.c_points)
  528. def prune_thresh(self, threshold):
  529. """Remove points in threshold, using the ``Vect_line_prune_thresh``
  530. C funtion. ::
  531. >>> line = Line([(0, 0), (1.0, 1.0), (1.2, 0.9), (2, 2)])
  532. >>> line.prune_thresh(0.5)
  533. >>> line #doctest: +SKIP +NORMALIZE_WHITESPACE
  534. Line([Point(0.000000, 0.000000),
  535. Point(1.000000, 1.000000),
  536. Point(2.000000, 2.000000)])
  537. .. warning ::
  538. prune_thresh is not working yet.
  539. """
  540. libvect.Vect_line_prune(self.c_points, ctypes.c_double(threshold))
  541. def remove(self, pnt):
  542. """Delete point at given index and move all points above down, using
  543. `Vect_line_delete_point` C function. ::
  544. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  545. >>> line.remove((2, 2))
  546. >>> line[-1]
  547. Point(1.000000, 1.000000)
  548. ..
  549. """
  550. for indx, point in enumerate(self.__iter__()):
  551. if pnt == point:
  552. libvect.Vect_line_delete_point(self.c_points, indx)
  553. return
  554. raise ValueError('list.remove(x): x not in list')
  555. def reverse(self):
  556. """Reverse the order of vertices, using `Vect_line_reverse`
  557. C function. ::
  558. >>> line = Line([(0, 0), (1, 1), (2, 2)])
  559. >>> line.reverse()
  560. >>> line #doctest: +NORMALIZE_WHITESPACE
  561. Line([Point(2.000000, 2.000000),
  562. Point(1.000000, 1.000000),
  563. Point(0.000000, 0.000000)])
  564. ..
  565. """
  566. libvect.Vect_line_reverse(self.c_points)
  567. def segment(self, start, end):
  568. """Create line segment. using the ``Vect_line_segment`` C function."""
  569. line = Line()
  570. libvect.Vect_line_segment(self.c_points, start, end, line.c_points)
  571. return line
  572. def tolist(self):
  573. """Return a list of tuple. ::
  574. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  575. >>> line.tolist()
  576. [(0.0, 0.0), (1.0, 1.0), (2.0, 0.0), (1.0, -1.0)]
  577. ..
  578. """
  579. return [pnt.coords() for pnt in self.__iter__()]
  580. def toarray(self):
  581. """Return an array of coordinates. ::
  582. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  583. >>> line.toarray() #doctest: +NORMALIZE_WHITESPACE
  584. array([[ 0., 0.],
  585. [ 1., 1.],
  586. [ 2., 0.],
  587. [ 1., -1.]])
  588. ..
  589. """
  590. return np.array(self.tolist())
  591. def get_wkt(self):
  592. """Return a Well Known Text string of the line. ::
  593. >>> line = Line([(0, 0), (1, 1), (1, 2)])
  594. >>> line.get_wkt() #doctest: +ELLIPSIS
  595. 'LINESTRING(0.000000 0.000000, ..., 1.000000 2.000000)'
  596. ..
  597. """
  598. return "LINESTRING(%s)" % ', '.join([
  599. ' '.join(['%f' % coord for coord in pnt.coords()])
  600. for pnt in self.__iter__()])
  601. def from_wkt(self, wkt):
  602. """Read a WKT string. ::
  603. >>> line = Line()
  604. >>> line.from_wkt("LINESTRING(0 0,1 1,1 2)")
  605. >>> line #doctest: +NORMALIZE_WHITESPACE
  606. Line([Point(0.000000, 0.000000),
  607. Point(1.000000, 1.000000),
  608. Point(1.000000, 2.000000)])
  609. ..
  610. """
  611. match = re.match('LINESTRING\((.*)\)', wkt)
  612. if match:
  613. self.reset()
  614. for coord in match.groups()[0].strip().split(','):
  615. self.append(tuple([float(e) for e in coord.split(' ')]))
  616. else:
  617. return None
  618. def get_wkb(self):
  619. """Return a WKB buffer.
  620. .. warning::
  621. Not implemented yet.
  622. """
  623. pass
  624. def buffer(self, dist=None, dist_x=None, dist_y=None,
  625. angle=0, round_=True, tol=0.1):
  626. """Return the buffer area around the line, using the
  627. ``Vect_line_buffer2`` C function.
  628. .. warning::
  629. Not implemented yet.
  630. """
  631. if dist is not None:
  632. dist_x = dist
  633. dist_y = dist
  634. area = Area()
  635. libvect.Vect_line_buffer2(self.c_points,
  636. dist_x, dist_y,
  637. angle, int(round_), tol,
  638. area.boundary.c_points,
  639. area.isles.c_points,
  640. area.num_isles)
  641. return area
  642. def reset(self):
  643. """Reset line, using `Vect_reset_line` C function. ::
  644. >>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
  645. >>> len(line)
  646. 4
  647. >>> line.reset()
  648. >>> len(line)
  649. 0
  650. >>> line
  651. Line([])
  652. ..
  653. """
  654. libvect.Vect_reset_line(self.c_points)
  655. class Node(object):
  656. pass
  657. class Boundary(Line):
  658. """
  659. """
  660. def __init__(self, area_id=None, lines=None, left=None, right=None,
  661. **kargs):
  662. super(Boundary, self).__init__(**kargs)
  663. self.area_id = area_id
  664. self.ilist = Ilist()
  665. self.lines = lines
  666. if lines:
  667. if len(lines) != len(left) or len(lines) != len(right):
  668. str_err = "Left and right must have the same length of lines"
  669. raise ValueError(str_err)
  670. self.left = Ilist()
  671. self.right = Ilist()
  672. # geometry type
  673. self.gtype = libvect.GV_BOUNDARY
  674. def __repr__(self):
  675. return "Boundary(v_id=%r)" % self.id
  676. def boundaries(self):
  677. """Returna Ilist object with the line id"""
  678. bounds = Ilist()
  679. libvect.Vect_get_area_boundaries(self.c_mapinfo, self.area_id,
  680. bounds.c_ilist)
  681. return bounds
  682. def get_left_right(self):
  683. """Return left and right value"""
  684. left = ctypes.poiter(ctypes.c_int())
  685. right = ctypes.poiter(ctypes.c_int())
  686. libvect.Vect_get_line_areas(self.c_mapinfo, self.id,
  687. left, right)
  688. return left.contents.value, right.contents.value
  689. class Centroid(Point):
  690. """The Centroid class inherit from the Point class.
  691. Centroid contains an attribute with the C Map_info struct, and attributes
  692. with the id of the Area. ::
  693. >>> centroid = Centroid(x=0, y=10)
  694. >>> centroid
  695. Centoid(0.000000, 10.000000)
  696. >>> import pygrass
  697. >>> mun = pygrass.vector.VectorTopo('boundary_municp_sqlite')
  698. >>> mun.open()
  699. >>> centroid = Centroid(v_id=5129, c_mapinfo=mun.c_mapinfo)
  700. >>> centroid
  701. Centoid(463784.493822, 311023.913274)
  702. ..
  703. """
  704. def __init__(self, area_id=None, **kargs):
  705. super(Centroid, self).__init__(**kargs)
  706. self.area_id = area_id
  707. if self.id and self.c_mapinfo and self.area_id is None:
  708. self.area_id = self.get_area_id()
  709. elif self.c_mapinfo and self.area_id and self.id is None:
  710. self.id = self.get_centroid_id()
  711. if self.area_id is not None:
  712. self.cats = Cats(c_mapinfo=self.c_mapinfo, v_id=self.area_id)
  713. #TODO: why not pass the self.id?
  714. self.read()
  715. # geometry type
  716. self.gtype = libvect.GV_CENTROID
  717. #self.c_pline = ctypes.pointer(libvect.P_line()) if topology else None
  718. def __repr__(self):
  719. return "Centoid(%s)" % ', '.join(['%f' % co for co in self.coords()])
  720. def get_centroid_id(self):
  721. """Return the centroid_id, using the c_mapinfo and an area_id
  722. attributes of the class, and calling the Vect_get_area_centroid
  723. C function, if no centroid_id were found return None"""
  724. centroid_id = libvect.Vect_get_area_centroid(self.c_mapinfo,
  725. self.area_id)
  726. return centroid_id if centroid_id != 0 else None
  727. def get_area_id(self):
  728. """Return the area_id, using the c_mapinfo and an centroid_id
  729. attributes of the class, and calling the Vect_get_centroid_area
  730. C function, if no area_id were found return None"""
  731. area_id = libvect.Vect_get_centroid_area(self.c_mapinfo,
  732. self.id)
  733. return area_id if area_id != 0 else None
  734. class Isle(Geo):
  735. """An Isle is an area contained by another area.
  736. """
  737. def __init__(self, **kargs):
  738. super(Isle, self).__init__(**kargs)
  739. #self.area_id = area_id
  740. def __repr__(self):
  741. return "Isle(%d)" % (self.id)
  742. def boundaries(self):
  743. ilist = Ilist()
  744. libvect.Vect_get_isle_boundaries(self.c_mapinfo, self.id,
  745. ilist.c_ilist)
  746. return ilist
  747. def bbox(self):
  748. bbox = Bbox()
  749. libvect.Vect_get_isle_box(self.c_mapinfo, self.id, bbox.c_bbox)
  750. return bbox
  751. def points(self):
  752. """Return a Line object with the outer ring points"""
  753. line = Line()
  754. libvect.Vect_get_isle_points(self.c_mapinfo, self.id, line.c_points)
  755. return line
  756. def points_geos(self):
  757. """Return a Line object with the outer ring points
  758. """
  759. return libvect.Vect_get_isle_points_geos(self.c_mapinfo, self.id)
  760. def area_id(self):
  761. """Returns area id for isle."""
  762. return libvect.Vect_get_isle_area(self.c_mapinfo, self.id)
  763. def alive(self):
  764. """Check if isle is alive or dead (topology required)"""
  765. return bool(libvect.Vect_isle_alive(self.c_mapinfo, self.id))
  766. def contain_pnt(self, pnt):
  767. """Check if point is in area."""
  768. bbox = self.bbox()
  769. return bool(libvect.Vect_point_in_island(pnt.x, pnt.y,
  770. self.c_mapinfo, self.id,
  771. bbox.c_bbox.contents))
  772. def area(self):
  773. """Return the area value of an Isle"""
  774. border = self.points()
  775. return libgis.G_area_of_polygon(border.c_points.contents.x,
  776. border.c_points.contents.y,
  777. border.c_points.contents.n_points)
  778. def perimeter(self):
  779. """Return the perimeter value of an Isle.
  780. ::
  781. double Vect_area_perimeter()
  782. """
  783. border = self.points()
  784. return libvect.Vect_area_perimeter(border.c_points)
  785. class Isles(object):
  786. def __init__(self, c_mapinfo, area_id):
  787. self.c_mapinfo = c_mapinfo
  788. self.area_id = area_id
  789. self._isles_id = self.get_isles_id()
  790. self._isles = self.get_isles()
  791. def __len__(self):
  792. return libvect.Vect_get_area_num_isles(self.c_mapinfo, self.area_id)
  793. def __repr__(self):
  794. return "Isles(%r)" % self._isles
  795. def __getitem__(self, key):
  796. return self._isles[key]
  797. def get_isles_id(self):
  798. return [libvect.Vect_get_area_isle(self.c_mapinfo, self.area_id, i)
  799. for i in range(self.__len__())]
  800. def get_isles(self):
  801. return [Isle(v_id=isle_id, c_mapinfo=self.c_mapinfo)
  802. for isle_id in self._isles_id]
  803. def select_by_bbox(self, bbox):
  804. """Vect_select_isles_by_box"""
  805. pass
  806. class Area(Geo):
  807. """
  808. 'Vect_build_line_area',
  809. 'Vect_find_area',
  810. 'Vect_get_area_box',
  811. 'Vect_get_area_points_geos',
  812. 'Vect_get_centroid_area',
  813. 'Vect_get_isle_area',
  814. 'Vect_get_line_areas',
  815. 'Vect_get_num_areas',
  816. 'Vect_get_point_in_area',
  817. 'Vect_isle_find_area',
  818. 'Vect_point_in_area',
  819. 'Vect_point_in_area_outer_ring',
  820. 'Vect_read_area_geos',
  821. 'Vect_remove_small_areas',
  822. 'Vect_select_areas_by_box',
  823. 'Vect_select_areas_by_polygon']
  824. """
  825. def __init__(self, boundary=None, centroid=None, isles=[], **kargs):
  826. super(Area, self).__init__(**kargs)
  827. self.boundary = self.points()
  828. self.centroid = self.centroid()
  829. self.isles = self.get_isles()
  830. # geometry type
  831. self.gtype = libvect.GV_AREA
  832. def __repr__(self):
  833. return "Area(%d)" % self.id
  834. def init_from_id(self, area_id=None):
  835. """Return an Area object"""
  836. if area_id is None and self.id is None:
  837. raise ValueError("You need to give or set the area_id")
  838. self.id = area_id if area_id is not None else self.id
  839. # get boundary
  840. self.get_boundary()
  841. # get isles
  842. self.get_isles()
  843. pass
  844. def points(self):
  845. """Return a Line object with the outer ring"""
  846. line = Line()
  847. libvect.Vect_get_area_points(self.c_mapinfo, self.id, line.c_points)
  848. return line
  849. def centroid(self):
  850. centroid_id = libvect.Vect_get_area_centroid(self.c_mapinfo, self.id)
  851. #import pdb; pdb.set_trace()
  852. return Centroid(v_id=centroid_id, c_mapinfo=self.c_mapinfo,
  853. area_id=self.id)
  854. def num_isles(self):
  855. return libvect.Vect_get_area_num_isles(self.c_mapinfo, self.id)
  856. def get_isles(self):
  857. """Instantiate the boundary attribute reading area_id"""
  858. return Isles(self.c_mapinfo, self.id)
  859. def area(self):
  860. """Returns area of area without areas of isles.
  861. double Vect_get_area_area (const struct Map_info *Map, int area)
  862. """
  863. return libvect.Vect_get_area_area(self.c_mapinfo, self.id)
  864. def alive(self):
  865. """Check if area is alive or dead (topology required)
  866. """
  867. return bool(libvect.Vect_area_alive(self.c_mapinfo, self.id))
  868. def bbox(self):
  869. """
  870. Vect_get_area_box
  871. """
  872. bbox = Bbox()
  873. libvect.Vect_get_area_box(self.c_mapinfo, self.id, bbox.c_bbox)
  874. return bbox
  875. def buffer(self):
  876. """Creates buffer around area.
  877. Parameters:
  878. Map vector map
  879. area area id
  880. da distance along major axis
  881. db distance along minor axis
  882. dalpha angle between 0x and major axis
  883. round make corners round
  884. caps add caps at line ends
  885. tol maximum distance between theoretical arc and output segments
  886. [out] oPoints output polygon outer border (ccw order)
  887. [out] inner_count number of holes
  888. [out] iPoints array of output polygon's holes (cw order)
  889. void Vect_area_buffer2(const struct Map_info * Map,
  890. int area,
  891. double da,
  892. double db,
  893. double dalpha,
  894. int round,
  895. int caps,
  896. double tol,
  897. struct line_pnts ** oPoints,
  898. struct line_pnts *** iPoints,
  899. int * inner_count)
  900. """
  901. pass
  902. def boundaries(self):
  903. """Creates list of boundaries for given area.
  904. int Vect_get_area_boundaries(const struct Map_info *Map,
  905. int area, struct ilist *List)
  906. """
  907. ilist = Ilist()
  908. libvect.Vect_get_area_boundaries(self.c_mapinfo, self.id,
  909. ilist.c_ilist)
  910. return ilist
  911. def cats(self):
  912. """Get area categories.
  913. int Vect_get_area_cats (const struct Map_info *Map,
  914. int area, struct line_cats *Cats)
  915. """
  916. return Cats(self.c_mapinfo, self.id)
  917. def get_first_cat(self):
  918. """Find FIRST category of given field and area.
  919. int Vect_get_area_cat(const struct Map_info *Map, int area, int field)
  920. """
  921. pass
  922. def contain_pnt(self, pnt):
  923. """Check if point is in area.
  924. int Vect_point_in_area(double x, double y,
  925. const struct Map_info *Map,
  926. int area, struct bound_box box)
  927. """
  928. bbox = self.bbox()
  929. libvect.Vect_point_in_area(pnt.x, pnt.y, self.c_mapinfo, self.id,
  930. bbox.c_bbox)
  931. return bbox
  932. def perimeter(self):
  933. """Calculate area perimeter.
  934. double Vect_area_perimeter (const struct line_pnts *Points)
  935. """
  936. border = self.points()
  937. return libvect.Vect_area_perimeter(border.c_points)