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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
##[Example scripts]=group
##Input_raster=raster
##Input_vector=vector
##Transform_vector_to_raster_CRS=boolean
##Output_layer=output vector
import os
from osgeo import gdal, ogr, osr
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()
drv = ogr.GetDriverByName('ESRI Shapefile')
if drv is None:
raise GeoAlgorithmExecutionException(
"'ESRI Shapefile' driver is not available.")
outputDataset = drv.CreateDataSource(Output_layer)
if outputDataset is None:
raise GeoAlgorithmExecutionException('Creation of output file failed.')
outputLayer = outputDataset.CreateLayer(
str(os.path.splitext(os.path.basename(Output_layer))[0]),
vectorCRS, ogr.wkbPoint)
if outputLayer is None:
raise GeoAlgorithmExecutionException('Layer creation failed.')
featureDefn = layer.GetLayerDefn()
for i in xrange(featureDefn.GetFieldCount()):
fieldDefn = featureDefn.GetFieldDefn(i)
if outputLayer.CreateField(fieldDefn) != 0:
raise GeoAlgorithmExecutionException("Can't create field '%s'."
% fieldDefn.GetNameRef())
columnName = str(rasterBaseName[:8])
for i in xrange(bandCount):
fieldDefn = ogr.FieldDefn(columnName + '_' + str(i + 1), ogr.OFTReal)
fieldDefn.SetWidth(18)
fieldDefn.SetPrecision(8)
if outputLayer.CreateField(fieldDefn) != 0:
raise GeoAlgorithmExecutionException("Can't create field '%s'."
% fieldDefn.GetNameRef())
outputFeature = ogr.Feature(outputLayer.GetLayerDefn())
current = 0
total = bandCount + featureCount * bandCount + featureCount
layer.ResetReading()
feature = layer.GetNextFeature()
while feature is not None:
current += 1
progress.setPercentage(int(current * total))
outputFeature.SetFrom(feature)
if outputLayer.CreateFeature(outputFeature) != 0:
raise GeoAlgorithmExecutionException('Failed to add feature.')
feature = layer.GetNextFeature()
vector.Destroy()
outputFeature.Destroy()
outputDataset.Destroy()
vector = ogr.Open(Output_layer, True)
layer = vector.GetLayer(0)
if Transform_vector_to_raster_CRS:
coordTransform = osr.CoordinateTransformation(vectorCRS, rasterCRS)
if coordTransform is None:
raise GeoAlgorithmExecutionException(
'Error while creating coordinate transformation.')
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()
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) = mapToPixel(x, y, geoTransform)
if rX >= rasterXSize or rY >= rasterYSize:
feature = layer.GetNextFeature()
continue
value = data[rY, rX]
feature.SetField(columnName + '_' + str(i + 1), float(value))
if layer.SetFeature(feature) != 0:
raise GeoAlgorithmExecutionException('Failed to update feature.')
feature = layer.GetNextFeature()
rasterBand = None
raster = None
vector.Destroy()
|
