|
1
|
|
|
2
|
|
|
3
|
|
|
4
|
|
|
5
|
|
|
6
|
|
|
7
|
import os
|
|
8
|
from osgeo import gdal, ogr, osr
|
|
9
|
from processing.core.GeoAlgorithmExecutionException import \
|
|
10
|
GeoAlgorithmExecutionException
|
|
11
|
from processing.tools.raster import *
|
|
12
|
|
|
13
|
raster = gdal.Open(Input_raster)
|
|
14
|
|
|
15
|
rasterBaseName = os.path.splitext(os.path.basename(Input_raster))[0]
|
|
16
|
|
|
17
|
bandCount = raster.RasterCount
|
|
18
|
rasterXSize = raster.RasterXSize
|
|
19
|
rasterYSize = raster.RasterYSize
|
|
20
|
geoTransform = raster.GetGeoTransform()
|
|
21
|
rasterCRS = osr.SpatialReference()
|
|
22
|
rasterCRS.ImportFromWkt(raster.GetProjectionRef())
|
|
23
|
|
|
24
|
vector = ogr.Open(Input_vector, False)
|
|
25
|
layer = vector.GetLayer(0)
|
|
26
|
featureCount = layer.GetFeatureCount()
|
|
27
|
if featureCount == 0:
|
|
28
|
raise GeoAlgorithmExecutionException(
|
|
29
|
'There are no features in input vector.')
|
|
30
|
|
|
31
|
vectorCRS = layer.GetSpatialRef()
|
|
32
|
|
|
33
|
drv = ogr.GetDriverByName('ESRI Shapefile')
|
|
34
|
if drv is None:
|
|
35
|
raise GeoAlgorithmExecutionException(
|
|
36
|
"'ESRI Shapefile' driver is not available.")
|
|
37
|
|
|
38
|
outputDataset = drv.CreateDataSource(Output_layer)
|
|
39
|
if outputDataset is None:
|
|
40
|
raise GeoAlgorithmExecutionException('Creation of output file failed.')
|
|
41
|
|
|
42
|
outputLayer = outputDataset.CreateLayer(
|
|
43
|
str(os.path.splitext(os.path.basename(Output_layer))[0]),
|
|
44
|
vectorCRS, ogr.wkbPoint)
|
|
45
|
if outputLayer is None:
|
|
46
|
raise GeoAlgorithmExecutionException('Layer creation failed.')
|
|
47
|
|
|
48
|
featureDefn = layer.GetLayerDefn()
|
|
49
|
for i in xrange(featureDefn.GetFieldCount()):
|
|
50
|
fieldDefn = featureDefn.GetFieldDefn(i)
|
|
51
|
if outputLayer.CreateField(fieldDefn) != 0:
|
|
52
|
raise GeoAlgorithmExecutionException("Can't create field '%s'."
|
|
53
|
% fieldDefn.GetNameRef())
|
|
54
|
|
|
55
|
columnName = str(rasterBaseName[:8])
|
|
56
|
for i in xrange(bandCount):
|
|
57
|
fieldDefn = ogr.FieldDefn(columnName + '_' + str(i + 1), ogr.OFTReal)
|
|
58
|
fieldDefn.SetWidth(18)
|
|
59
|
fieldDefn.SetPrecision(8)
|
|
60
|
if outputLayer.CreateField(fieldDefn) != 0:
|
|
61
|
raise GeoAlgorithmExecutionException("Can't create field '%s'."
|
|
62
|
% fieldDefn.GetNameRef())
|
|
63
|
|
|
64
|
outputFeature = ogr.Feature(outputLayer.GetLayerDefn())
|
|
65
|
|
|
66
|
current = 0
|
|
67
|
total = bandCount + featureCount * bandCount + featureCount
|
|
68
|
|
|
69
|
layer.ResetReading()
|
|
70
|
feature = layer.GetNextFeature()
|
|
71
|
while feature is not None:
|
|
72
|
current += 1
|
|
73
|
progress.setPercentage(int(current * total))
|
|
74
|
|
|
75
|
outputFeature.SetFrom(feature)
|
|
76
|
if outputLayer.CreateFeature(outputFeature) != 0:
|
|
77
|
raise GeoAlgorithmExecutionException('Failed to add feature.')
|
|
78
|
feature = layer.GetNextFeature()
|
|
79
|
|
|
80
|
vector.Destroy()
|
|
81
|
outputFeature.Destroy()
|
|
82
|
outputDataset.Destroy()
|
|
83
|
|
|
84
|
vector = ogr.Open(Output_layer, True)
|
|
85
|
layer = vector.GetLayer(0)
|
|
86
|
|
|
87
|
if Transform_vector_to_raster_CRS:
|
|
88
|
coordTransform = osr.CoordinateTransformation(vectorCRS, rasterCRS)
|
|
89
|
if coordTransform is None:
|
|
90
|
raise GeoAlgorithmExecutionException(
|
|
91
|
'Error while creating coordinate transformation.')
|
|
92
|
|
|
93
|
for i in xrange(bandCount):
|
|
94
|
current += 1
|
|
95
|
progress.setPercentage(int(current * total))
|
|
96
|
|
|
97
|
rasterBand = raster.GetRasterBand(i + 1)
|
|
98
|
data = rasterBand.ReadAsArray()
|
|
99
|
layer.ResetReading()
|
|
100
|
feature = layer.GetNextFeature()
|
|
101
|
while feature is not None:
|
|
102
|
current += 1
|
|
103
|
progress.setPercentage(int(current * total))
|
|
104
|
|
|
105
|
geometry = feature.GetGeometryRef()
|
|
106
|
x = geometry.GetX()
|
|
107
|
y = geometry.GetY()
|
|
108
|
if Transform_vector_to_raster_CRS:
|
|
109
|
pnt = coordTransform.TransformPoint(x, y, 0)
|
|
110
|
x = pnt[0]
|
|
111
|
y = pnt[1]
|
|
112
|
(rX, rY) = mapToPixel(x, y, geoTransform)
|
|
113
|
if rX >= rasterXSize or rY >= rasterYSize:
|
|
114
|
feature = layer.GetNextFeature()
|
|
115
|
continue
|
|
116
|
value = data[rY, rX]
|
|
117
|
|
|
118
|
feature.SetField(columnName + '_' + str(i + 1), float(value))
|
|
119
|
if layer.SetFeature(feature) != 0:
|
|
120
|
raise GeoAlgorithmExecutionException('Failed to update feature.')
|
|
121
|
|
|
122
|
feature = layer.GetNextFeature()
|
|
123
|
|
|
124
|
rasterBand = None
|
|
125
|
|
|
126
|
raster = None
|
|
127
|
vector.Destroy()
|