root / org.gvsig.toolbox / trunk / org.gvsig.toolbox / org.gvsig.toolbox.algorithm / src / main / java / es / unex / sextante / hydrology / channelNetwork / ChannelNetworkAlgorithm.java @ 59
History | View | Annotate | Download (11.3 KB)
1 | 59 | nbrodin | package es.unex.sextante.hydrology.channelNetwork; |
---|---|---|---|
2 | |||
3 | import java.awt.geom.Point2D; |
||
4 | import java.util.ArrayList; |
||
5 | import java.util.Arrays; |
||
6 | |||
7 | import com.vividsolutions.jts.geom.Coordinate; |
||
8 | import com.vividsolutions.jts.geom.Geometry; |
||
9 | import com.vividsolutions.jts.geom.GeometryFactory; |
||
10 | |||
11 | import es.unex.sextante.additionalInfo.AdditionalInfoNumericalValue; |
||
12 | import es.unex.sextante.core.AnalysisExtent; |
||
13 | import es.unex.sextante.core.GeoAlgorithm; |
||
14 | import es.unex.sextante.core.Sextante; |
||
15 | import es.unex.sextante.dataObjects.IRasterLayer; |
||
16 | import es.unex.sextante.dataObjects.IVectorLayer; |
||
17 | import es.unex.sextante.exceptions.GeoAlgorithmExecutionException; |
||
18 | import es.unex.sextante.exceptions.RepeatedParameterNameException; |
||
19 | import es.unex.sextante.outputs.OutputVectorLayer; |
||
20 | import es.unex.sextante.rasterWrappers.GridCell; |
||
21 | |||
22 | public class ChannelNetworkAlgorithm |
||
23 | extends
|
||
24 | GeoAlgorithm { |
||
25 | |||
26 | private final static int m_iOffsetX[] = { 0, 1, 1, 1, 0, -1, -1, -1 }; |
||
27 | private final static int m_iOffsetY[] = { 1, 1, 0, -1, -1, -1, 0, 1 }; |
||
28 | |||
29 | public static final String METHOD = "METHOD"; |
||
30 | public static final String DEM = "DEM"; |
||
31 | public static final String THRESHOLDLAYER = "THRESHOLDLAYER"; |
||
32 | public static final String THRESHOLD = "THRESHOLD"; |
||
33 | public static final String NETWORK = "NETWORK"; |
||
34 | public static final String NETWORKVECT = "NETWORKVECT"; |
||
35 | |||
36 | private int m_iMethod; |
||
37 | private int m_iNX, m_iNY; |
||
38 | private double m_dThreshold; |
||
39 | |||
40 | private IRasterLayer m_DEM = null; |
||
41 | private IRasterLayer m_Threshold = null; |
||
42 | private IRasterLayer m_Network;
|
||
43 | |||
44 | private ArrayList m_HeadersAndJunctions; |
||
45 | |||
46 | |||
47 | @Override
|
||
48 | public boolean processAlgorithm() throws GeoAlgorithmExecutionException { |
||
49 | |||
50 | m_iMethod = m_Parameters.getParameterValueAsInt(METHOD); |
||
51 | m_DEM = m_Parameters.getParameterValueAsRasterLayer(DEM); |
||
52 | m_Threshold = m_Parameters.getParameterValueAsRasterLayer(THRESHOLDLAYER); |
||
53 | m_dThreshold = m_Parameters.getParameterValueAsDouble(THRESHOLD); |
||
54 | |||
55 | m_Network = getNewRasterLayer(NETWORK, Sextante.getText("Channel_network"), IRasterLayer.RASTER_DATA_TYPE_INT);
|
||
56 | |||
57 | m_Network.assign(0.0);
|
||
58 | |||
59 | final AnalysisExtent extent = m_Network.getWindowGridExtent();
|
||
60 | m_DEM.setWindowExtent(extent); |
||
61 | m_Threshold.setWindowExtent(extent); |
||
62 | |||
63 | m_iNX = m_DEM.getNX(); |
||
64 | m_iNY = m_DEM.getNY(); |
||
65 | |||
66 | calculateChannelNetwork(); |
||
67 | |||
68 | m_Network.setNoDataValue(0.0);
|
||
69 | |||
70 | |||
71 | return !m_Task.isCanceled();
|
||
72 | |||
73 | } |
||
74 | |||
75 | |||
76 | @Override
|
||
77 | public void defineCharacteristics() { |
||
78 | |||
79 | final String sMethod[] = { Sextante.getText("Greater_than"), Sextante.getText("Lower_than") }; |
||
80 | |||
81 | setName(Sextante.getText("Channel_network"));
|
||
82 | setGroup(Sextante.getText("Basic_hydrological_analysis"));
|
||
83 | setUserCanDefineAnalysisExtent(true);
|
||
84 | |||
85 | try {
|
||
86 | m_Parameters.addInputRasterLayer(DEM, Sextante.getText("Elevation"), true); |
||
87 | m_Parameters.addInputRasterLayer(THRESHOLDLAYER, Sextante.getText("Threshold_layer"), true); |
||
88 | m_Parameters.addSelection(METHOD, Sextante.getText("Criteria"), sMethod);
|
||
89 | m_Parameters.addNumericalValue(THRESHOLD, Sextante.getText("Threshold"), 10000, |
||
90 | AdditionalInfoNumericalValue.NUMERICAL_VALUE_DOUBLE); |
||
91 | addOutputRasterLayer(NETWORK, Sextante.getText("Channel_network"));
|
||
92 | addOutputVectorLayer(NETWORKVECT, Sextante.getText("Channel_network"), OutputVectorLayer.SHAPE_TYPE_LINE);
|
||
93 | } |
||
94 | catch (final RepeatedParameterNameException e) { |
||
95 | Sextante.addErrorToLog(e); |
||
96 | } |
||
97 | |||
98 | } |
||
99 | |||
100 | |||
101 | private void calculateChannelNetwork() throws GeoAlgorithmExecutionException { |
||
102 | |||
103 | int i;
|
||
104 | final ArrayList alHeaders = getHeaders(); |
||
105 | |||
106 | if (alHeaders != null) { |
||
107 | m_HeadersAndJunctions = alHeaders; |
||
108 | final Object[] headers = alHeaders.toArray(); |
||
109 | Arrays.sort(headers);
|
||
110 | |||
111 | setProgressText(Sextante.getText("Delineating_channel_network"));
|
||
112 | if (headers.length > 0) { |
||
113 | for (i = 0; (i < headers.length) && setProgress(i, headers.length); i++) { |
||
114 | traceChannel((GridCell) headers[i]); |
||
115 | } |
||
116 | } |
||
117 | |||
118 | if (!m_Task.isCanceled()) {
|
||
119 | calculateOrderAndAddJunctions(); |
||
120 | createVectorLayer(); |
||
121 | } |
||
122 | } |
||
123 | |||
124 | } |
||
125 | |||
126 | |||
127 | private ArrayList getHeaders() { |
||
128 | |||
129 | int iDirection;
|
||
130 | int x, y;
|
||
131 | int ix, iy;
|
||
132 | double dValue;
|
||
133 | double dHeight1, dHeight2;
|
||
134 | boolean bIsHeader;
|
||
135 | final ArrayList headers = new ArrayList(); |
||
136 | |||
137 | for (y = 0; (y < m_iNY) && setProgress(y, m_iNY); y++) { |
||
138 | for (x = 0; x < m_iNX; x++) { |
||
139 | dValue = m_Threshold.getCellValueAsDouble(x, y); |
||
140 | dHeight1 = m_DEM.getCellValueAsDouble(x, y); |
||
141 | if (meetsChannelConditions(dValue)) {
|
||
142 | bIsHeader = true;
|
||
143 | dHeight1 = m_DEM.getCellValueAsDouble(x, y); |
||
144 | if (!m_DEM.isNoDataValue(dHeight1)) {
|
||
145 | for (iDirection = 0; iDirection < 8; iDirection++) { |
||
146 | ix = x + m_iOffsetX[iDirection]; |
||
147 | iy = y + m_iOffsetY[iDirection]; |
||
148 | dValue = m_Threshold.getCellValueAsDouble(ix, iy); |
||
149 | if (meetsChannelConditions(dValue)) {
|
||
150 | dHeight2 = m_DEM.getCellValueAsDouble(ix, iy); |
||
151 | if (dHeight2 >= dHeight1) {
|
||
152 | bIsHeader = false;
|
||
153 | break;
|
||
154 | } |
||
155 | } |
||
156 | } |
||
157 | if (bIsHeader) {
|
||
158 | headers.add(new GridCell(x, y, m_DEM.getCellValueAsDouble(x, y)));
|
||
159 | } |
||
160 | } |
||
161 | } |
||
162 | } |
||
163 | } |
||
164 | |||
165 | if (m_Task.isCanceled()) {
|
||
166 | return null; |
||
167 | } |
||
168 | else {
|
||
169 | return headers;
|
||
170 | } |
||
171 | |||
172 | } |
||
173 | |||
174 | |||
175 | private boolean meetsChannelConditions(final double dValue) { |
||
176 | |||
177 | if (m_iMethod == 0) { |
||
178 | return (dValue > m_dThreshold);
|
||
179 | } |
||
180 | else if (m_iMethod == 1) { |
||
181 | return (dValue < m_dThreshold);
|
||
182 | } |
||
183 | |||
184 | return false; |
||
185 | |||
186 | } |
||
187 | |||
188 | |||
189 | private void traceChannel(final GridCell cell) { |
||
190 | |||
191 | int iDirection;
|
||
192 | int x, y;
|
||
193 | boolean bContinue = true; |
||
194 | |||
195 | x = cell.getX(); |
||
196 | y = cell.getY(); |
||
197 | |||
198 | do {
|
||
199 | m_Network.setCellValue(x, y, -1);
|
||
200 | iDirection = m_DEM.getDirToNextDownslopeCell(x, y); |
||
201 | if (iDirection >= 0) { |
||
202 | x = x + m_iOffsetX[iDirection]; |
||
203 | y = y + m_iOffsetY[iDirection]; |
||
204 | } |
||
205 | else {
|
||
206 | bContinue = false;
|
||
207 | } |
||
208 | } |
||
209 | while (bContinue && !m_Task.isCanceled());
|
||
210 | |||
211 | } |
||
212 | |||
213 | |||
214 | private void calculateOrderAndAddJunctions() { |
||
215 | |||
216 | int x, y;
|
||
217 | |||
218 | setProgressText(Sextante.getText("Calculating_orders"));
|
||
219 | |||
220 | for (y = 0; (y < m_iNY) && setProgress(y, m_iNY); y++) { |
||
221 | for (x = 0; x < m_iNX; x++) { |
||
222 | getStrahlerOrder(x, y); |
||
223 | } |
||
224 | } |
||
225 | |||
226 | } |
||
227 | |||
228 | |||
229 | private int getStrahlerOrder(final int x, |
||
230 | final int y) { |
||
231 | |||
232 | int i;
|
||
233 | int ix, iy;
|
||
234 | int iDirection;
|
||
235 | int iMaxOrder = 1; |
||
236 | int iOrder = 1; |
||
237 | int iMaxOrderCells = 0; |
||
238 | int iUpslopeChannelCells = 0; |
||
239 | |||
240 | if (m_Network.getCellValueAsInt(x, y) == -1) { |
||
241 | ; |
||
242 | for (i = 0; i < 8; i++) { |
||
243 | ix = x + m_iOffsetX[i]; |
||
244 | iy = y + m_iOffsetY[i]; |
||
245 | iDirection = m_DEM.getDirToNextDownslopeCell(ix, iy); |
||
246 | if ((iDirection == (i + 4) % 8) && ((iOrder = m_Network.getCellValueAsInt(ix, iy)) != 0)) { |
||
247 | iUpslopeChannelCells++; |
||
248 | iOrder = m_Network.getCellValueAsInt(ix, iy); |
||
249 | if (iOrder == -1) { |
||
250 | iOrder = getStrahlerOrder(ix, iy); |
||
251 | } |
||
252 | if (iOrder > iMaxOrder) {
|
||
253 | iMaxOrder = iOrder; |
||
254 | iMaxOrderCells = 1;
|
||
255 | } |
||
256 | else if (iOrder == iMaxOrder) { |
||
257 | iMaxOrderCells++; |
||
258 | } |
||
259 | } |
||
260 | } |
||
261 | |||
262 | if (iMaxOrderCells > 1) { |
||
263 | iMaxOrder++; |
||
264 | } |
||
265 | |||
266 | if (iUpslopeChannelCells > 1) { |
||
267 | m_HeadersAndJunctions.add(new GridCell(x, y, m_DEM.getCellValueAsDouble(x, y)));
|
||
268 | } |
||
269 | |||
270 | m_Network.setCellValue(x, y, iMaxOrder); |
||
271 | |||
272 | } |
||
273 | |||
274 | return iMaxOrder;
|
||
275 | |||
276 | } |
||
277 | |||
278 | |||
279 | private void createVectorLayer() throws GeoAlgorithmExecutionException { |
||
280 | |||
281 | int i;
|
||
282 | int x, y;
|
||
283 | int ix, iy;
|
||
284 | int iDirection;
|
||
285 | int iIndexDownslope = -1; |
||
286 | int iOrder;
|
||
287 | boolean bContinue;
|
||
288 | double dLength;
|
||
289 | Point2D pt;
|
||
290 | GridCell cell; |
||
291 | ArrayList coordsList;
|
||
292 | final AnalysisExtent extent = m_DEM.getWindowGridExtent();
|
||
293 | final Object[] values = new Object[4]; |
||
294 | |||
295 | final String sNames[] = { Sextante.getText("ID"), Sextante.getText("Length"), Sextante.getText("Order"), |
||
296 | Sextante.getText("Next") };
|
||
297 | final Class[] types = { Integer.class, Double.class, Integer.class, Integer.class }; |
||
298 | |||
299 | final IVectorLayer network = getNewVectorLayer(NETWORKVECT, Sextante.getText("Channel_network"), |
||
300 | IVectorLayer.SHAPE_TYPE_LINE, types, sNames); |
||
301 | final Object[] headers = m_HeadersAndJunctions.toArray(); |
||
302 | Arrays.sort(headers);
|
||
303 | |||
304 | setProgressText(Sextante.getText("Creating_vector_layer"));
|
||
305 | for (i = headers.length - 1; (i > -1) && setProgress(headers.length - i, headers.length); i--) { |
||
306 | cell = (GridCell) headers[i]; |
||
307 | x = cell.getX(); |
||
308 | y = cell.getY(); |
||
309 | coordsList = new ArrayList(); |
||
310 | pt = extent.getWorldCoordsFromGridCoords(cell); |
||
311 | coordsList.add(new Coordinate(pt.getX(), pt.getY()));
|
||
312 | dLength = 0;
|
||
313 | iOrder = m_Network.getCellValueAsInt(x, y); |
||
314 | bContinue = true;
|
||
315 | do {
|
||
316 | iDirection = m_DEM.getDirToNextDownslopeCell(x, y); |
||
317 | if (iDirection >= 0) { |
||
318 | ix = x + m_iOffsetX[iDirection]; |
||
319 | iy = y + m_iOffsetY[iDirection]; |
||
320 | cell = new GridCell(ix, iy, m_DEM.getCellValueAsDouble(ix, iy));
|
||
321 | pt = extent.getWorldCoordsFromGridCoords(cell); |
||
322 | coordsList.add(new Coordinate(pt.getX(), pt.getY()));
|
||
323 | dLength += m_DEM.getDistToNeighborInDir(iDirection); |
||
324 | iIndexDownslope = m_HeadersAndJunctions.indexOf(cell); |
||
325 | if (iIndexDownslope != -1) { |
||
326 | bContinue = false;
|
||
327 | } |
||
328 | x = ix; |
||
329 | y = iy; |
||
330 | } |
||
331 | else {
|
||
332 | bContinue = false;
|
||
333 | } |
||
334 | |||
335 | } |
||
336 | while (bContinue && !m_Task.isCanceled());
|
||
337 | |||
338 | values[0] = new Integer(i); |
||
339 | values[1] = new Double(dLength); |
||
340 | values[2] = new Integer(iOrder); |
||
341 | values[3] = new Integer(iIndexDownslope); |
||
342 | |||
343 | final Coordinate coords[] = new Coordinate[coordsList.size()]; |
||
344 | for (int j = 0; j < coords.length; j++) { |
||
345 | coords[j] = (Coordinate) coordsList.get(j); |
||
346 | } |
||
347 | final Geometry geom = new GeometryFactory().createLineString(coords); |
||
348 | network.addFeature(geom, values); |
||
349 | } |
||
350 | } |
||
351 | } |