Skip to content

Commit 1c6aa9b

Browse files
committedNov 19, 2014
[processing] add hypsometric curves calculation algorithm
1 parent 474e667 commit 1c6aa9b

File tree

3 files changed

+247
-1
lines changed

3 files changed

+247
-1
lines changed
 
Lines changed: 206 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,206 @@
1+
# -*- coding: utf-8 -*-
2+
3+
"""
4+
***************************************************************************
5+
HypsometricCurves.py
6+
---------------------
7+
Date : November 2014
8+
Copyright : (C) 2014 by Alexander Bruy
9+
Email : alexander dot bruy at gmail dot com
10+
***************************************************************************
11+
* *
12+
* This program is free software; you can redistribute it and/or modify *
13+
* it under the terms of the GNU General Public License as published by *
14+
* the Free Software Foundation; either version 2 of the License, or *
15+
* (at your option) any later version. *
16+
* *
17+
***************************************************************************
18+
"""
19+
20+
__author__ = 'Alexander Bruy'
21+
__date__ = 'November 2014'
22+
__copyright__ = '(C) 2014, Alexander Bruy'
23+
24+
# This will get replaced with a git SHA1 when you do a git archive
25+
26+
__revision__ = '$Format:%H$'
27+
28+
29+
import os
30+
31+
import numpy
32+
from osgeo import gdal, ogr, osr
33+
34+
from qgis.core import *
35+
36+
from processing.core.GeoAlgorithm import GeoAlgorithm
37+
from processing.core.parameters import ParameterRaster
38+
from processing.core.parameters import ParameterVector
39+
from processing.core.parameters import ParameterNumber
40+
from processing.core.parameters import ParameterBoolean
41+
from processing.core.outputs import OutputDirectory
42+
43+
from processing.tools import raster
44+
from processing.tools import dataobjects
45+
from processing.tools import vector
46+
47+
48+
class HypsometricCurves(GeoAlgorithm):
49+
50+
INPUT_DEM = 'INPUT_DEM'
51+
BOUNDARY_LAYER = 'BOUNDARY_LAYER'
52+
STEP = 'STEP'
53+
USE_PERCENTAGE = 'USE_PERCENTAGE'
54+
OUTPUT_DIRECTORY = 'OUTPUT_DIRECTORY'
55+
56+
def defineCharacteristics(self):
57+
self.name = 'Hypsometric curves'
58+
self.group = 'Raster tools'
59+
60+
self.addParameter(ParameterRaster(self.INPUT_DEM, 'DEM to analyze'))
61+
self.addParameter(ParameterVector(self.BOUNDARY_LAYER, 'Boundary layer',
62+
ParameterVector.VECTOR_TYPE_POLYGON))
63+
self.addParameter(ParameterNumber(self.STEP, 'Step',
64+
0.0, 999999999.999999, 100.0))
65+
self.addParameter(ParameterBoolean(self.USE_PERCENTAGE,
66+
'Use % of area instead of absolute value', False))
67+
68+
self.addOutput(OutputDirectory(self.OUTPUT_DIRECTORY, 'Output directory'))
69+
70+
def processAlgorithm(self, progress):
71+
rasterPath = self.getParameterValue(self.INPUT_DEM)
72+
layer = dataobjects.getObjectFromUri(
73+
self.getParameterValue(self.BOUNDARY_LAYER))
74+
step = self.getParameterValue(self.STEP)
75+
percentage = self.getParameterValue(self.USE_PERCENTAGE)
76+
77+
outputPath = self.getOutputValue(self.OUTPUT_DIRECTORY)
78+
79+
rasterDS = gdal.Open(rasterPath, gdal.GA_ReadOnly)
80+
geoTransform = rasterDS.GetGeoTransform()
81+
rasterBand = rasterDS.GetRasterBand(1)
82+
noData = rasterBand.GetNoDataValue()
83+
84+
cellXSize = abs(geoTransform[1])
85+
cellYSize = abs(geoTransform[5])
86+
rasterXSize = rasterDS.RasterXSize
87+
rasterYSize = rasterDS.RasterYSize
88+
89+
rasterBBox = QgsRectangle(geoTransform[0], geoTransform[3] - cellYSize
90+
* rasterYSize, geoTransform[0] + cellXSize
91+
* rasterXSize, geoTransform[3])
92+
rasterGeom = QgsGeometry.fromRect(rasterBBox)
93+
94+
crs = osr.SpatialReference()
95+
crs.ImportFromProj4(str(layer.crs().toProj4()))
96+
97+
memVectorDriver = ogr.GetDriverByName('Memory')
98+
memRasterDriver = gdal.GetDriverByName('MEM')
99+
100+
features = vector.features(layer)
101+
count = len(features)
102+
total = 100.0 / float(count)
103+
104+
for count, f in enumerate(features):
105+
geom = f.geometry()
106+
intersectedGeom = rasterGeom.intersection(geom)
107+
108+
if intersectedGeom.isGeosEmpty():
109+
progress.setInfo('Feature %d does not intersect raster or '
110+
'entirely located in NODATA area' % f.id())
111+
continue
112+
113+
fName = os.path.join(
114+
outputPath, 'hystogram_%s_%s.csv' % (layer.name(), f.id()))
115+
116+
ogrGeom = ogr.CreateGeometryFromWkt(intersectedGeom.exportToWkt())
117+
bbox = intersectedGeom.boundingBox()
118+
xMin = bbox.xMinimum()
119+
xMax = bbox.xMaximum()
120+
yMin = bbox.yMinimum()
121+
yMax = bbox.yMaximum()
122+
123+
(startColumn, startRow) = raster.mapToPixel(xMin, yMax, geoTransform)
124+
(endColumn, endRow) = raster.mapToPixel(xMax, yMin, geoTransform)
125+
126+
width = endColumn - startColumn
127+
height = endRow - startRow
128+
129+
srcOffset = (startColumn, startRow, width, height)
130+
srcArray = rasterBand.ReadAsArray(*srcOffset)
131+
132+
newGeoTransform = (
133+
geoTransform[0] + srcOffset[0] * geoTransform[1],
134+
geoTransform[1],
135+
0.0,
136+
geoTransform[3] + srcOffset[1] * geoTransform[5],
137+
0.0,
138+
geoTransform[5]
139+
)
140+
141+
memVDS = memVectorDriver.CreateDataSource('out')
142+
memLayer = memVDS.CreateLayer('poly', crs, ogr.wkbPolygon)
143+
144+
ft = ogr.Feature(memLayer.GetLayerDefn())
145+
ft.SetGeometry(ogrGeom)
146+
memLayer.CreateFeature(ft)
147+
ft.Destroy()
148+
149+
rasterizedDS = memRasterDriver.Create('', srcOffset[2],
150+
srcOffset[3], 1, gdal.GDT_Byte)
151+
rasterizedDS.SetGeoTransform(newGeoTransform)
152+
gdal.RasterizeLayer(rasterizedDS, [1], memLayer, burn_values=[1])
153+
rasterizedArray = rasterizedDS.ReadAsArray()
154+
155+
srcArray = numpy.nan_to_num(srcArray)
156+
masked = numpy.ma.MaskedArray(srcArray,
157+
mask=numpy.logical_or(srcArray == noData,
158+
numpy.logical_not(rasterizedArray)))
159+
160+
self.calculateHypsometry(f.id(), fName, progress, masked,
161+
cellXSize, cellYSize, percentage, step)
162+
163+
memVDS = None
164+
rasterizedDS = None
165+
progress.setPercentage(int(count * total))
166+
167+
rasterDS = None
168+
169+
def calculateHypsometry(self, fid, fName, progress, data, pX, pY,
170+
percentage, step):
171+
out = dict()
172+
d = data.compressed()
173+
if d.size == 0:
174+
progress.setInfo('Feature %d does not intersect raster or '
175+
'entirely located in NODATA area' % fid)
176+
return
177+
178+
minValue = d.min()
179+
maxValue = d.max()
180+
startValue = minValue
181+
tmpValue = minValue + step
182+
while startValue < maxValue:
183+
out[tmpValue] = ((startValue <= d) & (d < tmpValue)).sum()
184+
startValue = tmpValue
185+
tmpValue += step
186+
187+
if percentage:
188+
multiplier = 100.0 / float(len(d.flat))
189+
else:
190+
multiplier = pX * pY
191+
192+
for k, v in out.iteritems():
193+
out[k] = v * multiplier
194+
195+
prev = None
196+
for i in sorted(out.items()):
197+
if prev is None:
198+
out[i[0]] = i[1]
199+
else:
200+
out[i[0]] = i[1] + out[prev]
201+
prev = i[0]
202+
203+
writer = vector.TableWriter(fName, 'utf-8', ['Area', 'Elevation'])
204+
for i in sorted(out.items()):
205+
writer.addRecord([i[1], i[0]])
206+
del writer

