1 : /******************************************************************************
2 : * $Id: gdaldem.cpp 19812 2010-06-06 22:37:13Z rouault $
3 : *
4 : * Project: GDAL DEM Utilities
5 : * Purpose:
6 : * Authors: Matthew Perry, perrygeo at gmail.com
7 : * Even Rouault, even dot rouault at mines dash paris dot org
8 : * Howard Butler, hobu.inc at gmail.com
9 : * Chris Yesson, chris dot yesson at ioz dot ac dot uk
10 : *
11 : ******************************************************************************
12 : * Copyright (c) 2006, 2009 Matthew Perry
13 : * Copyright (c) 2009 Even Rouault
14 : * Portions derived from GRASS 4.1 (public domain) See
15 : * http://trac.osgeo.org/gdal/ticket/2975 for more information regarding
16 : * history of this code
17 : *
18 : * Permission is hereby granted, free of charge, to any person obtaining a
19 : * copy of this software and associated documentation files (the "Software"),
20 : * to deal in the Software without restriction, including without limitation
21 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
22 : * and/or sell copies of the Software, and to permit persons to whom the
23 : * Software is furnished to do so, subject to the following conditions:
24 : *
25 : * The above copyright notice and this permission notice shall be included
26 : * in all copies or substantial portions of the Software.
27 : *
28 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
29 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
30 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
31 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
32 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
33 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
34 : * DEALINGS IN THE SOFTWARE.
35 : ****************************************************************************
36 : *
37 : * Slope and aspect calculations based on original method for GRASS GIS 4.1
38 : * by Michael Shapiro, U.S.Army Construction Engineering Research Laboratory
39 : * Olga Waupotitsch, U.S.Army Construction Engineering Research Laboratory
40 : * Marjorie Larson, U.S.Army Construction Engineering Research Laboratory
41 : * as found in GRASS's r.slope.aspect module.
42 : *
43 : * Horn's formula is used to find the first order derivatives in x and y directions
44 : * for slope and aspect calculations: Horn, B. K. P. (1981).
45 : * "Hill Shading and the Reflectance Map", Proceedings of the IEEE, 69(1):14-47.
46 : *
47 : * Other reference :
48 : * Burrough, P.A. and McDonell, R.A., 1998. Principles of Geographical Information
49 : * Systems. p. 190.
50 : *
51 : * Shaded relief based on original method for GRASS GIS 4.1 by Jim Westervelt,
52 : * U.S. Army Construction Engineering Research Laboratory
53 : * as found in GRASS's r.shaded.relief (formerly shade.rel.sh) module.
54 : * ref: "r.mapcalc: An Algebra for GIS and Image Processing",
55 : * by Michael Shapiro and Jim Westervelt, U.S. Army Construction Engineering
56 : * Research Laboratory (March/1991)
57 : *
58 : * Color table of named colors and lookup code derived from src/libes/gis/named_colr.c
59 : * of GRASS 4.1
60 : *
61 : * TRI - Terrain Ruggedness Index is as descibed in Wilson et al (2007)
62 : * this is based on the method of Valentine et al. (2004)
63 : *
64 : * TPI - Topographic Position Index follows the description in Wilson et al (2007), following Weiss (2001)
65 : * The radius is fixed at 1 cell width/height
66 : *
67 : * Roughness - follows the definition in Wilson et al. (2007), which follows Dartnell (2000)
68 : *
69 : * References for TRI/TPI/Roughness:
70 : * Dartnell, P. 2000. Applying Remote Sensing Techniques to map Seafloor
71 : * Geology/Habitat Relationships. Masters Thesis, San Francisco State
72 : * University, pp. 108.
73 : * Valentine, P. C., S. J. Fuller, L. A. Scully. 2004. Terrain Ruggedness
74 : * Analysis and Distribution of Boulder Ridges in the Stellwagen Bank National
75 : * Marine Sanctuary Region (poster). Galway, Ireland: 5th International
76 : * Symposium on Marine Geological and Biological Habitat Mapping (GeoHAB),
77 : * May 2004.
78 : * Weiss, A. D. 2001. Topographic Positions and Landforms Analysis (poster),
79 : * ESRI International User Conference, July 2001. San Diego, CA: ESRI.
80 : * Wilson, M. F. J.; O'Connell, B.; Brown, C.; Guinan, J. C. & Grehan, A. J.
81 : * Multiscale terrain analysis of multibeam bathymetry data for habitat mapping
82 : * on the continental slope Marine Geodesy, 2007, 30, 3-35
83 : ****************************************************************************/
84 :
85 : #include <stdlib.h>
86 : #include <math.h>
87 :
88 : #include "cpl_conv.h"
89 : #include "cpl_string.h"
90 : #include "gdal.h"
91 : #include "gdal_priv.h"
92 :
93 : CPL_CVSID("$Id: gdaldem.cpp 19812 2010-06-06 22:37:13Z rouault $");
94 :
95 : #ifndef M_PI
96 : # define M_PI 3.1415926535897932384626433832795
97 : #endif
98 :
99 : #define INTERPOL(a,b) ((bSrcHasNoData && (EQUAL_TO_NODATA(a, fSrcNoDataValue) || EQUAL_TO_NODATA(b, fSrcNoDataValue))) ? fSrcNoDataValue : 2 * (a) - (b))
100 :
101 : /************************************************************************/
102 : /* Usage() */
103 : /************************************************************************/
104 :
105 0 : static void Usage()
106 :
107 : {
108 : printf( " Usage: \n"
109 : " - To generate a shaded relief map from any GDAL-supported elevation raster : \n\n"
110 : " gdaldem hillshade input_dem output_hillshade \n"
111 : " [-z ZFactor (default=1)] [-s scale* (default=1)] \n"
112 : " [-az Azimuth (default=315)] [-alt Altitude (default=45)]\n"
113 : " [-alg ZevenbergenThorne]\n"
114 : " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
115 : "\n"
116 : " - To generates a slope map from any GDAL-supported elevation raster :\n\n"
117 : " gdaldem slope input_dem output_slope_map \n"
118 : " [-p use percent slope (default=degrees)] [-s scale* (default=1)]\n"
119 : " [-alg ZevenbergenThorne]\n"
120 : " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
121 : "\n"
122 : " - To generate an aspect map from any GDAL-supported elevation raster\n"
123 : " Outputs a 32-bit float tiff with pixel values from 0-360 indicating azimuth :\n\n"
124 : " gdaldem aspect input_dem output_aspect_map \n"
125 : " [-trigonometric] [-zero_for_flat]\n"
126 : " [-alg ZevenbergenThorne]\n"
127 : " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
128 : "\n"
129 : " - To generate a color relief map from any GDAL-supported elevation raster\n"
130 : " gdaldem color-relief input_dem color_text_file output_color_relief_map\n"
131 : " [-alpha] [-exact_color_entry | -nearest_color_entry]\n"
132 : " [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
133 : " where color_text_file contains lines of the format \"elevation_value red green blue\"\n"
134 : "\n"
135 : " - To generate a Terrain Ruggedness Index (TRI) map from any GDAL-supported elevation raster\n"
136 : " gdaldem TRI input_dem output_TRI_map\n"
137 : " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
138 : "\n"
139 : " - To generate a Topographic Position Index (TPI) map from any GDAL-supported elevation raster\n"
140 : " gdaldem TPI input_dem output_TPI_map\n"
141 : " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
142 : "\n"
143 : " - To generate a roughness map from any GDAL-supported elevation raster\n"
144 : " gdaldem roughness input_dem output_roughness_map\n"
145 : " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
146 : "\n"
147 : " Notes : \n"
148 : " Scale is the ratio of vertical units to horizontal\n"
149 0 : " for Feet:Latlong use scale=370400, for Meters:LatLong use scale=111120 \n\n");
150 0 : exit( 1 );
151 : }
152 :
153 : /************************************************************************/
154 : /* ComputeVal() */
155 : /************************************************************************/
156 :
157 : typedef float (*GDALGeneric3x3ProcessingAlg) (float* pafWindow, float fDstNoDataValue, void* pData);
158 :
159 : static float ComputeVal(int bSrcHasNoData, float fSrcNoDataValue,
160 : float* afWin, float fDstNoDataValue,
161 : GDALGeneric3x3ProcessingAlg pfnAlg,
162 : void* pData,
163 95530 : int bComputeAtEdges)
164 : {
165 95530 : if (bSrcHasNoData && EQUAL_TO_NODATA(afWin[4], fSrcNoDataValue))
166 : {
167 0 : return fDstNoDataValue;
168 : }
169 95530 : else if (bSrcHasNoData)
170 : {
171 : int k;
172 859260 : for(k=0;k<9;k++)
173 : {
174 773334 : if (EQUAL_TO_NODATA(afWin[k], fSrcNoDataValue))
175 : {
176 0 : if (bComputeAtEdges)
177 0 : afWin[k] = afWin[4];
178 : else
179 0 : return fDstNoDataValue;
180 : }
181 : }
182 : }
183 :
184 95530 : return pfnAlg(afWin, fDstNoDataValue, pData);
185 : }
186 :
187 : /************************************************************************/
188 : /* GDALGeneric3x3Processing() */
189 : /************************************************************************/
190 :
191 : CPLErr GDALGeneric3x3Processing ( GDALRasterBandH hSrcBand,
192 : GDALRasterBandH hDstBand,
193 : GDALGeneric3x3ProcessingAlg pfnAlg,
194 : void* pData,
195 : int bComputeAtEdges,
196 : GDALProgressFunc pfnProgress,
197 5 : void * pProgressData)
198 : {
199 : CPLErr eErr;
200 : float *pafThreeLineWin; /* 3 line rotating source buffer */
201 : float *pafOutputBuf; /* 1 line destination buffer */
202 : int i, j;
203 :
204 : int bSrcHasNoData, bDstHasNoData;
205 5 : float fSrcNoDataValue = 0.0, fDstNoDataValue = 0.0;
206 :
207 5 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
208 5 : int nYSize = GDALGetRasterBandYSize(hSrcBand);
209 :
210 5 : if (pfnProgress == NULL)
211 0 : pfnProgress = GDALDummyProgress;
212 :
213 : /* -------------------------------------------------------------------- */
214 : /* Initialize progress counter. */
215 : /* -------------------------------------------------------------------- */
216 5 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
217 : {
218 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
219 0 : return CE_Failure;
220 : }
221 :
222 5 : pafOutputBuf = (float *) CPLMalloc(sizeof(float)*nXSize);
223 5 : pafThreeLineWin = (float *) CPLMalloc(3*sizeof(float)*(nXSize+1));
224 :
225 5 : fSrcNoDataValue = (float) GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
226 5 : fDstNoDataValue = (float) GDALGetRasterNoDataValue(hDstBand, &bDstHasNoData);
227 5 : if (!bDstHasNoData)
228 0 : fDstNoDataValue = 0.0;
229 :
230 : // Move a 3x3 pafWindow over each cell
231 : // (where the cell in question is #4)
232 : //
233 : // 0 1 2
234 : // 3 4 5
235 : // 6 7 8
236 :
237 : /* Preload the first 2 lines */
238 15 : for ( i = 0; i < 2 && i < nYSize; i++)
239 : {
240 : GDALRasterIO( hSrcBand,
241 : GF_Read,
242 : 0, i,
243 : nXSize, 1,
244 : pafThreeLineWin + i * nXSize,
245 : nXSize, 1,
246 : GDT_Float32,
247 10 : 0, 0);
248 : }
249 :
250 6 : if (bComputeAtEdges && nXSize >= 2 && nYSize >= 2)
251 : {
252 122 : for (j = 0; j < nXSize; j++)
253 : {
254 : float afWin[9];
255 121 : int jmin = (j == 0) ? j : j - 1;
256 121 : int jmax = (j == nXSize - 1) ? j : j + 1;
257 :
258 121 : afWin[0] = INTERPOL(pafThreeLineWin[jmin], pafThreeLineWin[nXSize + jmin]);
259 121 : afWin[1] = INTERPOL(pafThreeLineWin[j], pafThreeLineWin[nXSize + j]);
260 121 : afWin[2] = INTERPOL(pafThreeLineWin[jmax], pafThreeLineWin[nXSize + jmax]);
261 121 : afWin[3] = pafThreeLineWin[jmin];
262 121 : afWin[4] = pafThreeLineWin[j];
263 121 : afWin[5] = pafThreeLineWin[jmax];
264 121 : afWin[6] = pafThreeLineWin[nXSize + jmin];
265 121 : afWin[7] = pafThreeLineWin[nXSize + j];
266 121 : afWin[8] = pafThreeLineWin[nXSize + jmax];
267 :
268 : pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
269 : afWin, fDstNoDataValue,
270 121 : pfnAlg, pData, bComputeAtEdges);
271 : }
272 : GDALRasterIO(hDstBand, GF_Write,
273 : 0, 0, nXSize, 1,
274 1 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
275 : }
276 : else
277 : {
278 : // Exclude the edges
279 467 : for (j = 0; j < nXSize; j++)
280 : {
281 463 : pafOutputBuf[j] = fDstNoDataValue;
282 : }
283 : GDALRasterIO(hDstBand, GF_Write,
284 : 0, 0, nXSize, 1,
285 4 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
286 :
287 4 : if (nYSize > 1)
288 : {
289 : GDALRasterIO(hDstBand, GF_Write,
290 : 0, nYSize - 1, nXSize, 1,
291 4 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
292 : }
293 : }
294 :
295 5 : int nLine1Off = 0*nXSize;
296 5 : int nLine2Off = 1*nXSize;
297 5 : int nLine3Off = 2*nXSize;
298 :
299 579 : for ( i = 1; i < nYSize-1; i++)
300 : {
301 : /* Read third line of the line buffer */
302 : eErr = GDALRasterIO( hSrcBand,
303 : GF_Read,
304 : 0, i+1,
305 : nXSize, 1,
306 : pafThreeLineWin + nLine3Off,
307 : nXSize, 1,
308 : GDT_Float32,
309 574 : 0, 0);
310 574 : if (eErr != CE_None)
311 0 : goto end;
312 :
313 693 : if (bComputeAtEdges && nXSize >= 2)
314 : {
315 : float afWin[9];
316 :
317 119 : j = 0;
318 119 : afWin[0] = INTERPOL(pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j+1]);
319 119 : afWin[1] = pafThreeLineWin[nLine1Off + j];
320 119 : afWin[2] = pafThreeLineWin[nLine1Off + j+1];
321 119 : afWin[3] = INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j+1]);
322 119 : afWin[4] = pafThreeLineWin[nLine2Off + j];
323 119 : afWin[5] = pafThreeLineWin[nLine2Off + j+1];
324 119 : afWin[6] = INTERPOL(pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j+1]);
325 119 : afWin[7] = pafThreeLineWin[nLine3Off + j];
326 119 : afWin[8] = pafThreeLineWin[nLine3Off + j+1];
327 :
328 : pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
329 : afWin, fDstNoDataValue,
330 119 : pfnAlg, pData, bComputeAtEdges);
331 119 : j = nXSize - 1;
332 :
333 119 : afWin[0] = pafThreeLineWin[nLine1Off + j-1];
334 119 : afWin[1] = pafThreeLineWin[nLine1Off + j];
335 119 : afWin[2] = INTERPOL(pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j-1]);
336 119 : afWin[3] = pafThreeLineWin[nLine2Off + j-1];
337 119 : afWin[4] = pafThreeLineWin[nLine2Off + j];
338 119 : afWin[5] = INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j-1]);
339 119 : afWin[6] = pafThreeLineWin[nLine3Off + j-1];
340 119 : afWin[7] = pafThreeLineWin[nLine3Off + j];
341 119 : afWin[8] = INTERPOL(pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j-1]);
342 :
343 : pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
344 : afWin, fDstNoDataValue,
345 119 : pfnAlg, pData, bComputeAtEdges);
346 : }
347 : else
348 : {
349 : // Exclude the edges
350 455 : pafOutputBuf[0] = fDstNoDataValue;
351 455 : if (nXSize > 1)
352 455 : pafOutputBuf[nXSize - 1] = fDstNoDataValue;
353 : }
354 :
355 66822 : for (j = 1; j < nXSize - 1; j++)
356 : {
357 : float afWin[9];
358 66248 : afWin[0] = pafThreeLineWin[nLine1Off + j-1];
359 66248 : afWin[1] = pafThreeLineWin[nLine1Off + j];
360 66248 : afWin[2] = pafThreeLineWin[nLine1Off + j+1];
361 66248 : afWin[3] = pafThreeLineWin[nLine2Off + j-1];
362 66248 : afWin[4] = pafThreeLineWin[nLine2Off + j];
363 66248 : afWin[5] = pafThreeLineWin[nLine2Off + j+1];
364 66248 : afWin[6] = pafThreeLineWin[nLine3Off + j-1];
365 66248 : afWin[7] = pafThreeLineWin[nLine3Off + j];
366 66248 : afWin[8] = pafThreeLineWin[nLine3Off + j+1];
367 :
368 : pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
369 : afWin, fDstNoDataValue,
370 66248 : pfnAlg, pData, bComputeAtEdges);
371 : }
372 :
373 : /* -----------------------------------------
374 : * Write Line to Raster
375 : */
376 : eErr = GDALRasterIO(hDstBand, GF_Write, 0, i, nXSize, 1,
377 574 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
378 574 : if (eErr != CE_None)
379 0 : goto end;
380 :
381 574 : if( !pfnProgress( 1.0 * (i+1) / nYSize, NULL, pProgressData ) )
382 : {
383 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
384 0 : eErr = CE_Failure;
385 0 : goto end;
386 : }
387 :
388 574 : int nTemp = nLine1Off;
389 574 : nLine1Off = nLine2Off;
390 574 : nLine2Off = nLine3Off;
391 574 : nLine3Off = nTemp;
392 : }
393 :
394 5 : if (bComputeAtEdges && nXSize >= 2 && nYSize >= 2)
395 : {
396 122 : for (j = 0; j < nXSize; j++)
397 : {
398 : float afWin[9];
399 121 : int jmin = (j == 0) ? j : j - 1;
400 121 : int jmax = (j == nXSize - 1) ? j : j + 1;
401 :
402 121 : afWin[0] = pafThreeLineWin[nLine1Off + jmin];
403 121 : afWin[1] = pafThreeLineWin[nLine1Off + j];
404 121 : afWin[2] = pafThreeLineWin[nLine1Off + jmax];
405 121 : afWin[3] = pafThreeLineWin[nLine2Off + jmin];
406 121 : afWin[4] = pafThreeLineWin[nLine2Off + j];
407 121 : afWin[5] = pafThreeLineWin[nLine2Off + jmax];
408 121 : afWin[6] = INTERPOL(pafThreeLineWin[nLine2Off + jmin], pafThreeLineWin[nLine1Off + jmin]);
409 121 : afWin[7] = INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine1Off + j]);
410 121 : afWin[8] = INTERPOL(pafThreeLineWin[nLine2Off + jmax], pafThreeLineWin[nLine1Off + jmax]);
411 :
412 : pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
413 : afWin, fDstNoDataValue,
414 121 : pfnAlg, pData, bComputeAtEdges);
415 : }
416 : GDALRasterIO(hDstBand, GF_Write,
417 : 0, i, nXSize, 1,
418 1 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
419 : }
420 :
421 5 : pfnProgress( 1.0, NULL, pProgressData );
422 5 : eErr = CE_None;
423 :
424 5 : end:
425 5 : CPLFree(pafOutputBuf);
426 5 : CPLFree(pafThreeLineWin);
427 :
428 5 : return eErr;
429 : }
430 :
431 :
432 : /************************************************************************/
433 : /* GDALHillshade() */
434 : /************************************************************************/
435 :
436 : typedef struct
437 : {
438 : double nsres;
439 : double ewres;
440 : double sin_altRadians;
441 : double cos_altRadians_mul_z_scale_factor;
442 : double azRadians;
443 : double square_z_scale_factor;
444 : } GDALHillshadeAlgData;
445 :
446 : /* Unoptimized formulas are :
447 : x = psData->z*((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
448 : (afWin[2] + afWin[5] + afWin[5] + afWin[8])) /
449 : (8.0 * psData->ewres * psData->scale);
450 :
451 : y = psData->z*((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
452 : (afWin[0] + afWin[1] + afWin[1] + afWin[2])) /
453 : (8.0 * psData->nsres * psData->scale);
454 :
455 : slope = M_PI / 2 - atan(sqrt(x*x + y*y));
456 :
457 : aspect = atan2(y,x);
458 :
459 : cang = sin(alt * degreesToRadians) * sin(slope) +
460 : cos(alt * degreesToRadians) * cos(slope) *
461 : cos(az * degreesToRadians - M_PI/2 - aspect);
462 : */
463 :
464 67208 : float GDALHillshadeAlg (float* afWin, float fDstNoDataValue, void* pData)
465 : {
466 67208 : GDALHillshadeAlgData* psData = (GDALHillshadeAlgData*)pData;
467 : double x, y, aspect, xx_plus_yy, cang;
468 :
469 : // First Slope ...
470 : x = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
471 67208 : (afWin[2] + afWin[5] + afWin[5] + afWin[8])) / psData->ewres;
472 :
473 : y = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
474 67208 : (afWin[0] + afWin[1] + afWin[1] + afWin[2])) / psData->nsres;
475 :
476 67208 : xx_plus_yy = x * x + y * y;
477 :
478 : // ... then aspect...
479 67208 : aspect = atan2(y,x);
480 :
481 : // ... then the shade value
482 : cang = (psData->sin_altRadians -
483 : psData->cos_altRadians_mul_z_scale_factor * sqrt(xx_plus_yy) *
484 : sin(aspect - psData->azRadians)) /
485 67208 : sqrt(1 + psData->square_z_scale_factor * xx_plus_yy);
486 :
487 67208 : if (cang <= 0.0)
488 416 : cang = 1.0;
489 : else
490 66792 : cang = 1.0 + (254.0 * cang);
491 :
492 67208 : return cang;
493 : }
494 :
495 0 : float GDALHillshadeZevenbergenThorneAlg (float* afWin, float fDstNoDataValue, void* pData)
496 : {
497 0 : GDALHillshadeAlgData* psData = (GDALHillshadeAlgData*)pData;
498 : double x, y, aspect, xx_plus_yy, cang;
499 :
500 : // First Slope ...
501 0 : x = (afWin[3] - afWin[5]) / psData->ewres;
502 :
503 0 : y = (afWin[7] - afWin[1]) / psData->nsres;
504 :
505 0 : xx_plus_yy = x * x + y * y;
506 :
507 : // ... then aspect...
508 0 : aspect = atan2(y,x);
509 :
510 : // ... then the shade value
511 : cang = (psData->sin_altRadians -
512 : psData->cos_altRadians_mul_z_scale_factor * sqrt(xx_plus_yy) *
513 : sin(aspect - psData->azRadians)) /
514 0 : sqrt(1 + psData->square_z_scale_factor * xx_plus_yy);
515 :
516 0 : if (cang <= 0.0)
517 0 : cang = 1.0;
518 : else
519 0 : cang = 1.0 + (254.0 * cang);
520 :
521 0 : return cang;
522 : }
523 :
524 : void* GDALCreateHillshadeData(double* adfGeoTransform,
525 : double z,
526 : double scale,
527 : double alt,
528 : double az,
529 5 : int bZevenbergenThorne)
530 : {
531 : GDALHillshadeAlgData* pData =
532 5 : (GDALHillshadeAlgData*)CPLMalloc(sizeof(GDALHillshadeAlgData));
533 :
534 5 : const double degreesToRadians = M_PI / 180.0;
535 5 : pData->nsres = adfGeoTransform[5];
536 5 : pData->ewres = adfGeoTransform[1];
537 5 : pData->sin_altRadians = sin(alt * degreesToRadians);
538 5 : pData->azRadians = az * degreesToRadians;
539 5 : double z_scale_factor = z / (((bZevenbergenThorne) ? 2 : 8) * scale);
540 : pData->cos_altRadians_mul_z_scale_factor =
541 5 : cos(alt * degreesToRadians) * z_scale_factor;
542 5 : pData->square_z_scale_factor = z_scale_factor * z_scale_factor;
543 5 : return pData;
544 : }
545 :
546 : /************************************************************************/
547 : /* GDALSlope() */
548 : /************************************************************************/
549 :
550 : typedef struct
551 : {
552 : double nsres;
553 : double ewres;
554 : double scale;
555 : int slopeFormat;
556 : } GDALSlopeAlgData;
557 :
558 14161 : float GDALSlopeHornAlg (float* afWin, float fDstNoDataValue, void* pData)
559 : {
560 14161 : const double radiansToDegrees = 180.0 / M_PI;
561 14161 : GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
562 : double dx, dy, key;
563 :
564 : dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
565 14161 : (afWin[2] + afWin[5] + afWin[5] + afWin[8]))/psData->ewres;
566 :
567 : dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
568 14161 : (afWin[0] + afWin[1] + afWin[1] + afWin[2]))/psData->nsres;
569 :
570 14161 : key = (dx * dx + dy * dy);
571 :
572 14161 : if (psData->slopeFormat == 1)
573 14161 : return atan(sqrt(key) / (8*psData->scale)) * radiansToDegrees;
574 : else
575 0 : return 100*(sqrt(key) / (8*psData->scale));
576 : }
577 :
578 0 : float GDALSlopeZevenbergenThorneAlg (float* afWin, float fDstNoDataValue, void* pData)
579 : {
580 0 : const double radiansToDegrees = 180.0 / M_PI;
581 0 : GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
582 : double dx, dy, key;
583 :
584 0 : dx = (afWin[3] - afWin[5])/psData->ewres;
585 :
586 0 : dy = (afWin[7] - afWin[1])/psData->nsres;
587 :
588 0 : key = (dx * dx + dy * dy);
589 :
590 0 : if (psData->slopeFormat == 1)
591 0 : return atan(sqrt(key) / (2*psData->scale)) * radiansToDegrees;
592 : else
593 0 : return 100*(sqrt(key) / (2*psData->scale));
594 : }
595 :
596 : void* GDALCreateSlopeData(double* adfGeoTransform,
597 : double scale,
598 1 : int slopeFormat)
599 : {
600 : GDALSlopeAlgData* pData =
601 1 : (GDALSlopeAlgData*)CPLMalloc(sizeof(GDALSlopeAlgData));
602 :
603 1 : pData->nsres = adfGeoTransform[5];
604 1 : pData->ewres = adfGeoTransform[1];
605 1 : pData->scale = scale;
606 1 : pData->slopeFormat = slopeFormat;
607 1 : return pData;
608 : }
609 :
610 : /************************************************************************/
611 : /* GDALAspect() */
612 : /************************************************************************/
613 :
614 : typedef struct
615 : {
616 : int bAngleAsAzimuth;
617 : } GDALAspectAlgData;
618 :
619 14161 : float GDALAspectAlg (float* afWin, float fDstNoDataValue, void* pData)
620 : {
621 14161 : const double degreesToRadians = M_PI / 180.0;
622 14161 : GDALAspectAlgData* psData = (GDALAspectAlgData*)pData;
623 : double dx, dy;
624 : float aspect;
625 :
626 : dx = ((afWin[2] + afWin[5] + afWin[5] + afWin[8]) -
627 14161 : (afWin[0] + afWin[3] + afWin[3] + afWin[6]));
628 :
629 : dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
630 14161 : (afWin[0] + afWin[1] + afWin[1] + afWin[2]));
631 :
632 14161 : aspect = atan2(dy,-dx) / degreesToRadians;
633 :
634 18344 : if (dx == 0 && dy == 0)
635 : {
636 : /* Flat area */
637 4183 : aspect = fDstNoDataValue;
638 : }
639 9978 : else if ( psData->bAngleAsAzimuth )
640 : {
641 9978 : if (aspect > 90.0)
642 1243 : aspect = 450.0 - aspect;
643 : else
644 8735 : aspect = 90.0 - aspect;
645 : }
646 : else
647 : {
648 0 : if (aspect < 0)
649 0 : aspect += 360.0;
650 : }
651 :
652 14161 : if (aspect == 360.0)
653 0 : aspect = 0.0;
654 :
655 14161 : return aspect;
656 : }
657 :
658 0 : float GDALAspectZevenbergenThorneAlg (float* afWin, float fDstNoDataValue, void* pData)
659 : {
660 0 : const double degreesToRadians = M_PI / 180.0;
661 0 : GDALAspectAlgData* psData = (GDALAspectAlgData*)pData;
662 : double dx, dy;
663 : float aspect;
664 :
665 0 : dx = (afWin[5] - afWin[3]);
666 :
667 0 : dy = (afWin[7] - afWin[1]);
668 :
669 0 : aspect = atan2(dy,-dx) / degreesToRadians;
670 :
671 0 : if (dx == 0 && dy == 0)
672 : {
673 : /* Flat area */
674 0 : aspect = fDstNoDataValue;
675 : }
676 0 : else if ( psData->bAngleAsAzimuth )
677 : {
678 0 : if (aspect > 90.0)
679 0 : aspect = 450.0 - aspect;
680 : else
681 0 : aspect = 90.0 - aspect;
682 : }
683 : else
684 : {
685 0 : if (aspect < 0)
686 0 : aspect += 360.0;
687 : }
688 :
689 0 : if (aspect == 360.0)
690 0 : aspect = 0.0;
691 :
692 0 : return aspect;
693 : }
694 1 : void* GDALCreateAspectData(int bAngleAsAzimuth)
695 : {
696 : GDALAspectAlgData* pData =
697 1 : (GDALAspectAlgData*)CPLMalloc(sizeof(GDALAspectAlgData));
698 :
699 1 : pData->bAngleAsAzimuth = bAngleAsAzimuth;
700 1 : return pData;
701 : }
702 :
703 : /************************************************************************/
704 : /* GDALColorRelief() */
705 : /************************************************************************/
706 :
707 : typedef struct
708 : {
709 : double dfVal;
710 : int nR;
711 : int nG;
712 : int nB;
713 : int nA;
714 : } ColorAssociation;
715 :
716 77 : static int GDALColorReliefSortColors(const void* pA, const void* pB)
717 : {
718 77 : ColorAssociation* pC1 = (ColorAssociation*)pA;
719 77 : ColorAssociation* pC2 = (ColorAssociation*)pB;
720 : return (pC1->dfVal < pC2->dfVal) ? -1 :
721 77 : (pC1->dfVal == pC2->dfVal) ? 0 : 1;
722 : }
723 :
724 : typedef enum
725 : {
726 : COLOR_SELECTION_INTERPOLATE,
727 : COLOR_SELECTION_NEAREST_ENTRY,
728 : COLOR_SELECTION_EXACT_ENTRY
729 : } ColorSelectionMode;
730 :
731 : static int GDALColorReliefGetRGBA (ColorAssociation* pasColorAssociation,
732 : int nColorAssociation,
733 : double dfVal,
734 : ColorSelectionMode eColorSelectionMode,
735 : int* pnR,
736 : int* pnG,
737 : int* pnB,
738 131769 : int* pnA)
739 : {
740 : int i;
741 131769 : int lower = 0;
742 131769 : int upper = nColorAssociation - 1;
743 : int mid;
744 :
745 : /* Find the index of the first element in the LUT input array that */
746 : /* is not smaller than the dfVal value. */
747 288783 : while(TRUE)
748 : {
749 420552 : mid = (lower + upper) / 2;
750 420552 : if (upper - lower <= 1)
751 : {
752 131769 : if (dfVal < pasColorAssociation[lower].dfVal)
753 0 : i = lower;
754 131769 : else if (dfVal < pasColorAssociation[upper].dfVal)
755 89685 : i = upper;
756 : else
757 42084 : i = upper + 1;
758 : break;
759 : }
760 288783 : else if (pasColorAssociation[mid].dfVal >= dfVal)
761 : {
762 173367 : upper = mid;
763 : }
764 : else
765 : {
766 115416 : lower = mid;
767 : }
768 : }
769 :
770 131769 : if (i == 0)
771 : {
772 0 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
773 : pasColorAssociation[0].dfVal != dfVal)
774 : {
775 0 : *pnR = 0;
776 0 : *pnG = 0;
777 0 : *pnB = 0;
778 0 : *pnA = 0;
779 0 : return FALSE;
780 : }
781 : else
782 : {
783 0 : *pnR = pasColorAssociation[0].nR;
784 0 : *pnG = pasColorAssociation[0].nG;
785 0 : *pnB = pasColorAssociation[0].nB;
786 0 : *pnA = pasColorAssociation[0].nA;
787 0 : return TRUE;
788 : }
789 : }
790 131769 : else if (i == nColorAssociation)
791 : {
792 0 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
793 : pasColorAssociation[i-1].dfVal != dfVal)
794 : {
795 0 : *pnR = 0;
796 0 : *pnG = 0;
797 0 : *pnB = 0;
798 0 : *pnA = 0;
799 0 : return FALSE;
800 : }
801 : else
802 : {
803 0 : *pnR = pasColorAssociation[i-1].nR;
804 0 : *pnG = pasColorAssociation[i-1].nG;
805 0 : *pnB = pasColorAssociation[i-1].nB;
806 0 : *pnA = pasColorAssociation[i-1].nA;
807 0 : return TRUE;
808 : }
809 : }
810 : else
811 : {
812 131769 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
813 : pasColorAssociation[i-1].dfVal != dfVal)
814 : {
815 0 : *pnR = 0;
816 0 : *pnG = 0;
817 0 : *pnB = 0;
818 0 : *pnA = 0;
819 0 : return FALSE;
820 : }
821 :
822 131769 : if (eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY &&
823 : pasColorAssociation[i-1].dfVal != dfVal)
824 : {
825 : int index;
826 9965 : if (dfVal - pasColorAssociation[i-1].dfVal <
827 : pasColorAssociation[i].dfVal - dfVal)
828 6716 : index = i -1;
829 : else
830 3249 : index = i;
831 :
832 9965 : *pnR = pasColorAssociation[index].nR;
833 9965 : *pnG = pasColorAssociation[index].nG;
834 9965 : *pnB = pasColorAssociation[index].nB;
835 9965 : *pnA = pasColorAssociation[index].nA;
836 9965 : return TRUE;
837 : }
838 :
839 : double dfRatio = (dfVal - pasColorAssociation[i-1].dfVal) /
840 121804 : (pasColorAssociation[i].dfVal - pasColorAssociation[i-1].dfVal);
841 : *pnR = (int)(0.45 + pasColorAssociation[i-1].nR + dfRatio *
842 121804 : (pasColorAssociation[i].nR - pasColorAssociation[i-1].nR));
843 121804 : if (*pnR < 0) *pnR = 0;
844 121804 : else if (*pnR > 255) *pnR = 255;
845 : *pnG = (int)(0.45 + pasColorAssociation[i-1].nG + dfRatio *
846 121804 : (pasColorAssociation[i].nG - pasColorAssociation[i-1].nG));
847 121804 : if (*pnG < 0) *pnG = 0;
848 121804 : else if (*pnG > 255) *pnG = 255;
849 : *pnB = (int)(0.45 + pasColorAssociation[i-1].nB + dfRatio *
850 121804 : (pasColorAssociation[i].nB - pasColorAssociation[i-1].nB));
851 121804 : if (*pnB < 0) *pnB = 0;
852 121804 : else if (*pnB > 255) *pnB = 255;
853 : *pnA = (int)(0.45 + pasColorAssociation[i-1].nA + dfRatio *
854 121804 : (pasColorAssociation[i].nA - pasColorAssociation[i-1].nA));
855 121804 : if (*pnA < 0) *pnA = 0;
856 121804 : else if (*pnA > 255) *pnA = 255;
857 :
858 121804 : return TRUE;
859 : }
860 : }
861 :
862 : /* dfPct : percentage between 0 and 1 */
863 : static double GDALColorReliefGetAbsoluteValFromPct(GDALRasterBandH hSrcBand,
864 0 : double dfPct)
865 : {
866 : double dfMin, dfMax;
867 : int bSuccessMin, bSuccessMax;
868 0 : dfMin = GDALGetRasterMinimum(hSrcBand, &bSuccessMin);
869 0 : dfMax = GDALGetRasterMaximum(hSrcBand, &bSuccessMax);
870 0 : if (!bSuccessMin || !bSuccessMax)
871 : {
872 : double dfMean, dfStdDev;
873 0 : fprintf(stderr, "Computing source raster statistics...\n");
874 : GDALComputeRasterStatistics(hSrcBand, FALSE, &dfMin, &dfMax,
875 0 : &dfMean, &dfStdDev, NULL, NULL);
876 : }
877 0 : return dfMin + dfPct * (dfMax - dfMin);
878 : }
879 :
880 : typedef struct
881 : {
882 : const char *name;
883 : float r, g, b;
884 : } NamedColor;
885 :
886 : static const NamedColor namedColors[] = {
887 : { "white", 1.00, 1.00, 1.00 },
888 : { "black", 0.00, 0.00, 0.00 },
889 : { "red", 1.00, 0.00, 0.00 },
890 : { "green", 0.00, 1.00, 0.00 },
891 : { "blue", 0.00, 0.00, 1.00 },
892 : { "yellow", 1.00, 1.00, 0.00 },
893 : { "magenta",1.00, 0.00, 1.00 },
894 : { "cyan", 0.00, 1.00, 1.00 },
895 : { "aqua", 0.00, 0.75, 0.75 },
896 : { "grey", 0.75, 0.75, 0.75 },
897 : { "gray", 0.75, 0.75, 0.75 },
898 : { "orange", 1.00, 0.50, 0.00 },
899 : { "brown", 0.75, 0.50, 0.25 },
900 : { "purple", 0.50, 0.00, 1.00 },
901 : { "violet", 0.50, 0.00, 1.00 },
902 : { "indigo", 0.00, 0.50, 1.00 },
903 : };
904 :
905 : static
906 0 : int GDALColorReliefFindNamedColor(const char *pszColorName, int *pnR, int *pnG, int *pnB)
907 : {
908 : unsigned int i;
909 :
910 0 : *pnR = *pnG = *pnB = 0;
911 0 : for (i = 0; i < sizeof(namedColors) / sizeof(namedColors[0]); i++)
912 : {
913 0 : if (EQUAL(pszColorName, namedColors[i].name))
914 : {
915 0 : *pnR = (int)(255. * namedColors[i].r);
916 0 : *pnG = (int)(255. * namedColors[i].g);
917 0 : *pnB = (int)(255. * namedColors[i].b);
918 0 : return TRUE;
919 : }
920 : }
921 0 : return FALSE;
922 : }
923 :
924 : static
925 : ColorAssociation* GDALColorReliefParseColorFile(GDALRasterBandH hSrcBand,
926 : const char* pszColorFilename,
927 7 : int* pnColors)
928 : {
929 7 : FILE* fpColorFile = VSIFOpenL(pszColorFilename, "rt");
930 7 : if (fpColorFile == NULL)
931 : {
932 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszColorFilename);
933 0 : *pnColors = 0;
934 0 : return NULL;
935 : }
936 :
937 7 : ColorAssociation* pasColorAssociation = NULL;
938 7 : int nColorAssociation = 0;
939 :
940 7 : int bSrcHasNoData = FALSE;
941 7 : double dfSrcNoDataValue = GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
942 :
943 : const char* pszLine;
944 63 : while ((pszLine = CPLReadLineL(fpColorFile)) != NULL)
945 : {
946 : char** papszFields = CSLTokenizeStringComplex(pszLine, " ,\t:",
947 49 : FALSE, FALSE );
948 : /* Skip comment lines */
949 49 : int nTokens = CSLCount(papszFields);
950 49 : if (nTokens >= 2 &&
951 : papszFields[0][0] != '#' &&
952 : papszFields[0][0] != '/')
953 : {
954 : pasColorAssociation =
955 : (ColorAssociation*)CPLRealloc(pasColorAssociation,
956 49 : (nColorAssociation + 1) * sizeof(ColorAssociation));
957 49 : if (EQUAL(papszFields[0], "nv") && bSrcHasNoData)
958 0 : pasColorAssociation[nColorAssociation].dfVal = dfSrcNoDataValue;
959 49 : else if (strlen(papszFields[0]) > 1 && papszFields[0][strlen(papszFields[0])-1] == '%')
960 : {
961 0 : double dfPct = atof(papszFields[0]) / 100.;
962 0 : if (dfPct < 0.0 || dfPct > 1.0)
963 : {
964 : CPLError(CE_Failure, CPLE_AppDefined,
965 0 : "Wrong value for a percentage : %s", papszFields[0]);
966 0 : CSLDestroy(papszFields);
967 0 : VSIFCloseL(fpColorFile);
968 0 : CPLFree(pasColorAssociation);
969 0 : *pnColors = 0;
970 0 : return NULL;
971 : }
972 : pasColorAssociation[nColorAssociation].dfVal =
973 0 : GDALColorReliefGetAbsoluteValFromPct(hSrcBand, dfPct);
974 : }
975 : else
976 49 : pasColorAssociation[nColorAssociation].dfVal = atof(papszFields[0]);
977 :
978 49 : if (nTokens >= 4)
979 : {
980 49 : pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]);
981 49 : pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]);
982 49 : pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]);
983 : pasColorAssociation[nColorAssociation].nA =
984 49 : (CSLCount(papszFields) >= 5 ) ? atoi(papszFields[4]) : 255;
985 : }
986 : else
987 : {
988 : int nR, nG, nB;
989 0 : if (!GDALColorReliefFindNamedColor(papszFields[1], &nR, &nG, &nB))
990 : {
991 : CPLError(CE_Failure, CPLE_AppDefined,
992 0 : "Unknown color : %s", papszFields[1]);
993 0 : CSLDestroy(papszFields);
994 0 : VSIFCloseL(fpColorFile);
995 0 : CPLFree(pasColorAssociation);
996 0 : *pnColors = 0;
997 0 : return NULL;
998 : }
999 0 : pasColorAssociation[nColorAssociation].nR = nR;
1000 0 : pasColorAssociation[nColorAssociation].nG = nG;
1001 0 : pasColorAssociation[nColorAssociation].nB = nB;
1002 : pasColorAssociation[nColorAssociation].nA =
1003 0 : (CSLCount(papszFields) >= 3 ) ? atoi(papszFields[2]) : 255;
1004 : }
1005 :
1006 49 : nColorAssociation ++;
1007 : }
1008 49 : CSLDestroy(papszFields);
1009 : }
1010 7 : VSIFCloseL(fpColorFile);
1011 :
1012 7 : if (nColorAssociation == 0)
1013 : {
1014 : CPLError(CE_Failure, CPLE_AppDefined,
1015 0 : "No color association found in %s", pszColorFilename);
1016 0 : *pnColors = 0;
1017 0 : return NULL;
1018 : }
1019 :
1020 : qsort(pasColorAssociation, nColorAssociation,
1021 7 : sizeof(ColorAssociation), GDALColorReliefSortColors);
1022 :
1023 7 : *pnColors = nColorAssociation;
1024 7 : return pasColorAssociation;
1025 : }
1026 :
1027 : static
1028 : GByte* GDALColorReliefPrecompute(GDALRasterBandH hSrcBand,
1029 : ColorAssociation* pasColorAssociation,
1030 : int nColorAssociation,
1031 : ColorSelectionMode eColorSelectionMode,
1032 5 : int* pnIndexOffset)
1033 : {
1034 5 : GDALDataType eDT = GDALGetRasterDataType(hSrcBand);
1035 5 : GByte* pabyPrecomputed = NULL;
1036 5 : int nIndexOffset = (eDT == GDT_Int16) ? 32768 : 0;
1037 5 : *pnIndexOffset = nIndexOffset;
1038 5 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
1039 5 : int nYSize = GDALGetRasterBandXSize(hSrcBand);
1040 5 : if (eDT == GDT_Byte ||
1041 : ((eDT == GDT_Int16 || eDT == GDT_UInt16) && nXSize * nYSize > 65536))
1042 : {
1043 0 : int iMax = (eDT == GDT_Byte) ? 256: 65536;
1044 0 : pabyPrecomputed = (GByte*) VSIMalloc(4 * iMax);
1045 0 : if (pabyPrecomputed)
1046 : {
1047 : int i;
1048 0 : for(i=0;i<iMax;i++)
1049 : {
1050 : int nR, nG, nB, nA;
1051 : GDALColorReliefGetRGBA (pasColorAssociation,
1052 : nColorAssociation,
1053 : i - nIndexOffset,
1054 : eColorSelectionMode,
1055 0 : &nR, &nG, &nB, &nA);
1056 0 : pabyPrecomputed[4 * i] = (GByte) nR;
1057 0 : pabyPrecomputed[4 * i + 1] = (GByte) nG;
1058 0 : pabyPrecomputed[4 * i + 2] = (GByte) nB;
1059 0 : pabyPrecomputed[4 * i + 3] = (GByte) nA;
1060 : }
1061 : }
1062 : }
1063 5 : return pabyPrecomputed;
1064 : }
1065 :
1066 : /************************************************************************/
1067 : /* ==================================================================== */
1068 : /* GDALColorReliefDataset */
1069 : /* ==================================================================== */
1070 : /************************************************************************/
1071 :
1072 : class GDALColorReliefRasterBand;
1073 :
1074 : class GDALColorReliefDataset : public GDALDataset
1075 : {
1076 : friend class GDALColorReliefRasterBand;
1077 :
1078 : GDALDatasetH hSrcDS;
1079 : GDALRasterBandH hSrcBand;
1080 : int nColorAssociation;
1081 : ColorAssociation* pasColorAssociation;
1082 : ColorSelectionMode eColorSelectionMode;
1083 : GByte* pabyPrecomputed;
1084 : int nIndexOffset;
1085 : float* pafSourceBuf;
1086 : int* panSourceBuf;
1087 : int nCurBlockXOff;
1088 : int nCurBlockYOff;
1089 :
1090 : public:
1091 : GDALColorReliefDataset(GDALDatasetH hSrcDS,
1092 : GDALRasterBandH hSrcBand,
1093 : const char* pszColorFilename,
1094 : ColorSelectionMode eColorSelectionMode,
1095 : int bAlpha);
1096 : ~GDALColorReliefDataset();
1097 :
1098 : CPLErr GetGeoTransform( double * padfGeoTransform );
1099 : const char *GetProjectionRef();
1100 : };
1101 :
1102 : /************************************************************************/
1103 : /* ==================================================================== */
1104 : /* GDALColorReliefRasterBand */
1105 : /* ==================================================================== */
1106 : /************************************************************************/
1107 :
1108 : class GDALColorReliefRasterBand : public GDALRasterBand
1109 6 : {
1110 : friend class GDALColorReliefDataset;
1111 :
1112 :
1113 : public:
1114 : GDALColorReliefRasterBand( GDALColorReliefDataset *, int );
1115 :
1116 : virtual CPLErr IReadBlock( int, int, void * );
1117 : virtual GDALColorInterp GetColorInterpretation();
1118 : };
1119 :
1120 : GDALColorReliefDataset::GDALColorReliefDataset(
1121 : GDALDatasetH hSrcDS,
1122 : GDALRasterBandH hSrcBand,
1123 : const char* pszColorFilename,
1124 : ColorSelectionMode eColorSelectionMode,
1125 2 : int bAlpha)
1126 : {
1127 2 : this->hSrcDS = hSrcDS;
1128 2 : this->hSrcBand = hSrcBand;
1129 2 : nColorAssociation = 0;
1130 : pasColorAssociation =
1131 : GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
1132 2 : &nColorAssociation);
1133 2 : this->eColorSelectionMode = eColorSelectionMode;
1134 :
1135 2 : nRasterXSize = GDALGetRasterXSize(hSrcDS);
1136 2 : nRasterYSize = GDALGetRasterYSize(hSrcDS);
1137 :
1138 : int nBlockXSize, nBlockYSize;
1139 2 : GDALGetBlockSize( hSrcBand, &nBlockXSize, &nBlockYSize);
1140 :
1141 2 : nIndexOffset = 0;
1142 : pabyPrecomputed =
1143 : GDALColorReliefPrecompute(hSrcBand,
1144 : pasColorAssociation,
1145 : nColorAssociation,
1146 : eColorSelectionMode,
1147 2 : &nIndexOffset);
1148 :
1149 : int i;
1150 8 : for(i=0;i<((bAlpha) ? 4 : 3);i++)
1151 : {
1152 6 : SetBand(i + 1, new GDALColorReliefRasterBand(this, i+1));
1153 : }
1154 :
1155 2 : pafSourceBuf = NULL;
1156 2 : panSourceBuf = NULL;
1157 2 : if (pabyPrecomputed)
1158 0 : panSourceBuf = (int *) CPLMalloc(sizeof(int)*nBlockXSize*nBlockYSize);
1159 : else
1160 2 : pafSourceBuf = (float *) CPLMalloc(sizeof(float)*nBlockXSize*nBlockYSize);
1161 2 : nCurBlockXOff = -1;
1162 2 : nCurBlockYOff = -1;
1163 2 : }
1164 :
1165 2 : GDALColorReliefDataset::~GDALColorReliefDataset()
1166 : {
1167 2 : CPLFree(pasColorAssociation);
1168 2 : CPLFree(pabyPrecomputed);
1169 2 : CPLFree(panSourceBuf);
1170 2 : CPLFree(pafSourceBuf);
1171 2 : }
1172 :
1173 2 : CPLErr GDALColorReliefDataset::GetGeoTransform( double * padfGeoTransform )
1174 : {
1175 2 : return GDALGetGeoTransform(hSrcDS, padfGeoTransform);
1176 : }
1177 :
1178 2 : const char *GDALColorReliefDataset::GetProjectionRef()
1179 : {
1180 2 : return GDALGetProjectionRef(hSrcDS);
1181 : }
1182 :
1183 : GDALColorReliefRasterBand::GDALColorReliefRasterBand(
1184 6 : GDALColorReliefDataset * poDS, int nBand)
1185 : {
1186 6 : this->poDS = poDS;
1187 6 : this->nBand = nBand;
1188 6 : eDataType = GDT_Byte;
1189 6 : GDALGetBlockSize( poDS->hSrcBand, &nBlockXSize, &nBlockYSize);
1190 6 : }
1191 :
1192 : CPLErr GDALColorReliefRasterBand::IReadBlock( int nBlockXOff,
1193 : int nBlockYOff,
1194 387 : void *pImage )
1195 : {
1196 387 : GDALColorReliefDataset * poGDS = (GDALColorReliefDataset *) poDS;
1197 : int nReqXSize, nReqYSize;
1198 :
1199 387 : if ((nBlockXOff + 1) * nBlockXSize >= nRasterXSize)
1200 27 : nReqXSize = nRasterXSize - nBlockXOff * nBlockXSize;
1201 : else
1202 360 : nReqXSize = nBlockXSize;
1203 :
1204 387 : if ((nBlockYOff + 1) * nBlockYSize >= nRasterYSize)
1205 366 : nReqYSize = nRasterYSize - nBlockYOff * nBlockYSize;
1206 : else
1207 21 : nReqYSize = nBlockYSize;
1208 :
1209 387 : int nCount = nReqXSize * nReqYSize;
1210 387 : if ( poGDS->nCurBlockXOff != nBlockXOff ||
1211 : poGDS->nCurBlockYOff != nBlockYOff )
1212 : {
1213 371 : poGDS->nCurBlockXOff = nBlockXOff;
1214 371 : poGDS->nCurBlockYOff = nBlockYOff;
1215 :
1216 : CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1217 : GF_Read,
1218 : nBlockXOff * nBlockXSize,
1219 : nBlockYOff * nBlockYSize,
1220 : nReqXSize, nReqYSize,
1221 : (poGDS->panSourceBuf) ?
1222 : (void*) poGDS->panSourceBuf :
1223 : (void* )poGDS->pafSourceBuf,
1224 : nReqXSize, nReqYSize,
1225 : (poGDS->panSourceBuf) ? GDT_Int32 : GDT_Float32,
1226 371 : 0, 0);
1227 371 : if (eErr != CE_None)
1228 : {
1229 0 : memset(pImage, 0, nCount);
1230 0 : return eErr;
1231 : }
1232 : }
1233 :
1234 : int j;
1235 387 : if (poGDS->panSourceBuf)
1236 : {
1237 0 : for ( j = 0; j < nCount; j++)
1238 : {
1239 0 : int nIndex = poGDS->panSourceBuf[j] + poGDS->nIndexOffset;
1240 0 : ((GByte*)pImage)[j] = poGDS->pabyPrecomputed[4*nIndex + nBand-1];
1241 : }
1242 : }
1243 : else
1244 : {
1245 : int anComponents[4];
1246 88233 : for ( j = 0; j < nCount; j++)
1247 : {
1248 : GDALColorReliefGetRGBA (poGDS->pasColorAssociation,
1249 : poGDS->nColorAssociation,
1250 : poGDS->pafSourceBuf[j],
1251 : poGDS->eColorSelectionMode,
1252 : &anComponents[0],
1253 : &anComponents[1],
1254 : &anComponents[2],
1255 87846 : &anComponents[3]);
1256 87846 : ((GByte*)pImage)[j] = (GByte) anComponents[nBand-1];
1257 : }
1258 : }
1259 :
1260 387 : return CE_None;
1261 : }
1262 :
1263 12 : GDALColorInterp GDALColorReliefRasterBand::GetColorInterpretation()
1264 : {
1265 12 : return (GDALColorInterp)(GCI_RedBand + nBand - 1);
1266 : }
1267 :
1268 :
1269 : CPLErr GDALColorRelief (GDALRasterBandH hSrcBand,
1270 : GDALRasterBandH hDstBand1,
1271 : GDALRasterBandH hDstBand2,
1272 : GDALRasterBandH hDstBand3,
1273 : GDALRasterBandH hDstBand4,
1274 : const char* pszColorFilename,
1275 : ColorSelectionMode eColorSelectionMode,
1276 : GDALProgressFunc pfnProgress,
1277 3 : void * pProgressData)
1278 : {
1279 : CPLErr eErr;
1280 :
1281 3 : if (hSrcBand == NULL || hDstBand1 == NULL || hDstBand2 == NULL ||
1282 : hDstBand3 == NULL)
1283 0 : return CE_Failure;
1284 :
1285 3 : int nColorAssociation = 0;
1286 : ColorAssociation* pasColorAssociation =
1287 : GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
1288 3 : &nColorAssociation);
1289 3 : if (pasColorAssociation == NULL)
1290 0 : return CE_Failure;
1291 :
1292 3 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
1293 3 : int nYSize = GDALGetRasterBandYSize(hSrcBand);
1294 :
1295 3 : if (pfnProgress == NULL)
1296 0 : pfnProgress = GDALDummyProgress;
1297 :
1298 3 : int nR = 0, nG = 0, nB = 0, nA = 0;
1299 :
1300 : /* -------------------------------------------------------------------- */
1301 : /* Precompute the map from values to RGBA quadruplets */
1302 : /* for GDT_Byte, GDT_Int16 or GDT_UInt16 */
1303 : /* -------------------------------------------------------------------- */
1304 3 : int nIndexOffset = 0;
1305 : GByte* pabyPrecomputed =
1306 : GDALColorReliefPrecompute(hSrcBand,
1307 : pasColorAssociation,
1308 : nColorAssociation,
1309 : eColorSelectionMode,
1310 3 : &nIndexOffset);
1311 :
1312 : /* -------------------------------------------------------------------- */
1313 : /* Initialize progress counter. */
1314 : /* -------------------------------------------------------------------- */
1315 :
1316 3 : float* pafSourceBuf = NULL;
1317 3 : int* panSourceBuf = NULL;
1318 3 : if (pabyPrecomputed)
1319 0 : panSourceBuf = (int *) CPLMalloc(sizeof(int)*nXSize);
1320 : else
1321 3 : pafSourceBuf = (float *) CPLMalloc(sizeof(float)*nXSize);
1322 3 : GByte* pabyDestBuf1 = (GByte*) CPLMalloc( 4 * nXSize );
1323 3 : GByte* pabyDestBuf2 = pabyDestBuf1 + nXSize;
1324 3 : GByte* pabyDestBuf3 = pabyDestBuf2 + nXSize;
1325 3 : GByte* pabyDestBuf4 = pabyDestBuf3 + nXSize;
1326 : int i, j;
1327 :
1328 3 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
1329 : {
1330 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1331 0 : eErr = CE_Failure;
1332 0 : goto end;
1333 : }
1334 :
1335 366 : for ( i = 0; i < nYSize; i++)
1336 : {
1337 : /* Read source buffer */
1338 : eErr = GDALRasterIO( hSrcBand,
1339 : GF_Read,
1340 : 0, i,
1341 : nXSize, 1,
1342 : (panSourceBuf) ? (void*) panSourceBuf : (void* )pafSourceBuf,
1343 : nXSize, 1,
1344 : (panSourceBuf) ? GDT_Int32 : GDT_Float32,
1345 363 : 0, 0);
1346 363 : if (eErr != CE_None)
1347 0 : goto end;
1348 :
1349 363 : if (panSourceBuf)
1350 : {
1351 0 : for ( j = 0; j < nXSize; j++)
1352 : {
1353 0 : int nIndex = panSourceBuf[j] + nIndexOffset;
1354 0 : pabyDestBuf1[j] = pabyPrecomputed[4 * nIndex];
1355 0 : pabyDestBuf2[j] = pabyPrecomputed[4 * nIndex + 1];
1356 0 : pabyDestBuf3[j] = pabyPrecomputed[4 * nIndex + 2];
1357 0 : pabyDestBuf4[j] = pabyPrecomputed[4 * nIndex + 3];
1358 : }
1359 : }
1360 : else
1361 : {
1362 44286 : for ( j = 0; j < nXSize; j++)
1363 : {
1364 : GDALColorReliefGetRGBA (pasColorAssociation,
1365 : nColorAssociation,
1366 : pafSourceBuf[j],
1367 : eColorSelectionMode,
1368 : &nR,
1369 : &nG,
1370 : &nB,
1371 43923 : &nA);
1372 43923 : pabyDestBuf1[j] = (GByte) nR;
1373 43923 : pabyDestBuf2[j] = (GByte) nG;
1374 43923 : pabyDestBuf3[j] = (GByte) nB;
1375 43923 : pabyDestBuf4[j] = (GByte) nA;
1376 : }
1377 : }
1378 :
1379 : /* -----------------------------------------
1380 : * Write Line to Raster
1381 : */
1382 : eErr = GDALRasterIO(hDstBand1,
1383 : GF_Write,
1384 : 0, i, nXSize,
1385 363 : 1, pabyDestBuf1, nXSize, 1, GDT_Byte, 0, 0);
1386 363 : if (eErr != CE_None)
1387 0 : goto end;
1388 :
1389 : eErr = GDALRasterIO(hDstBand2,
1390 : GF_Write,
1391 : 0, i, nXSize,
1392 363 : 1, pabyDestBuf2, nXSize, 1, GDT_Byte, 0, 0);
1393 363 : if (eErr != CE_None)
1394 0 : goto end;
1395 :
1396 : eErr = GDALRasterIO(hDstBand3,
1397 : GF_Write,
1398 : 0, i, nXSize,
1399 363 : 1, pabyDestBuf3, nXSize, 1, GDT_Byte, 0, 0);
1400 363 : if (eErr != CE_None)
1401 0 : goto end;
1402 :
1403 363 : if (hDstBand4)
1404 : {
1405 : eErr = GDALRasterIO(hDstBand4,
1406 : GF_Write,
1407 : 0, i, nXSize,
1408 0 : 1, pabyDestBuf4, nXSize, 1, GDT_Byte, 0, 0);
1409 0 : if (eErr != CE_None)
1410 0 : goto end;
1411 : }
1412 :
1413 363 : if( !pfnProgress( 1.0 * (i+1) / nYSize, NULL, pProgressData ) )
1414 : {
1415 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1416 0 : eErr = CE_Failure;
1417 0 : goto end;
1418 : }
1419 : }
1420 3 : pfnProgress( 1.0, NULL, pProgressData );
1421 3 : eErr = CE_None;
1422 :
1423 3 : end:
1424 3 : VSIFree(pabyPrecomputed);
1425 3 : CPLFree(pafSourceBuf);
1426 3 : CPLFree(panSourceBuf);
1427 3 : CPLFree(pabyDestBuf1);
1428 3 : CPLFree(pasColorAssociation);
1429 :
1430 3 : return eErr;
1431 : }
1432 :
1433 : /************************************************************************/
1434 : /* GDALGenerateVRTColorRelief() */
1435 : /************************************************************************/
1436 :
1437 : CPLErr GDALGenerateVRTColorRelief(const char* pszDstFilename,
1438 : GDALDatasetH hSrcDataset,
1439 : GDALRasterBandH hSrcBand,
1440 : const char* pszColorFilename,
1441 : ColorSelectionMode eColorSelectionMode,
1442 2 : int bAddAlpha)
1443 : {
1444 :
1445 2 : int nColorAssociation = 0;
1446 : ColorAssociation* pasColorAssociation =
1447 : GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
1448 2 : &nColorAssociation);
1449 2 : if (pasColorAssociation == NULL)
1450 0 : return CE_Failure;
1451 :
1452 2 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
1453 2 : int nYSize = GDALGetRasterBandYSize(hSrcBand);
1454 :
1455 2 : FILE* fp = VSIFOpenL(pszDstFilename, "wt");
1456 2 : if (fp == NULL)
1457 : {
1458 0 : CPLFree(pasColorAssociation);
1459 0 : return CE_Failure;
1460 : }
1461 :
1462 2 : VSIFPrintfL(fp, "<VRTDataset rasterXSize=\"%d\" rasterYSize=\"%d\">\n", nXSize, nYSize);
1463 2 : const char* pszProjectionRef = GDALGetProjectionRef(hSrcDataset);
1464 2 : if (pszProjectionRef && pszProjectionRef[0] != '\0')
1465 : {
1466 2 : char* pszEscapedString = CPLEscapeString(pszProjectionRef, -1, CPLES_XML);
1467 2 : VSIFPrintfL(fp, " <SRS>%s</SRS>\n", pszEscapedString);
1468 2 : VSIFree(pszEscapedString);
1469 : }
1470 : double adfGT[6];
1471 2 : if (GDALGetGeoTransform(hSrcDataset, adfGT) == CE_None)
1472 : {
1473 : VSIFPrintfL(fp, " <GeoTransform> %.16g, %.16g, %.16g, "
1474 : "%.16g, %.16g, %.16g</GeoTransform>\n",
1475 2 : adfGT[0], adfGT[1], adfGT[2], adfGT[3], adfGT[4], adfGT[5]);
1476 : }
1477 2 : int nBands = 3 + (bAddAlpha ? 1 : 0);
1478 : int iBand;
1479 :
1480 : int nBlockXSize, nBlockYSize;
1481 2 : GDALGetBlockSize(hSrcBand, &nBlockXSize, &nBlockYSize);
1482 :
1483 : int bRelativeToVRT;
1484 2 : CPLString osPath = CPLGetPath(pszDstFilename);
1485 : char* pszSourceFilename = CPLStrdup(
1486 : CPLExtractRelativePath( osPath.c_str(), GDALGetDescription(hSrcDataset),
1487 2 : &bRelativeToVRT ));
1488 :
1489 8 : for(iBand = 0; iBand < nBands; iBand++)
1490 : {
1491 6 : VSIFPrintfL(fp, " <VRTRasterBand dataType=\"Byte\" band=\"%d\">\n", iBand + 1);
1492 : VSIFPrintfL(fp, " <ColorInterp>%s</ColorInterp>\n",
1493 6 : GDALGetColorInterpretationName((GDALColorInterp)(GCI_RedBand + iBand)));
1494 6 : VSIFPrintfL(fp, " <ComplexSource>\n");
1495 : VSIFPrintfL(fp, " <SourceFilename relativeToVRT=\"%d\">%s</SourceFilename>\n",
1496 6 : bRelativeToVRT, pszSourceFilename);
1497 6 : VSIFPrintfL(fp, " <SourceBand>%d</SourceBand>\n", GDALGetBandNumber(hSrcBand));
1498 : VSIFPrintfL(fp, " <SourceProperties RasterXSize=\"%d\" "
1499 : "RasterYSize=\"%d\" DataType=\"%s\" "
1500 : "BlockXSize=\"%d\" BlockYSize=\"%d\"/>\n",
1501 : nXSize, nYSize,
1502 : GDALGetDataTypeName(GDALGetRasterDataType(hSrcBand)),
1503 6 : nBlockXSize, nBlockYSize);
1504 : VSIFPrintfL(fp, " <SrcRect xOff=\"0\" yOff=\"0\" xSize=\"%d\" ySize=\"%d\"/>\n",
1505 6 : nXSize, nYSize);
1506 : VSIFPrintfL(fp, " <DstRect xOff=\"0\" yOff=\"0\" xSize=\"%d\" ySize=\"%d\"/>\n",
1507 6 : nXSize, nYSize);
1508 :
1509 6 : VSIFPrintfL(fp, " <LUT>");
1510 : int iColor;
1511 : #define EPSILON 1e-8
1512 48 : for(iColor=0;iColor<nColorAssociation;iColor++)
1513 : {
1514 42 : if (eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY)
1515 : {
1516 21 : if (iColor > 1)
1517 15 : VSIFPrintfL(fp, ",");
1518 : }
1519 21 : else if (iColor > 0)
1520 18 : VSIFPrintfL(fp, ",");
1521 :
1522 42 : double dfVal = pasColorAssociation[iColor].dfVal;
1523 :
1524 42 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
1525 : {
1526 0 : VSIFPrintfL(fp, "%.12g:0,", dfVal - EPSILON);
1527 : }
1528 42 : else if (iColor > 0 &&
1529 : eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY)
1530 : {
1531 18 : double dfMidVal = (dfVal + pasColorAssociation[iColor-1].dfVal) / 2;
1532 : VSIFPrintfL(fp, "%.12g:%d", dfMidVal - EPSILON,
1533 : (iBand == 0) ? pasColorAssociation[iColor-1].nR :
1534 : (iBand == 1) ? pasColorAssociation[iColor-1].nG :
1535 : (iBand == 2) ? pasColorAssociation[iColor-1].nB :
1536 18 : pasColorAssociation[iColor-1].nA);
1537 : VSIFPrintfL(fp, ",%.12g:%d", dfMidVal ,
1538 : (iBand == 0) ? pasColorAssociation[iColor].nR :
1539 : (iBand == 1) ? pasColorAssociation[iColor].nG :
1540 : (iBand == 2) ? pasColorAssociation[iColor].nB :
1541 18 : pasColorAssociation[iColor].nA);
1542 :
1543 : }
1544 :
1545 42 : if (eColorSelectionMode != COLOR_SELECTION_NEAREST_ENTRY)
1546 : {
1547 21 : if (dfVal != (double)(int)dfVal)
1548 0 : VSIFPrintfL(fp, "%.12g", dfVal);
1549 : else
1550 21 : VSIFPrintfL(fp, "%d", (int)dfVal);
1551 : VSIFPrintfL(fp, ":%d",
1552 : (iBand == 0) ? pasColorAssociation[iColor].nR :
1553 : (iBand == 1) ? pasColorAssociation[iColor].nG :
1554 : (iBand == 2) ? pasColorAssociation[iColor].nB :
1555 21 : pasColorAssociation[iColor].nA);
1556 : }
1557 :
1558 42 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
1559 : {
1560 0 : VSIFPrintfL(fp, ",%.12g:0", dfVal + EPSILON);
1561 : }
1562 :
1563 : }
1564 6 : VSIFPrintfL(fp, "</LUT>\n");
1565 :
1566 6 : VSIFPrintfL(fp, " </ComplexSource>\n");
1567 6 : VSIFPrintfL(fp, " </VRTRasterBand>\n");
1568 : }
1569 :
1570 2 : CPLFree(pszSourceFilename);
1571 :
1572 2 : VSIFPrintfL(fp, "</VRTDataset>\n");
1573 :
1574 2 : VSIFCloseL(fp);
1575 :
1576 2 : CPLFree(pasColorAssociation);
1577 :
1578 2 : return CE_None;
1579 : }
1580 :
1581 :
1582 : /************************************************************************/
1583 : /* GDALTRIAlg() */
1584 : /************************************************************************/
1585 :
1586 0 : float GDALTRIAlg (float* afWin, float fDstNoDataValue, void* pData)
1587 : {
1588 : // Terrain Ruggedness is average difference in height
1589 : return (fabs(afWin[0]-afWin[4]) +
1590 : fabs(afWin[1]-afWin[4]) +
1591 : fabs(afWin[2]-afWin[4]) +
1592 : fabs(afWin[3]-afWin[4]) +
1593 : fabs(afWin[5]-afWin[4]) +
1594 : fabs(afWin[6]-afWin[4]) +
1595 : fabs(afWin[7]-afWin[4]) +
1596 0 : fabs(afWin[8]-afWin[4]))/8;
1597 : }
1598 :
1599 :
1600 : /************************************************************************/
1601 : /* GDALTPIAlg() */
1602 : /************************************************************************/
1603 :
1604 0 : float GDALTPIAlg (float* afWin, float fDstNoDataValue, void* pData)
1605 : {
1606 : // Terrain Position is the difference between
1607 : // The central cell and the mean of the surrounding cells
1608 : return afWin[4] -
1609 : ((afWin[0]+
1610 : afWin[1]+
1611 : afWin[2]+
1612 : afWin[3]+
1613 : afWin[5]+
1614 : afWin[6]+
1615 : afWin[7]+
1616 0 : afWin[8])/8);
1617 : }
1618 :
1619 : /************************************************************************/
1620 : /* GDALRoughnessAlg() */
1621 : /************************************************************************/
1622 :
1623 0 : float GDALRoughnessAlg (float* afWin, float fDstNoDataValue, void* pData)
1624 : {
1625 : // Roughness is the largest difference
1626 : // between any two cells
1627 :
1628 0 : float pafRoughnessMin = afWin[0];
1629 0 : float pafRoughnessMax = afWin[0];
1630 :
1631 0 : for ( int k = 1; k < 9; k++)
1632 : {
1633 0 : if (afWin[k] > pafRoughnessMax)
1634 : {
1635 0 : pafRoughnessMax=afWin[k];
1636 : }
1637 0 : if (afWin[k] < pafRoughnessMin)
1638 : {
1639 0 : pafRoughnessMin=afWin[k];
1640 : }
1641 : }
1642 0 : return pafRoughnessMax - pafRoughnessMin;
1643 : }
1644 :
1645 : /************************************************************************/
1646 : /* ==================================================================== */
1647 : /* GDALGeneric3x3Dataset */
1648 : /* ==================================================================== */
1649 : /************************************************************************/
1650 :
1651 : class GDALGeneric3x3RasterBand;
1652 :
1653 : class GDALGeneric3x3Dataset : public GDALDataset
1654 : {
1655 : friend class GDALGeneric3x3RasterBand;
1656 :
1657 : GDALGeneric3x3ProcessingAlg pfnAlg;
1658 : void* pAlgData;
1659 : GDALDatasetH hSrcDS;
1660 : GDALRasterBandH hSrcBand;
1661 : float* apafSourceBuf[3];
1662 : int bDstHasNoData;
1663 : double dfDstNoDataValue;
1664 : int nCurLine;
1665 : int bComputeAtEdges;
1666 :
1667 : public:
1668 : GDALGeneric3x3Dataset(GDALDatasetH hSrcDS,
1669 : GDALRasterBandH hSrcBand,
1670 : GDALDataType eDstDataType,
1671 : int bDstHasNoData,
1672 : double dfDstNoDataValue,
1673 : GDALGeneric3x3ProcessingAlg pfnAlg,
1674 : void* pAlgData,
1675 : int bComputeAtEdges);
1676 : ~GDALGeneric3x3Dataset();
1677 :
1678 : CPLErr GetGeoTransform( double * padfGeoTransform );
1679 : const char *GetProjectionRef();
1680 : };
1681 :
1682 : /************************************************************************/
1683 : /* ==================================================================== */
1684 : /* GDALGeneric3x3RasterBand */
1685 : /* ==================================================================== */
1686 : /************************************************************************/
1687 :
1688 : class GDALGeneric3x3RasterBand : public GDALRasterBand
1689 2 : {
1690 : friend class GDALGeneric3x3Dataset;
1691 : int bSrcHasNoData;
1692 : float fSrcNoDataValue;
1693 :
1694 : void InitWidthNoData(void* pImage);
1695 :
1696 : public:
1697 : GDALGeneric3x3RasterBand( GDALGeneric3x3Dataset *poDS,
1698 : GDALDataType eDstDataType );
1699 :
1700 : virtual CPLErr IReadBlock( int, int, void * );
1701 : virtual double GetNoDataValue( int* pbHasNoData );
1702 : };
1703 :
1704 : GDALGeneric3x3Dataset::GDALGeneric3x3Dataset(
1705 : GDALDatasetH hSrcDS,
1706 : GDALRasterBandH hSrcBand,
1707 : GDALDataType eDstDataType,
1708 : int bDstHasNoData,
1709 : double dfDstNoDataValue,
1710 : GDALGeneric3x3ProcessingAlg pfnAlg,
1711 : void* pAlgData,
1712 2 : int bComputeAtEdges)
1713 : {
1714 2 : this->hSrcDS = hSrcDS;
1715 2 : this->hSrcBand = hSrcBand;
1716 2 : this->pfnAlg = pfnAlg;
1717 2 : this->pAlgData = pAlgData;
1718 2 : this->bDstHasNoData = bDstHasNoData;
1719 2 : this->dfDstNoDataValue = dfDstNoDataValue;
1720 2 : this->bComputeAtEdges = bComputeAtEdges;
1721 :
1722 2 : CPLAssert(eDstDataType == GDT_Byte || eDstDataType == GDT_Float32);
1723 :
1724 2 : nRasterXSize = GDALGetRasterXSize(hSrcDS);
1725 2 : nRasterYSize = GDALGetRasterYSize(hSrcDS);
1726 :
1727 2 : SetBand(1, new GDALGeneric3x3RasterBand(this, eDstDataType));
1728 :
1729 2 : apafSourceBuf[0] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1730 2 : apafSourceBuf[1] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1731 2 : apafSourceBuf[2] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1732 :
1733 2 : nCurLine = -1;
1734 2 : }
1735 :
1736 2 : GDALGeneric3x3Dataset::~GDALGeneric3x3Dataset()
1737 : {
1738 2 : CPLFree(apafSourceBuf[0]);
1739 2 : CPLFree(apafSourceBuf[1]);
1740 2 : CPLFree(apafSourceBuf[2]);
1741 2 : }
1742 :
1743 2 : CPLErr GDALGeneric3x3Dataset::GetGeoTransform( double * padfGeoTransform )
1744 : {
1745 2 : return GDALGetGeoTransform(hSrcDS, padfGeoTransform);
1746 : }
1747 :
1748 2 : const char *GDALGeneric3x3Dataset::GetProjectionRef()
1749 : {
1750 2 : return GDALGetProjectionRef(hSrcDS);
1751 : }
1752 :
1753 : GDALGeneric3x3RasterBand::GDALGeneric3x3RasterBand(GDALGeneric3x3Dataset *poDS,
1754 2 : GDALDataType eDstDataType)
1755 : {
1756 2 : this->poDS = poDS;
1757 2 : this->nBand = 1;
1758 2 : eDataType = eDstDataType;
1759 2 : nBlockXSize = poDS->GetRasterXSize();
1760 2 : nBlockYSize = 1;
1761 :
1762 2 : bSrcHasNoData = FALSE;
1763 : fSrcNoDataValue = (float)GDALGetRasterNoDataValue(poDS->hSrcBand,
1764 2 : &bSrcHasNoData);
1765 2 : }
1766 :
1767 2 : void GDALGeneric3x3RasterBand::InitWidthNoData(void* pImage)
1768 : {
1769 : int j;
1770 2 : GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
1771 2 : if (eDataType == GDT_Byte)
1772 : {
1773 244 : for(j=0;j<nBlockXSize;j++)
1774 242 : ((GByte*)pImage)[j] = (GByte) poGDS->dfDstNoDataValue;
1775 : }
1776 : else
1777 : {
1778 0 : for(j=0;j<nBlockXSize;j++)
1779 0 : ((float*)pImage)[j] = (float) poGDS->dfDstNoDataValue;
1780 : }
1781 2 : }
1782 :
1783 : CPLErr GDALGeneric3x3RasterBand::IReadBlock( int nBlockXOff,
1784 : int nBlockYOff,
1785 242 : void *pImage )
1786 : {
1787 : int i, j;
1788 : float fVal;
1789 242 : GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
1790 :
1791 361 : if (poGDS->bComputeAtEdges && nRasterXSize >= 2 && nRasterYSize >= 2)
1792 : {
1793 121 : if (nBlockYOff == 0)
1794 : {
1795 3 : for(i=0;i<2;i++)
1796 : {
1797 : CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1798 : GF_Read,
1799 : 0, i, nBlockXSize, 1,
1800 : poGDS->apafSourceBuf[i+1],
1801 : nBlockXSize, 1,
1802 : GDT_Float32,
1803 2 : 0, 0);
1804 2 : if (eErr != CE_None)
1805 : {
1806 0 : InitWidthNoData(pImage);
1807 0 : return eErr;
1808 : }
1809 : }
1810 1 : poGDS->nCurLine = 0;
1811 :
1812 122 : for (j = 0; j < nRasterXSize; j++)
1813 : {
1814 : float afWin[9];
1815 121 : int jmin = (j == 0) ? j : j - 1;
1816 121 : int jmax = (j == nRasterXSize - 1) ? j : j + 1;
1817 :
1818 121 : afWin[0] = INTERPOL(poGDS->apafSourceBuf[1][jmin], poGDS->apafSourceBuf[2][jmin]);
1819 121 : afWin[1] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[2][j]);
1820 121 : afWin[2] = INTERPOL(poGDS->apafSourceBuf[1][jmax], poGDS->apafSourceBuf[2][jmax]);
1821 121 : afWin[3] = poGDS->apafSourceBuf[1][jmin];
1822 121 : afWin[4] = poGDS->apafSourceBuf[1][j];
1823 121 : afWin[5] = poGDS->apafSourceBuf[1][jmax];
1824 121 : afWin[6] = poGDS->apafSourceBuf[2][jmin];
1825 121 : afWin[7] = poGDS->apafSourceBuf[2][j];
1826 121 : afWin[8] = poGDS->apafSourceBuf[2][jmax];
1827 :
1828 : fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
1829 : afWin, poGDS->dfDstNoDataValue,
1830 : poGDS->pfnAlg,
1831 : poGDS->pAlgData,
1832 121 : poGDS->bComputeAtEdges);
1833 :
1834 121 : if (eDataType == GDT_Byte)
1835 121 : ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
1836 : else
1837 0 : ((float*)pImage)[j] = fVal;
1838 : }
1839 :
1840 1 : return CE_None;
1841 : }
1842 120 : else if (nBlockYOff == nRasterYSize - 1)
1843 : {
1844 1 : if (poGDS->nCurLine != nRasterYSize - 2)
1845 : {
1846 0 : for(i=0;i<2;i++)
1847 : {
1848 : CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1849 : GF_Read,
1850 : 0, nRasterYSize - 2 + i, nBlockXSize, 1,
1851 : poGDS->apafSourceBuf[i+1],
1852 : nBlockXSize, 1,
1853 : GDT_Float32,
1854 0 : 0, 0);
1855 0 : if (eErr != CE_None)
1856 : {
1857 0 : InitWidthNoData(pImage);
1858 0 : return eErr;
1859 : }
1860 : }
1861 : }
1862 :
1863 122 : for (j = 0; j < nRasterXSize; j++)
1864 : {
1865 : float afWin[9];
1866 121 : int jmin = (j == 0) ? j : j - 1;
1867 121 : int jmax = (j == nRasterXSize - 1) ? j : j + 1;
1868 :
1869 121 : afWin[0] = poGDS->apafSourceBuf[1][jmin];
1870 121 : afWin[1] = poGDS->apafSourceBuf[1][j];
1871 121 : afWin[2] = poGDS->apafSourceBuf[1][jmax];
1872 121 : afWin[3] = poGDS->apafSourceBuf[2][jmin];
1873 121 : afWin[4] = poGDS->apafSourceBuf[2][j];
1874 121 : afWin[5] = poGDS->apafSourceBuf[2][jmax];
1875 121 : afWin[6] = INTERPOL(poGDS->apafSourceBuf[2][jmin], poGDS->apafSourceBuf[1][jmin]);
1876 121 : afWin[7] = INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[1][j]);
1877 121 : afWin[8] = INTERPOL(poGDS->apafSourceBuf[2][jmax], poGDS->apafSourceBuf[1][jmax]);
1878 :
1879 : fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
1880 : afWin, poGDS->dfDstNoDataValue,
1881 : poGDS->pfnAlg,
1882 : poGDS->pAlgData,
1883 121 : poGDS->bComputeAtEdges);
1884 :
1885 121 : if (eDataType == GDT_Byte)
1886 121 : ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
1887 : else
1888 0 : ((float*)pImage)[j] = fVal;
1889 : }
1890 :
1891 1 : return CE_None;
1892 : }
1893 : }
1894 121 : else if ( nBlockYOff == 0 || nBlockYOff == nRasterYSize - 1)
1895 : {
1896 2 : InitWidthNoData(pImage);
1897 2 : return CE_None;
1898 : }
1899 :
1900 238 : if ( poGDS->nCurLine != nBlockYOff )
1901 : {
1902 238 : if (poGDS->nCurLine + 1 == nBlockYOff)
1903 : {
1904 237 : float* pafTmp = poGDS->apafSourceBuf[0];
1905 237 : poGDS->apafSourceBuf[0] = poGDS->apafSourceBuf[1];
1906 237 : poGDS->apafSourceBuf[1] = poGDS->apafSourceBuf[2];
1907 237 : poGDS->apafSourceBuf[2] = pafTmp;
1908 :
1909 : CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1910 : GF_Read,
1911 : 0, nBlockYOff + 1, nBlockXSize, 1,
1912 : poGDS->apafSourceBuf[2],
1913 : nBlockXSize, 1,
1914 : GDT_Float32,
1915 237 : 0, 0);
1916 :
1917 237 : if (eErr != CE_None)
1918 : {
1919 0 : InitWidthNoData(pImage);
1920 0 : return eErr;
1921 : }
1922 : }
1923 : else
1924 : {
1925 4 : for(i=0;i<3;i++)
1926 : {
1927 : CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1928 : GF_Read,
1929 : 0, nBlockYOff + i - 1, nBlockXSize, 1,
1930 : poGDS->apafSourceBuf[i],
1931 : nBlockXSize, 1,
1932 : GDT_Float32,
1933 3 : 0, 0);
1934 3 : if (eErr != CE_None)
1935 : {
1936 0 : InitWidthNoData(pImage);
1937 0 : return eErr;
1938 : }
1939 : }
1940 : }
1941 :
1942 238 : poGDS->nCurLine = nBlockYOff;
1943 : }
1944 :
1945 357 : if (poGDS->bComputeAtEdges && nRasterXSize >= 2)
1946 : {
1947 : float afWin[9];
1948 :
1949 119 : j = 0;
1950 119 : afWin[0] = INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j+1]);
1951 119 : afWin[1] = poGDS->apafSourceBuf[0][j];
1952 119 : afWin[2] = poGDS->apafSourceBuf[0][j+1];
1953 119 : afWin[3] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j+1]);
1954 119 : afWin[4] = poGDS->apafSourceBuf[1][j];
1955 119 : afWin[5] = poGDS->apafSourceBuf[1][j+1];
1956 119 : afWin[6] = INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j+1]);
1957 119 : afWin[7] = poGDS->apafSourceBuf[2][j];
1958 119 : afWin[8] = poGDS->apafSourceBuf[2][j+1];
1959 :
1960 : fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
1961 : afWin, poGDS->dfDstNoDataValue,
1962 : poGDS->pfnAlg,
1963 : poGDS->pAlgData,
1964 119 : poGDS->bComputeAtEdges);
1965 :
1966 119 : if (eDataType == GDT_Byte)
1967 119 : ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
1968 : else
1969 0 : ((float*)pImage)[j] = fVal;
1970 :
1971 119 : j = nRasterXSize - 1;
1972 :
1973 119 : afWin[0] = poGDS->apafSourceBuf[0][j-1];
1974 119 : afWin[1] = poGDS->apafSourceBuf[0][j];
1975 119 : afWin[2] = INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j-1]);
1976 119 : afWin[3] = poGDS->apafSourceBuf[1][j-1];
1977 119 : afWin[4] = poGDS->apafSourceBuf[1][j];
1978 119 : afWin[5] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j-1]);
1979 119 : afWin[6] = poGDS->apafSourceBuf[2][j-1];
1980 119 : afWin[7] = poGDS->apafSourceBuf[2][j];
1981 119 : afWin[8] = INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j-1]);
1982 :
1983 : fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
1984 : afWin, poGDS->dfDstNoDataValue,
1985 : poGDS->pfnAlg,
1986 : poGDS->pAlgData,
1987 119 : poGDS->bComputeAtEdges);
1988 :
1989 119 : if (eDataType == GDT_Byte)
1990 119 : ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
1991 : else
1992 0 : ((float*)pImage)[j] = fVal;
1993 : }
1994 : else
1995 : {
1996 119 : if (eDataType == GDT_Byte)
1997 : {
1998 119 : ((GByte*)pImage)[0] = (GByte) poGDS->dfDstNoDataValue;
1999 119 : if (nBlockXSize > 1)
2000 119 : ((GByte*)pImage)[nBlockXSize - 1] = (GByte) poGDS->dfDstNoDataValue;
2001 : }
2002 : else
2003 : {
2004 0 : ((float*)pImage)[0] = (float) poGDS->dfDstNoDataValue;
2005 0 : if (nBlockXSize > 1)
2006 0 : ((float*)pImage)[nBlockXSize - 1] = (float) poGDS->dfDstNoDataValue;
2007 : }
2008 : }
2009 :
2010 :
2011 28560 : for(j=1;j<nBlockXSize - 1;j++)
2012 : {
2013 : float afWin[9];
2014 28322 : afWin[0] = poGDS->apafSourceBuf[0][j-1];
2015 28322 : afWin[1] = poGDS->apafSourceBuf[0][j];
2016 28322 : afWin[2] = poGDS->apafSourceBuf[0][j+1];
2017 28322 : afWin[3] = poGDS->apafSourceBuf[1][j-1];
2018 28322 : afWin[4] = poGDS->apafSourceBuf[1][j];
2019 28322 : afWin[5] = poGDS->apafSourceBuf[1][j+1];
2020 28322 : afWin[6] = poGDS->apafSourceBuf[2][j-1];
2021 28322 : afWin[7] = poGDS->apafSourceBuf[2][j];
2022 28322 : afWin[8] = poGDS->apafSourceBuf[2][j+1];
2023 :
2024 : fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
2025 : afWin, poGDS->dfDstNoDataValue,
2026 : poGDS->pfnAlg,
2027 : poGDS->pAlgData,
2028 28322 : poGDS->bComputeAtEdges);
2029 :
2030 28322 : if (eDataType == GDT_Byte)
2031 28322 : ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
2032 : else
2033 0 : ((float*)pImage)[j] = fVal;
2034 :
2035 : }
2036 :
2037 238 : return CE_None;
2038 : }
2039 :
2040 8 : double GDALGeneric3x3RasterBand::GetNoDataValue( int* pbHasNoData )
2041 : {
2042 8 : GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
2043 8 : if (pbHasNoData)
2044 6 : *pbHasNoData = poGDS->bDstHasNoData;
2045 8 : return poGDS->dfDstNoDataValue;
2046 : }
2047 :
2048 : /************************************************************************/
2049 : /* main() */
2050 : /************************************************************************/
2051 :
2052 : typedef enum
2053 : {
2054 : HILL_SHADE,
2055 : SLOPE,
2056 : ASPECT,
2057 : COLOR_RELIEF,
2058 : TRI,
2059 : TPI,
2060 : ROUGHNESS
2061 : } Algorithm;
2062 :
2063 15 : int main( int argc, char ** argv )
2064 :
2065 : {
2066 : Algorithm eUtilityMode;
2067 15 : double z = 1.0;
2068 15 : double scale = 1.0;
2069 15 : double az = 315.0;
2070 15 : double alt = 45.0;
2071 : // 0 = 'percent' or 1 = 'degrees'
2072 15 : int slopeFormat = 1;
2073 15 : int bAddAlpha = FALSE;
2074 15 : int bZeroForFlat = FALSE;
2075 15 : int bAngleAsAzimuth = TRUE;
2076 15 : ColorSelectionMode eColorSelectionMode = COLOR_SELECTION_INTERPOLATE;
2077 :
2078 15 : int nBand = 1;
2079 : double adfGeoTransform[6];
2080 :
2081 15 : const char *pszSrcFilename = NULL;
2082 15 : const char *pszDstFilename = NULL;
2083 15 : const char *pszColorFilename = NULL;
2084 15 : const char *pszFormat = "GTiff";
2085 15 : char **papszCreateOptions = NULL;
2086 :
2087 15 : GDALDatasetH hSrcDataset = NULL;
2088 15 : GDALDatasetH hDstDataset = NULL;
2089 15 : GDALRasterBandH hSrcBand = NULL;
2090 15 : GDALRasterBandH hDstBand = NULL;
2091 15 : GDALDriverH hDriver = NULL;
2092 :
2093 15 : GDALProgressFunc pfnProgress = GDALTermProgress;
2094 :
2095 15 : int nXSize = 0;
2096 15 : int nYSize = 0;
2097 :
2098 15 : int bComputeAtEdges = FALSE;
2099 15 : int bZevenbergenThorne = FALSE;
2100 :
2101 : /* Check strict compilation and runtime library version as we use C++ API */
2102 15 : if (! GDAL_CHECK_VERSION(argv[0]))
2103 0 : exit(1);
2104 :
2105 15 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
2106 15 : if( argc < 2 )
2107 : {
2108 0 : fprintf(stderr, "Not enough arguments\n");
2109 0 : Usage();
2110 0 : exit( 1 );
2111 : }
2112 :
2113 15 : if( EQUAL(argv[1], "--utility_version") || EQUAL(argv[1], "--utility-version") )
2114 : {
2115 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
2116 1 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
2117 1 : return 0;
2118 : }
2119 19 : else if ( EQUAL(argv[1], "shade") || EQUAL(argv[1], "hillshade") )
2120 : {
2121 5 : eUtilityMode = HILL_SHADE;
2122 : }
2123 9 : else if ( EQUAL(argv[1], "slope") )
2124 : {
2125 1 : eUtilityMode = SLOPE;
2126 : }
2127 8 : else if ( EQUAL(argv[1], "aspect") )
2128 : {
2129 1 : eUtilityMode = ASPECT;
2130 : }
2131 7 : else if ( EQUAL(argv[1], "color-relief") )
2132 : {
2133 7 : eUtilityMode = COLOR_RELIEF;
2134 : }
2135 0 : else if ( EQUAL(argv[1], "TRI") )
2136 : {
2137 0 : eUtilityMode = TRI;
2138 : }
2139 0 : else if ( EQUAL(argv[1], "TPI") )
2140 : {
2141 0 : eUtilityMode = TPI;
2142 : }
2143 0 : else if ( EQUAL(argv[1], "roughness") )
2144 : {
2145 0 : eUtilityMode = ROUGHNESS;
2146 : }
2147 : else
2148 : {
2149 0 : fprintf(stderr, "Missing valid sub-utility mention\n");
2150 0 : Usage();
2151 0 : exit( 1 );
2152 : }
2153 :
2154 : /* -------------------------------------------------------------------- */
2155 : /* Parse arguments. */
2156 : /* -------------------------------------------------------------------- */
2157 72 : for(int i = 2; i < argc; i++ )
2158 : {
2159 63 : if( eUtilityMode == HILL_SHADE && i + 1 < argc &&
2160 : (EQUAL(argv[i], "--z") || EQUAL(argv[i], "-z")))
2161 : {
2162 5 : z = atof(argv[++i]);
2163 : }
2164 53 : else if ( eUtilityMode == SLOPE && EQUAL(argv[i], "-p"))
2165 : {
2166 0 : slopeFormat = 0;
2167 : }
2168 53 : else if ( (eUtilityMode == HILL_SHADE || eUtilityMode == SLOPE ||
2169 : eUtilityMode == ASPECT) && EQUAL(argv[i], "-alg") && i + 1 < argc)
2170 : {
2171 0 : i ++;
2172 0 : if (EQUAL(argv[i], "ZevenbergenThorne"))
2173 0 : bZevenbergenThorne = TRUE;
2174 0 : else if (!EQUAL(argv[i], "Horn"))
2175 : {
2176 0 : fprintf(stderr, "Wrong value for alg : %s\n", argv[i]);
2177 0 : Usage();
2178 : }
2179 : }
2180 53 : else if ( eUtilityMode == ASPECT && EQUAL(argv[i], "-trigonometric"))
2181 : {
2182 0 : bAngleAsAzimuth = FALSE;
2183 : }
2184 53 : else if ( eUtilityMode == ASPECT && EQUAL(argv[i], "-zero_for_flat"))
2185 : {
2186 0 : bZeroForFlat = TRUE;
2187 : }
2188 53 : else if ( eUtilityMode == COLOR_RELIEF && EQUAL(argv[i], "-exact_color_entry"))
2189 : {
2190 0 : eColorSelectionMode = COLOR_SELECTION_EXACT_ENTRY;
2191 : }
2192 55 : else if ( eUtilityMode == COLOR_RELIEF && EQUAL(argv[i], "-nearest_color_entry"))
2193 : {
2194 2 : eColorSelectionMode = COLOR_SELECTION_NEAREST_ENTRY;
2195 : }
2196 57 : else if( i + 1 < argc &&
2197 : (EQUAL(argv[i], "--s") ||
2198 : EQUAL(argv[i], "-s") ||
2199 : EQUAL(argv[i], "--scale") ||
2200 : EQUAL(argv[i], "-scale"))
2201 : )
2202 : {
2203 6 : scale = atof(argv[++i]);
2204 : }
2205 46 : else if( eUtilityMode == HILL_SHADE && i + 1 < argc &&
2206 : (EQUAL(argv[i], "--az") ||
2207 : EQUAL(argv[i], "-az") ||
2208 : EQUAL(argv[i], "--azimuth") ||
2209 : EQUAL(argv[i], "-azimuth"))
2210 : )
2211 : {
2212 1 : az = atof(argv[++i]);
2213 : }
2214 44 : else if( eUtilityMode == HILL_SHADE && i + 1 < argc &&
2215 : (EQUAL(argv[i], "--alt") ||
2216 : EQUAL(argv[i], "-alt") ||
2217 : EQUAL(argv[i], "--alt") ||
2218 : EQUAL(argv[i], "-alt"))
2219 : )
2220 : {
2221 0 : alt = atof(argv[++i]);
2222 : }
2223 44 : else if( eUtilityMode == COLOR_RELIEF &&
2224 : EQUAL(argv[i], "-alpha"))
2225 : {
2226 0 : bAddAlpha = TRUE;
2227 : }
2228 46 : else if( eUtilityMode != COLOR_RELIEF &&
2229 : EQUAL(argv[i], "-compute_edges"))
2230 : {
2231 2 : bComputeAtEdges = TRUE;
2232 : }
2233 42 : else if( i + 1 < argc &&
2234 : (EQUAL(argv[i], "--b") ||
2235 : EQUAL(argv[i], "-b"))
2236 : )
2237 : {
2238 0 : nBand = atoi(argv[++i]);
2239 : }
2240 42 : else if ( EQUAL(argv[i], "-q") || EQUAL(argv[i], "-quiet") )
2241 : {
2242 0 : pfnProgress = GDALDummyProgress;
2243 : }
2244 43 : else if( EQUAL(argv[i],"-co") && i < argc-1 )
2245 : {
2246 1 : papszCreateOptions = CSLAddString( papszCreateOptions, argv[++i] );
2247 : }
2248 47 : else if( EQUAL(argv[i],"-of") && i < argc-1 )
2249 : {
2250 6 : pszFormat = argv[++i];
2251 : }
2252 35 : else if( argv[i][0] == '-' )
2253 : {
2254 : fprintf( stderr, "Option %s incomplete, or not recognised.\n\n",
2255 0 : argv[i] );
2256 0 : Usage();
2257 0 : GDALDestroyDriverManager();
2258 0 : exit( 2 );
2259 : }
2260 35 : else if( pszSrcFilename == NULL )
2261 : {
2262 14 : pszSrcFilename = argv[i];
2263 : }
2264 28 : else if( eUtilityMode == COLOR_RELIEF && pszColorFilename == NULL )
2265 : {
2266 7 : pszColorFilename = argv[i];
2267 : }
2268 14 : else if( pszDstFilename == NULL )
2269 : {
2270 14 : pszDstFilename = argv[i];
2271 : }
2272 : else
2273 0 : Usage();
2274 : }
2275 :
2276 14 : if( pszSrcFilename == NULL )
2277 : {
2278 0 : fprintf( stderr, "Missing source.\n\n" );
2279 0 : Usage();
2280 : }
2281 14 : if ( eUtilityMode == COLOR_RELIEF && pszColorFilename == NULL )
2282 : {
2283 0 : fprintf( stderr, "Missing color file.\n\n" );
2284 0 : Usage();
2285 : }
2286 14 : if( pszDstFilename == NULL )
2287 : {
2288 0 : fprintf( stderr, "Missing destination.\n\n" );
2289 0 : Usage();
2290 : }
2291 :
2292 14 : GDALAllRegister();
2293 :
2294 : // Open Dataset and get raster band
2295 14 : hSrcDataset = GDALOpen( pszSrcFilename, GA_ReadOnly );
2296 :
2297 14 : if( hSrcDataset == NULL )
2298 : {
2299 : fprintf( stderr,
2300 : "GDALOpen failed - %d\n%s\n",
2301 0 : CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
2302 0 : GDALDestroyDriverManager();
2303 0 : exit( 1 );
2304 : }
2305 :
2306 14 : nXSize = GDALGetRasterXSize(hSrcDataset);
2307 14 : nYSize = GDALGetRasterYSize(hSrcDataset);
2308 :
2309 14 : if( nBand <= 0 || nBand > GDALGetRasterCount(hSrcDataset) )
2310 : {
2311 : fprintf( stderr,
2312 0 : "Unable to fetch band #%d\n", nBand );
2313 0 : GDALDestroyDriverManager();
2314 0 : exit( 1 );
2315 : }
2316 14 : hSrcBand = GDALGetRasterBand( hSrcDataset, nBand );
2317 :
2318 14 : GDALGetGeoTransform(hSrcDataset, adfGeoTransform);
2319 :
2320 14 : hDriver = GDALGetDriverByName(pszFormat);
2321 14 : if( hDriver == NULL
2322 : || (GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL &&
2323 : GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, NULL ) == NULL))
2324 : {
2325 : int iDr;
2326 :
2327 : fprintf( stderr, "Output driver `%s' not recognised to have\n",
2328 0 : pszFormat );
2329 : fprintf( stderr, "output support. The following format drivers are configured\n"
2330 0 : "and support output:\n" );
2331 :
2332 0 : for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
2333 : {
2334 0 : GDALDriverH hDriver = GDALGetDriver(iDr);
2335 :
2336 0 : if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) != NULL ||
2337 : GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, NULL ) != NULL )
2338 : {
2339 : printf( " %s: %s\n",
2340 : GDALGetDriverShortName( hDriver ),
2341 0 : GDALGetDriverLongName( hDriver ) );
2342 : }
2343 : }
2344 0 : GDALDestroyDriverManager();
2345 0 : exit( 1 );
2346 : }
2347 :
2348 14 : double dfDstNoDataValue = 0;
2349 14 : int bDstHasNoData = FALSE;
2350 14 : void* pData = NULL;
2351 14 : GDALGeneric3x3ProcessingAlg pfnAlg = NULL;
2352 :
2353 14 : if (eUtilityMode == HILL_SHADE)
2354 : {
2355 5 : dfDstNoDataValue = 0;
2356 5 : bDstHasNoData = TRUE;
2357 : pData = GDALCreateHillshadeData (adfGeoTransform,
2358 : z,
2359 : scale,
2360 : alt,
2361 : az,
2362 5 : bZevenbergenThorne);
2363 5 : if (bZevenbergenThorne)
2364 0 : pfnAlg = GDALHillshadeZevenbergenThorneAlg;
2365 : else
2366 5 : pfnAlg = GDALHillshadeAlg;
2367 : }
2368 9 : else if (eUtilityMode == SLOPE)
2369 : {
2370 1 : dfDstNoDataValue = -9999;
2371 1 : bDstHasNoData = TRUE;
2372 :
2373 1 : pData = GDALCreateSlopeData(adfGeoTransform, scale, slopeFormat);
2374 1 : if (bZevenbergenThorne)
2375 0 : pfnAlg = GDALSlopeZevenbergenThorneAlg;
2376 : else
2377 1 : pfnAlg = GDALSlopeHornAlg;
2378 : }
2379 :
2380 8 : else if (eUtilityMode == ASPECT)
2381 : {
2382 1 : if (!bZeroForFlat)
2383 : {
2384 1 : dfDstNoDataValue = -9999;
2385 1 : bDstHasNoData = TRUE;
2386 : }
2387 :
2388 1 : pData = GDALCreateAspectData(bAngleAsAzimuth);
2389 1 : if (bZevenbergenThorne)
2390 0 : pfnAlg = GDALAspectZevenbergenThorneAlg;
2391 : else
2392 1 : pfnAlg = GDALAspectAlg;
2393 : }
2394 7 : else if (eUtilityMode == TRI)
2395 : {
2396 0 : dfDstNoDataValue = -9999;
2397 0 : bDstHasNoData = TRUE;
2398 0 : pfnAlg = GDALTRIAlg;
2399 : }
2400 7 : else if (eUtilityMode == TPI)
2401 : {
2402 0 : dfDstNoDataValue = -9999;
2403 0 : bDstHasNoData = TRUE;
2404 0 : pfnAlg = GDALTPIAlg;
2405 : }
2406 7 : else if (eUtilityMode == ROUGHNESS)
2407 : {
2408 0 : dfDstNoDataValue = -9999;
2409 0 : bDstHasNoData = TRUE;
2410 0 : pfnAlg = GDALRoughnessAlg;
2411 : }
2412 :
2413 : GDALDataType eDstDataType = (eUtilityMode == HILL_SHADE ||
2414 : eUtilityMode == COLOR_RELIEF) ? GDT_Byte :
2415 14 : GDT_Float32;
2416 :
2417 14 : if( EQUAL(pszFormat, "VRT") )
2418 : {
2419 2 : if (eUtilityMode == COLOR_RELIEF)
2420 : {
2421 : GDALGenerateVRTColorRelief(pszDstFilename,
2422 : hSrcDataset,
2423 : hSrcBand,
2424 : pszColorFilename,
2425 : eColorSelectionMode,
2426 2 : bAddAlpha);
2427 2 : GDALClose(hSrcDataset);
2428 :
2429 2 : CPLFree(pData);
2430 :
2431 2 : GDALDestroyDriverManager();
2432 2 : CSLDestroy( argv );
2433 2 : CSLDestroy( papszCreateOptions );
2434 :
2435 2 : return 0;
2436 : }
2437 : else
2438 : {
2439 0 : fprintf(stderr, "VRT driver can only be used with color-relief utility\n");
2440 0 : GDALDestroyDriverManager();
2441 :
2442 0 : exit(1);
2443 : }
2444 : }
2445 :
2446 12 : if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL &&
2447 : GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, NULL ) != NULL)
2448 : {
2449 : GDALDatasetH hIntermediateDataset;
2450 :
2451 4 : if (eUtilityMode == COLOR_RELIEF)
2452 : hIntermediateDataset = (GDALDatasetH)
2453 : new GDALColorReliefDataset (hSrcDataset,
2454 : hSrcBand,
2455 : pszColorFilename,
2456 : eColorSelectionMode,
2457 2 : bAddAlpha);
2458 : else
2459 : hIntermediateDataset = (GDALDatasetH)
2460 : new GDALGeneric3x3Dataset(hSrcDataset, hSrcBand,
2461 : eDstDataType,
2462 : bDstHasNoData,
2463 : dfDstNoDataValue,
2464 : pfnAlg,
2465 : pData,
2466 2 : bComputeAtEdges);
2467 :
2468 : GDALDatasetH hOutDS = GDALCreateCopy(
2469 : hDriver, pszDstFilename, hIntermediateDataset,
2470 : TRUE, papszCreateOptions,
2471 4 : pfnProgress, NULL );
2472 :
2473 4 : if( hOutDS != NULL )
2474 4 : GDALClose( hOutDS );
2475 4 : GDALClose(hIntermediateDataset);
2476 4 : GDALClose(hSrcDataset);
2477 :
2478 4 : CPLFree(pData);
2479 :
2480 4 : GDALDestroyDriverManager();
2481 4 : CSLDestroy( argv );
2482 4 : CSLDestroy( papszCreateOptions );
2483 :
2484 4 : return 0;
2485 : }
2486 :
2487 : int nDstBands;
2488 8 : if (eUtilityMode == COLOR_RELIEF)
2489 3 : nDstBands = (bAddAlpha) ? 4 : 3;
2490 : else
2491 5 : nDstBands = 1;
2492 :
2493 : hDstDataset = GDALCreate( hDriver,
2494 : pszDstFilename,
2495 : nXSize,
2496 : nYSize,
2497 : nDstBands,
2498 : eDstDataType,
2499 8 : papszCreateOptions);
2500 :
2501 8 : if( hDstDataset == NULL )
2502 : {
2503 : fprintf( stderr,
2504 : "Unable to create dataset %s %d\n%s\n", pszDstFilename,
2505 0 : CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
2506 0 : GDALDestroyDriverManager();
2507 0 : exit( 1 );
2508 : }
2509 :
2510 8 : hDstBand = GDALGetRasterBand( hDstDataset, 1 );
2511 :
2512 8 : GDALSetGeoTransform(hDstDataset, adfGeoTransform);
2513 8 : GDALSetProjection(hDstDataset, GDALGetProjectionRef(hSrcDataset));
2514 :
2515 8 : if (eUtilityMode == COLOR_RELIEF)
2516 : {
2517 : GDALColorRelief (hSrcBand,
2518 : GDALGetRasterBand(hDstDataset, 1),
2519 : GDALGetRasterBand(hDstDataset, 2),
2520 : GDALGetRasterBand(hDstDataset, 3),
2521 : (bAddAlpha) ? GDALGetRasterBand(hDstDataset, 4) : NULL,
2522 : pszColorFilename,
2523 : eColorSelectionMode,
2524 3 : pfnProgress, NULL);
2525 : }
2526 : else
2527 : {
2528 5 : if (bDstHasNoData)
2529 5 : GDALSetRasterNoDataValue(hDstBand, dfDstNoDataValue);
2530 :
2531 : GDALGeneric3x3Processing(hSrcBand, hDstBand,
2532 : pfnAlg, pData,
2533 : bComputeAtEdges,
2534 5 : pfnProgress, NULL);
2535 :
2536 : }
2537 :
2538 8 : GDALClose(hSrcDataset);
2539 8 : GDALClose(hDstDataset);
2540 8 : CPLFree(pData);
2541 :
2542 8 : GDALDestroyDriverManager();
2543 8 : CSLDestroy( argv );
2544 8 : CSLDestroy( papszCreateOptions );
2545 :
2546 8 : return 0;
2547 : }
2548 :
|
