|
| 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