‎python/plugins/processing/algs/qgis/QGISAlgorithmProvider.py

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -117,6 +117,7 @@
117117
from SetVectorStyle import SetVectorStyle
118118
from SetRasterStyle import SetRasterStyle
119119
from SelectByExpression import SelectByExpression
120+
from HypsometricCurves import HypsometricCurves
120121
# from VectorLayerHistogram import VectorLayerHistogram
121122
# from VectorLayerScatterplot import VectorLayerScatterplot
122123
# from MeanAndStdDevPlot import MeanAndStdDevPlot
@@ -166,7 +167,8 @@ def __init__(self):
166167
RandomPointsPolygonsVariable(),
167168
RandomPointsAlongLines(), PointsToPaths(),
168169
PostGISExecuteSQL(), ImportIntoPostGIS(),
169-
SetVectorStyle(), SetRasterStyle(), SelectByExpression()
170+
SetVectorStyle(), SetRasterStyle(),
171+
SelectByExpression(), HypsometricCurves(),
170172
# ------ raster ------
171173
# CreateConstantRaster(),
172174
# ------ graphics ------
Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
HYPSOMETRIC CURVES
2+
==================
3+
4+
Description
5+
-----------
6+
Calculate hypsometric curves for features of polygon layer and save them as
7+
CSV file for further processing.
8+
9+
Parameters
10+
----------
11+
12+
- ``DEM to analyze [Raster]``: DEM to use for calculating altitudes.
13+
- ``Boundary layer [Vector]``: Polygonal vector layer with boundaries of areas
14+
used to calculate hypsometric curves.
15+
- ``Step [Number]``: Distanse between curves. Default is 100.0
16+
- ``Use % of area instead of absolute value [Boolean]``: Write area percentage
17+
to "Area" field of the CSV file instead of absolute area value.
18+
19+
Outputs
20+
-------
21+
22+
- ``Output directory [Directory]``: Directory where output will be saved. For
23+
each feature from input vector layer CSV file with area and altitude values
24+
will be created.
25+
26+
File name consists of prefix "hystogram_" followed by layer name and feature
27+
ID.
28+
29+
See also
30+
--------
31+
32+
33+
Console usage
34+
-------------
35+
36+
::
37+
38+
processing.runalg('qgis:hypsometriccurves', path_to_DEM, path_to_vector, step, use_percentage, output_path)

0 commit comments

Comments
 (0)
Please sign in to comment.