소스 검색

add more capabilities to the function to query raster from vector point

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@55047 15284696-431f-4ddb-bdfa-cd5b030d7da7
Luca Delucchi 12 년 전
부모
커밋
df6d001302
1개의 변경된 파일48개의 추가작업 그리고 6개의 파일을 삭제
  1. 48 6
      lib/python/pygrass/functions.py

+ 48 - 6
lib/python/pygrass/functions.py

@@ -128,20 +128,62 @@ def pixel2coor((col, row), region):
             libraster.Rast_col_to_easting(col, region.c_region))
 
 
-def get_raster_for_points(point, raster):
+def get_raster_for_points(poi_vector, raster, column=None):
     """Query a raster map for each point feature of a vector
 
+    Example ::
+        
+        >>> from grass.pygrass.vector import VectorTopo
+        >>> from grass.pygrass.raster import RasterRow
+        >>> ele = RasterRow('elevation')
+        >>> copy('schools','myschools','vect')
+        >>> sch = VectorTopo('myschools')
+        >>> get_raster_for_points(sch, ele)               # doctest: +ELLIPSIS
+        [(1, 633649.2856743174, 221412.94434781274, 145.06602)...
+        >>> sch.table.columns.add('elevation','double precision')
+        >>> 'elevation' in sch.table.columns
+        True
+        >>> get_raster_for_points(sch, ele, 'elevation')
+        True
+        >>> sch.table.filters.select('NAMESHORT','elevation')
+        Filters('SELECT NAMESHORT, elevation FROM myschools;')
+        >>> cur = sch.table.execute()
+        >>> cur.fetchall()                                # doctest: +ELLIPSIS
+        [(u'SWIFT CREEK', 145.06602), ... (u'9TH GRADE CTR', None)]
+        >>> remove('myschools','vect')
+
+
     Parameters
     -------------
 
     point: point vector object
 
     raster: raster object
+
+    column: column name to update
     """
+    from math import isnan
+    if not column:
+        result = []
     reg = Region()
-    if not point.is_open():
-        point.open()
-    if point.num_primitive_of('point') == 0:
+    if not poi_vector.is_open():
+        poi_vector.open()
+    if not raster.is_open():
+        raster.open()
+    if poi_vector.num_primitive_of('point') == 0:
         raise GrassError(_("Vector doesn't contain points"))
-    values = [raster.get_value(poi.coords, reg) for poi in point.viter('point')]
-    return values
+    for poi in poi_vector.viter('points'):
+        val = raster.get_value(poi, reg)
+        if column:
+            if val is not None and not isnan(val):
+                poi.attrs[column] = val
+        else:
+            if val is not None and not isnan(val):
+                result.append((poi.id, poi.x, poi.y, val))
+            else:
+                result.append((poi.id, poi.x, poi.y, None))
+    if not column:
+        return result
+    else:
+        poi.attrs.commit()
+        return True