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