Skip to content

Commit 5a4eb16

Browse files
committedSep 13, 2013
[processing] add Extract raster values to CSV algorithm
1 parent 0bda3b6 commit 5a4eb16

File tree

1 file changed

+105
-0
lines changed

1 file changed

+105
-0
lines changed
 
Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
##[Example scripts]=group
2+
##Input_raster=raster
3+
##Input_vector=vector
4+
##Transform_vector_to_raster_CRS=boolean
5+
##Output_table=output table
6+
7+
import os
8+
9+
from osgeo import gdal, ogr, osr
10+
11+
from processing.core.TableWriter import TableWriter
12+
from processing.core.GeoAlgorithmExecutionException import GeoAlgorithmExecutionException
13+
14+
from processing.algs import QGISUtils as utils
15+
16+
raster = gdal.Open(Input_raster)
17+
18+
rasterBaseName = os.path.splitext(os.path.basename(Input_raster))[0]
19+
20+
bandCount = raster.RasterCount
21+
rasterXSize = raster.RasterXSize
22+
rasterYSize = raster.RasterYSize
23+
geoTransform = raster.GetGeoTransform()
24+
rasterCRS = osr.SpatialReference()
25+
rasterCRS.ImportFromWkt(raster.GetProjectionRef())
26+
27+
vector = ogr.Open(Input_vector, False)
28+
layer = vector.GetLayer(0)
29+
featureCount = layer.GetFeatureCount()
30+
if featureCount == 0:
31+
raise GeoAlgorithmExecutionException("There are no features in input vector.")
32+
33+
vectorCRS = layer.GetSpatialRef()
34+
35+
columns = []
36+
featureDefn = layer.GetLayerDefn()
37+
for i in xrange(featureDefn.GetFieldCount()):
38+
fieldDefn = featureDefn.GetFieldDefn(i)
39+
columns.append([fieldDefn.GetNameRef()])
40+
41+
layer.ResetReading()
42+
feature = layer.GetNextFeature()
43+
while feature is not None:
44+
for i in xrange(featureDefn.GetFieldCount()):
45+
fieldDefn = featureDefn.GetFieldDefn(i)
46+
if fieldDefn.GetType() == ogr.OFTInteger:
47+
columns[i].append(feature.GetFieldAsInteger(i))
48+
elif fieldDefn.GetType() == ogr.OFTReal:
49+
columns[i].append(feature.GetFieldAsDouble(i))
50+
else:
51+
columns[i].append(feature.GetFieldAsString(i))
52+
feature = layer.GetNextFeature()
53+
54+
current = 0
55+
total = bandCount + featureCount * bandCount
56+
57+
if Transform_vector_to_raster_CRS:
58+
coordTransform = osr.CoordinateTransformation(vectorCRS, rasterCRS)
59+
if coordTransform is None:
60+
raise GeoAlgorithmExecutionException("Error while creating coordinate transformation.")
61+
62+
columnName = rasterBaseName[:8]
63+
for i in xrange(bandCount):
64+
current += 1
65+
progress.setPercentage(int(current * total))
66+
67+
rasterBand = raster.GetRasterBand(i + 1)
68+
data = rasterBand.ReadAsArray()
69+
layer.ResetReading()
70+
feature = layer.GetNextFeature()
71+
col = []
72+
col.append(columnName + '_' + str(i + 1))
73+
while feature is not None:
74+
current += 1
75+
progress.setPercentage(int(current * total))
76+
77+
geometry = feature.GetGeometryRef()
78+
x = geometry.GetX()
79+
y = geometry.GetY()
80+
if Transform_vector_to_raster_CRS:
81+
pnt = coordTransform.TransformPoint(x, y, 0)
82+
x = pnt[0]
83+
y = pnt[1]
84+
rX, rY = utils.mapToPixel(x, y, geoTransform)
85+
if rX > rasterXSize or rY > rasterYSize:
86+
feature = layer.GetNextFeature()
87+
continue
88+
value = data[rY, rX]
89+
col.append(value)
90+
91+
feature = layer.GetNextFeature()
92+
93+
rasterBand = None
94+
columns.append(col)
95+
96+
raster = None
97+
vector.Destroy()
98+
99+
writer = TableWriter(Output_table, "utf-8", [])
100+
row = []
101+
for i in xrange(len(columns[0])):
102+
for col in columns:
103+
row.append(col[i])
104+
writer.addRecord(row)
105+
row[:] = []

0 commit comments

Comments
 (0)
Please sign in to comment.