1 : /******************************************************************************
2 : * $Id: gdaldem.cpp 18306 2009-12-15 18:57:11Z 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 18306 2009-12-15 18:57:11Z rouault $");
94 :
95 : #ifndef M_PI
96 : # define M_PI 3.1415926535897932384626433832795
97 : #endif
98 :
99 : /************************************************************************/
100 : /* Usage() */
101 : /************************************************************************/
102 :
103 0 : static void Usage()
104 :
105 : {
106 : printf( " Usage: \n"
107 : " - To generate a shaded relief map from any GDAL-supported elevation raster : \n\n"
108 : " gdaldem hillshade input_dem output_hillshade \n"
109 : " [-z ZFactor (default=1)] [-s scale* (default=1)] \n"
110 : " [-az Azimuth (default=315)] [-alt Altitude (default=45)]\n"
111 : " [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
112 : "\n"
113 : " - To generates a slope map from any GDAL-supported elevation raster :\n\n"
114 : " gdaldem slope input_dem output_slope_map \n"
115 : " [-p use percent slope (default=degrees)] [-s scale* (default=1)]\n"
116 : " [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
117 : "\n"
118 : " - To generate an aspect map from any GDAL-supported elevation raster\n"
119 : " Outputs a 32-bit float tiff with pixel values from 0-360 indicating azimuth :\n\n"
120 : " gdaldem aspect input_dem output_aspect_map \n"
121 : " [-trigonometric] [-zero_for_flat]\n"
122 : " [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
123 : "\n"
124 : " - To generate a color relief map from any GDAL-supported elevation raster\n"
125 : " gdaldem color-relief input_dem color_text_file output_color_relief_map\n"
126 : " [-alpha] [-exact_color_entry | -nearest_color_entry]\n"
127 : " [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
128 : " where color_text_file contains lines of the format \"elevation_value red green blue\"\n"
129 : "\n"
130 : " - To generate a Terrain Ruggedness Index (TRI) map from any GDAL-supported elevation raster\n"
131 : " gdaldem TRI input_dem output_TRI_map\n"
132 : " [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
133 : "\n"
134 : " - To generate a Topographic Position Index (TPI) map from any GDAL-supported elevation raster\n"
135 : " gdaldem TPI input_dem output_TPI_map\n"
136 : " [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
137 : "\n"
138 : " - To generate a roughness map from any GDAL-supported elevation raster\n"
139 : " gdaldem roughness input_dem output_roughness_map\n"
140 : " [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
141 : "\n"
142 : " Notes : \n"
143 : " Scale is the ratio of vertical units to horizontal\n"
144 0 : " for Feet:Latlong use scale=370400, for Meters:LatLong use scale=111120 \n\n");
145 0 : exit( 1 );
146 : }
147 :
148 : /************************************************************************/
149 : /* GDALGeneric3x3Processing() */
150 : /************************************************************************/
151 :
152 : typedef float (*GDALGeneric3x3ProcessingAlg) (float* pafWindow, float fDstNoDataValue, void* pData);
153 :
154 3 : CPLErr GDALGeneric3x3Processing ( GDALRasterBandH hSrcBand,
155 : GDALRasterBandH hDstBand,
156 : GDALGeneric3x3ProcessingAlg pfnAlg,
157 : void* pData,
158 : GDALProgressFunc pfnProgress,
159 : void * pProgressData)
160 : {
161 : CPLErr eErr;
162 : float *pafThreeLineWin; /* 3 line rotating source buffer */
163 : float *pafOutputBuf; /* 1 line destination buffer */
164 : int i, j;
165 :
166 : int bSrcHasNoData, bDstHasNoData;
167 3 : float fSrcNoDataValue = 0.0, fDstNoDataValue = 0.0;
168 :
169 3 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
170 3 : int nYSize = GDALGetRasterBandYSize(hSrcBand);
171 :
172 3 : if (pfnProgress == NULL)
173 0 : pfnProgress = GDALDummyProgress;
174 :
175 : /* -------------------------------------------------------------------- */
176 : /* Initialize progress counter. */
177 : /* -------------------------------------------------------------------- */
178 3 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
179 : {
180 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
181 0 : return CE_Failure;
182 : }
183 :
184 3 : pafOutputBuf = (float *) CPLMalloc(sizeof(float)*nXSize);
185 3 : pafThreeLineWin = (float *) CPLMalloc(3*sizeof(float)*(nXSize+1));
186 :
187 3 : fSrcNoDataValue = (float) GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
188 3 : fDstNoDataValue = (float) GDALGetRasterNoDataValue(hDstBand, &bDstHasNoData);
189 3 : if (!bDstHasNoData)
190 0 : fDstNoDataValue = 0.0;
191 :
192 : // Move a 3x3 pafWindow over each cell
193 : // (where the cell in question is #4)
194 : //
195 : // 0 1 2
196 : // 3 4 5
197 : // 6 7 8
198 :
199 : /* Preload the first 2 lines */
200 9 : for ( i = 0; i < 2 && i < nYSize; i++)
201 : {
202 : GDALRasterIO( hSrcBand,
203 : GF_Read,
204 : 0, i,
205 : nXSize, 1,
206 : pafThreeLineWin + i * nXSize,
207 : nXSize, 1,
208 : GDT_Float32,
209 6 : 0, 0);
210 : }
211 :
212 : // Exclude the edges
213 366 : for (j = 0; j < nXSize; j++)
214 : {
215 363 : pafOutputBuf[j] = fDstNoDataValue;
216 : }
217 : GDALRasterIO(hDstBand, GF_Write,
218 : 0, 0, nXSize, 1,
219 3 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
220 :
221 3 : if (nYSize > 1)
222 : {
223 : GDALRasterIO(hDstBand, GF_Write,
224 : 0, nYSize - 1, nXSize, 1,
225 3 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
226 : }
227 :
228 3 : int nLine1Off = 0*nXSize;
229 3 : int nLine2Off = 1*nXSize;
230 3 : int nLine3Off = 2*nXSize;
231 :
232 360 : for ( i = 1; i < nYSize-1; i++)
233 : {
234 : /* Read third line of the line buffer */
235 : eErr = GDALRasterIO( hSrcBand,
236 : GF_Read,
237 : 0, i+1,
238 : nXSize, 1,
239 : pafThreeLineWin + nLine3Off,
240 : nXSize, 1,
241 : GDT_Float32,
242 357 : 0, 0);
243 357 : if (eErr != CE_None)
244 0 : goto end;
245 :
246 : // Exclude the edges
247 357 : pafOutputBuf[0] = fDstNoDataValue;
248 357 : if (nXSize > 1)
249 357 : pafOutputBuf[nXSize - 1] = fDstNoDataValue;
250 :
251 42840 : for (j = 1; j < nXSize - 1; j++)
252 : {
253 : float afWin[9];
254 42483 : afWin[0] = pafThreeLineWin[nLine1Off + j-1];
255 42483 : afWin[1] = pafThreeLineWin[nLine1Off + j];
256 42483 : afWin[2] = pafThreeLineWin[nLine1Off + j+1];
257 42483 : afWin[3] = pafThreeLineWin[nLine2Off + j-1];
258 42483 : afWin[4] = pafThreeLineWin[nLine2Off + j];
259 42483 : afWin[5] = pafThreeLineWin[nLine2Off + j+1];
260 42483 : afWin[6] = pafThreeLineWin[nLine3Off + j-1];
261 42483 : afWin[7] = pafThreeLineWin[nLine3Off + j];
262 42483 : afWin[8] = pafThreeLineWin[nLine3Off + j+1];
263 :
264 424830 : if (bSrcHasNoData && (
265 42483 : afWin[0] == fSrcNoDataValue ||
266 42483 : afWin[1] == fSrcNoDataValue ||
267 42483 : afWin[2] == fSrcNoDataValue ||
268 42483 : afWin[3] == fSrcNoDataValue ||
269 42483 : afWin[4] == fSrcNoDataValue ||
270 42483 : afWin[5] == fSrcNoDataValue ||
271 42483 : afWin[6] == fSrcNoDataValue ||
272 42483 : afWin[7] == fSrcNoDataValue ||
273 42483 : afWin[8] == fSrcNoDataValue))
274 : {
275 : // We have nulls so write nullValue and move on
276 0 : pafOutputBuf[j] = fDstNoDataValue;
277 : }
278 : else
279 : {
280 : // We have a valid 3x3 window.
281 42483 : pafOutputBuf[j] = pfnAlg(afWin, fDstNoDataValue, pData);
282 : }
283 : }
284 :
285 : /* -----------------------------------------
286 : * Write Line to Raster
287 : */
288 : eErr = GDALRasterIO(hDstBand, GF_Write, 0, i, nXSize, 1,
289 357 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
290 357 : if (eErr != CE_None)
291 0 : goto end;
292 :
293 357 : if( !pfnProgress( 1.0 * (i+1) / nYSize, NULL, pProgressData ) )
294 : {
295 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
296 0 : eErr = CE_Failure;
297 0 : goto end;
298 : }
299 :
300 357 : int nTemp = nLine1Off;
301 357 : nLine1Off = nLine2Off;
302 357 : nLine2Off = nLine3Off;
303 357 : nLine3Off = nTemp;
304 : }
305 :
306 3 : pfnProgress( 1.0, NULL, pProgressData );
307 3 : eErr = CE_None;
308 :
309 : end:
310 3 : CPLFree(pafOutputBuf);
311 3 : CPLFree(pafThreeLineWin);
312 :
313 3 : return eErr;
314 : }
315 :
316 :
317 : /************************************************************************/
318 : /* GDALHillshade() */
319 : /************************************************************************/
320 :
321 : typedef struct
322 : {
323 : double nsres;
324 : double ewres;
325 : double sin_altRadians;
326 : double cos_altRadians_mul_z_scale_factor;
327 : double azRadians;
328 : double square_z_scale_factor;
329 : } GDALHillshadeAlgData;
330 :
331 : /* Unoptimized formulas are :
332 : x = psData->z*((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
333 : (afWin[2] + afWin[5] + afWin[5] + afWin[8])) /
334 : (8.0 * psData->ewres * psData->scale);
335 :
336 : y = psData->z*((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
337 : (afWin[0] + afWin[1] + afWin[1] + afWin[2])) /
338 : (8.0 * psData->nsres * psData->scale);
339 :
340 : slope = M_PI / 2 - atan(sqrt(x*x + y*y));
341 :
342 : aspect = atan2(x,y);
343 :
344 : cang = sin(alt * degreesToRadians) * sin(slope) +
345 : cos(alt * degreesToRadians) * cos(slope) *
346 : cos(az * degreesToRadians - M_PI/2 - aspect);
347 : */
348 :
349 28322 : float GDALHillshadeAlg (float* afWin, float fDstNoDataValue, void* pData)
350 : {
351 28322 : GDALHillshadeAlgData* psData = (GDALHillshadeAlgData*)pData;
352 : double x, y, aspect, xx_plus_yy, cang;
353 :
354 : // First Slope ...
355 84966 : x = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
356 84966 : (afWin[2] + afWin[5] + afWin[5] + afWin[8])) / psData->ewres;
357 :
358 84966 : y = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
359 84966 : (afWin[0] + afWin[1] + afWin[1] + afWin[2])) / psData->nsres;
360 :
361 28322 : xx_plus_yy = x * x + y * y;
362 :
363 : // ... then aspect...
364 28322 : aspect = atan2(x,y);
365 :
366 : // ... then the shade value
367 : cang = (psData->sin_altRadians -
368 : psData->cos_altRadians_mul_z_scale_factor * sqrt(xx_plus_yy) *
369 : sin(aspect - psData->azRadians)) /
370 28322 : sqrt(1 + psData->square_z_scale_factor * xx_plus_yy);
371 :
372 28322 : if (cang <= 0.0)
373 202 : cang = 1.0;
374 : else
375 28120 : cang = 1.0 + (254.0 * cang);
376 :
377 28322 : return cang;
378 : }
379 :
380 2 : void* GDALCreateHillshadeData(double* adfGeoTransform,
381 : double z,
382 : double scale,
383 : double alt,
384 : double az)
385 : {
386 : GDALHillshadeAlgData* pData =
387 2 : (GDALHillshadeAlgData*)CPLMalloc(sizeof(GDALHillshadeAlgData));
388 :
389 2 : const double degreesToRadians = M_PI / 180.0;
390 2 : pData->nsres = adfGeoTransform[5];
391 2 : pData->ewres = adfGeoTransform[1];
392 2 : pData->sin_altRadians = sin(alt * degreesToRadians);
393 2 : pData->azRadians = az * degreesToRadians;
394 2 : double z_scale_factor = z / (8 * scale);
395 : pData->cos_altRadians_mul_z_scale_factor =
396 2 : cos(alt * degreesToRadians) * z_scale_factor;
397 2 : pData->square_z_scale_factor = z_scale_factor * z_scale_factor;
398 2 : return pData;
399 : }
400 :
401 : /************************************************************************/
402 : /* GDALSlope() */
403 : /************************************************************************/
404 :
405 : typedef struct
406 : {
407 : double nsres;
408 : double ewres;
409 : double scale;
410 : int slopeFormat;
411 : } GDALSlopeAlgData;
412 :
413 14161 : float GDALSlopeAlg (float* afWin, float fDstNoDataValue, void* pData)
414 : {
415 14161 : const double radiansToDegrees = 180.0 / M_PI;
416 14161 : GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
417 : double dx, dy, key;
418 :
419 42483 : dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
420 42483 : (afWin[2] + afWin[5] + afWin[5] + afWin[8]))/psData->ewres;
421 :
422 42483 : dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
423 42483 : (afWin[0] + afWin[1] + afWin[1] + afWin[2]))/psData->nsres;
424 :
425 14161 : key = (dx * dx + dy * dy);
426 :
427 14161 : if (psData->slopeFormat == 1)
428 14161 : return atan(sqrt(key) / (8*psData->scale)) * radiansToDegrees;
429 : else
430 0 : return 100*(sqrt(key) / (8*psData->scale));
431 : }
432 :
433 1 : void* GDALCreateSlopeData(double* adfGeoTransform,
434 : double scale,
435 : int slopeFormat)
436 : {
437 : GDALSlopeAlgData* pData =
438 1 : (GDALSlopeAlgData*)CPLMalloc(sizeof(GDALSlopeAlgData));
439 :
440 1 : pData->nsres = adfGeoTransform[5];
441 1 : pData->ewres = adfGeoTransform[1];
442 1 : pData->scale = scale;
443 1 : pData->slopeFormat = slopeFormat;
444 1 : return pData;
445 : }
446 :
447 : /************************************************************************/
448 : /* GDALAspect() */
449 : /************************************************************************/
450 :
451 : typedef struct
452 : {
453 : int bAngleAsAzimuth;
454 : } GDALAspectAlgData;
455 :
456 14161 : float GDALAspectAlg (float* afWin, float fDstNoDataValue, void* pData)
457 : {
458 14161 : const double degreesToRadians = M_PI / 180.0;
459 14161 : GDALAspectAlgData* psData = (GDALAspectAlgData*)pData;
460 : double dx, dy;
461 : float aspect;
462 :
463 42483 : dx = ((afWin[2] + afWin[5] + afWin[5] + afWin[8]) -
464 42483 : (afWin[0] + afWin[3] + afWin[3] + afWin[6]));
465 :
466 42483 : dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
467 42483 : (afWin[0] + afWin[1] + afWin[1] + afWin[2]));
468 :
469 14161 : aspect = atan2(dy,-dx) / degreesToRadians;
470 :
471 18344 : if (dx == 0 && dy == 0)
472 : {
473 : /* Flat area */
474 4183 : aspect = fDstNoDataValue;
475 : }
476 9978 : else if ( psData->bAngleAsAzimuth )
477 : {
478 9978 : if (aspect > 90.0)
479 1243 : aspect = 450.0 - aspect;
480 : else
481 8735 : aspect = 90.0 - aspect;
482 : }
483 : else
484 : {
485 0 : if (aspect < 0)
486 0 : aspect += 360.0;
487 : }
488 :
489 14161 : if (aspect == 360.0)
490 0 : aspect = 0.0;
491 :
492 14161 : return aspect;
493 : }
494 :
495 1 : void* GDALCreateAspectData(int bAngleAsAzimuth)
496 : {
497 : GDALAspectAlgData* pData =
498 1 : (GDALAspectAlgData*)CPLMalloc(sizeof(GDALAspectAlgData));
499 :
500 1 : pData->bAngleAsAzimuth = bAngleAsAzimuth;
501 1 : return pData;
502 : }
503 :
504 : /************************************************************************/
505 : /* GDALColorRelief() */
506 : /************************************************************************/
507 :
508 : typedef struct
509 : {
510 : double dfVal;
511 : int nR;
512 : int nG;
513 : int nB;
514 : int nA;
515 : } ColorAssociation;
516 :
517 55 : static int GDALColorReliefSortColors(const void* pA, const void* pB)
518 : {
519 55 : ColorAssociation* pC1 = (ColorAssociation*)pA;
520 55 : ColorAssociation* pC2 = (ColorAssociation*)pB;
521 : return (pC1->dfVal < pC2->dfVal) ? -1 :
522 55 : (pC1->dfVal == pC2->dfVal) ? 0 : 1;
523 : }
524 :
525 : typedef enum
526 : {
527 : COLOR_SELECTION_INTERPOLATE,
528 : COLOR_SELECTION_NEAREST_ENTRY,
529 : COLOR_SELECTION_EXACT_ENTRY
530 : } ColorSelectionMode;
531 :
532 131769 : static int GDALColorReliefGetRGBA (ColorAssociation* pasColorAssociation,
533 : int nColorAssociation,
534 : double dfVal,
535 : ColorSelectionMode eColorSelectionMode,
536 : int* pnR,
537 : int* pnG,
538 : int* pnB,
539 : int* pnA)
540 : {
541 : int i;
542 131769 : int lower = 0;
543 131769 : int upper = nColorAssociation - 1;
544 : int mid;
545 :
546 : /* Find the index of the first element in the LUT input array that */
547 : /* is not smaller than the dfVal value. */
548 288783 : while(TRUE)
549 : {
550 420552 : mid = (lower + upper) / 2;
551 420552 : if (upper - lower <= 1)
552 : {
553 131769 : if (dfVal < pasColorAssociation[lower].dfVal)
554 0 : i = lower;
555 131769 : else if (dfVal < pasColorAssociation[upper].dfVal)
556 89685 : i = upper;
557 : else
558 42084 : i = upper + 1;
559 : break;
560 : }
561 288783 : else if (pasColorAssociation[mid].dfVal >= dfVal)
562 : {
563 173367 : upper = mid;
564 : }
565 : else
566 : {
567 115416 : lower = mid;
568 : }
569 : }
570 :
571 131769 : if (i == 0)
572 : {
573 0 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
574 : pasColorAssociation[0].dfVal != dfVal)
575 : {
576 0 : *pnR = 0;
577 0 : *pnG = 0;
578 0 : *pnB = 0;
579 0 : *pnA = 0;
580 0 : return FALSE;
581 : }
582 : else
583 : {
584 0 : *pnR = pasColorAssociation[0].nR;
585 0 : *pnG = pasColorAssociation[0].nG;
586 0 : *pnB = pasColorAssociation[0].nB;
587 0 : *pnA = pasColorAssociation[0].nA;
588 0 : return TRUE;
589 : }
590 : }
591 131769 : else if (i == nColorAssociation)
592 : {
593 0 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
594 0 : pasColorAssociation[i-1].dfVal != dfVal)
595 : {
596 0 : *pnR = 0;
597 0 : *pnG = 0;
598 0 : *pnB = 0;
599 0 : *pnA = 0;
600 0 : return FALSE;
601 : }
602 : else
603 : {
604 0 : *pnR = pasColorAssociation[i-1].nR;
605 0 : *pnG = pasColorAssociation[i-1].nG;
606 0 : *pnB = pasColorAssociation[i-1].nB;
607 0 : *pnA = pasColorAssociation[i-1].nA;
608 0 : return TRUE;
609 : }
610 : }
611 : else
612 : {
613 131769 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
614 0 : pasColorAssociation[i-1].dfVal != dfVal)
615 : {
616 0 : *pnR = 0;
617 0 : *pnG = 0;
618 0 : *pnB = 0;
619 0 : *pnA = 0;
620 0 : return FALSE;
621 : }
622 :
623 146410 : if (eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY &&
624 14641 : pasColorAssociation[i-1].dfVal != dfVal)
625 : {
626 : int index;
627 19930 : if (dfVal - pasColorAssociation[i-1].dfVal <
628 9965 : pasColorAssociation[i].dfVal - dfVal)
629 6716 : index = i -1;
630 : else
631 3249 : index = i;
632 :
633 9965 : *pnR = pasColorAssociation[index].nR;
634 9965 : *pnG = pasColorAssociation[index].nG;
635 9965 : *pnB = pasColorAssociation[index].nB;
636 9965 : *pnA = pasColorAssociation[index].nA;
637 9965 : return TRUE;
638 : }
639 :
640 121804 : double dfRatio = (dfVal - pasColorAssociation[i-1].dfVal) /
641 121804 : (pasColorAssociation[i].dfVal - pasColorAssociation[i-1].dfVal);
642 121804 : *pnR = (int)(0.45 + pasColorAssociation[i-1].nR + dfRatio *
643 121804 : (pasColorAssociation[i].nR - pasColorAssociation[i-1].nR));
644 121804 : if (*pnR < 0) *pnR = 0;
645 121804 : else if (*pnR > 255) *pnR = 255;
646 121804 : *pnG = (int)(0.45 + pasColorAssociation[i-1].nG + dfRatio *
647 121804 : (pasColorAssociation[i].nG - pasColorAssociation[i-1].nG));
648 121804 : if (*pnG < 0) *pnG = 0;
649 121804 : else if (*pnG > 255) *pnG = 255;
650 121804 : *pnB = (int)(0.45 + pasColorAssociation[i-1].nB + dfRatio *
651 121804 : (pasColorAssociation[i].nB - pasColorAssociation[i-1].nB));
652 121804 : if (*pnB < 0) *pnB = 0;
653 121804 : else if (*pnB > 255) *pnB = 255;
654 121804 : *pnA = (int)(0.45 + pasColorAssociation[i-1].nA + dfRatio *
655 121804 : (pasColorAssociation[i].nA - pasColorAssociation[i-1].nA));
656 121804 : if (*pnA < 0) *pnA = 0;
657 121804 : else if (*pnA > 255) *pnA = 255;
658 :
659 121804 : return TRUE;
660 : }
661 : }
662 :
663 : /* dfPct : percentage between 0 and 1 */
664 0 : static double GDALColorReliefGetAbsoluteValFromPct(GDALRasterBandH hSrcBand,
665 : double dfPct)
666 : {
667 : double dfMin, dfMax;
668 : int bSuccessMin, bSuccessMax;
669 0 : dfMin = GDALGetRasterMinimum(hSrcBand, &bSuccessMin);
670 0 : dfMax = GDALGetRasterMaximum(hSrcBand, &bSuccessMax);
671 0 : if (!bSuccessMin || !bSuccessMax)
672 : {
673 : double dfMean, dfStdDev;
674 0 : fprintf(stderr, "Computing source raster statistics...\n");
675 : GDALComputeRasterStatistics(hSrcBand, FALSE, &dfMin, &dfMax,
676 0 : &dfMean, &dfStdDev, NULL, NULL);
677 : }
678 0 : return dfMin + dfPct * (dfMax - dfMin);
679 : }
680 :
681 : typedef struct
682 : {
683 : const char *name;
684 : float r, g, b;
685 : } NamedColor;
686 :
687 : static const NamedColor namedColors[] = {
688 : { "white", 1.00, 1.00, 1.00 },
689 : { "black", 0.00, 0.00, 0.00 },
690 : { "red", 1.00, 0.00, 0.00 },
691 : { "green", 0.00, 1.00, 0.00 },
692 : { "blue", 0.00, 0.00, 1.00 },
693 : { "yellow", 1.00, 1.00, 0.00 },
694 : { "magenta",1.00, 0.00, 1.00 },
695 : { "cyan", 0.00, 1.00, 1.00 },
696 : { "aqua", 0.00, 0.75, 0.75 },
697 : { "grey", 0.75, 0.75, 0.75 },
698 : { "gray", 0.75, 0.75, 0.75 },
699 : { "orange", 1.00, 0.50, 0.00 },
700 : { "brown", 0.75, 0.50, 0.25 },
701 : { "purple", 0.50, 0.00, 1.00 },
702 : { "violet", 0.50, 0.00, 1.00 },
703 : { "indigo", 0.00, 0.50, 1.00 },
704 : };
705 :
706 : static
707 0 : int GDALColorReliefFindNamedColor(const char *pszColorName, int *pnR, int *pnG, int *pnB)
708 : {
709 : unsigned int i;
710 :
711 0 : *pnR = *pnG = *pnB = 0;
712 0 : for (i = 0; i < sizeof(namedColors) / sizeof(namedColors[0]); i++)
713 : {
714 0 : if (EQUAL(pszColorName, namedColors[i].name))
715 : {
716 0 : *pnR = (int)(255. * namedColors[i].r);
717 0 : *pnG = (int)(255. * namedColors[i].g);
718 0 : *pnB = (int)(255. * namedColors[i].b);
719 0 : return TRUE;
720 : }
721 : }
722 0 : return FALSE;
723 : }
724 :
725 : static
726 5 : ColorAssociation* GDALColorReliefParseColorFile(GDALRasterBandH hSrcBand,
727 : const char* pszColorFilename,
728 : int* pnColors)
729 : {
730 5 : FILE* fpColorFile = VSIFOpen(pszColorFilename, "rt");
731 5 : if (fpColorFile == NULL)
732 : {
733 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszColorFilename);
734 0 : *pnColors = 0;
735 0 : return NULL;
736 : }
737 :
738 5 : ColorAssociation* pasColorAssociation = NULL;
739 5 : int nColorAssociation = 0;
740 :
741 5 : int bSrcHasNoData = FALSE;
742 5 : double dfSrcNoDataValue = GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
743 :
744 : const char* pszLine;
745 45 : while ((pszLine = CPLReadLine(fpColorFile)) != NULL)
746 : {
747 : char** papszFields = CSLTokenizeStringComplex(pszLine, " ,\t:",
748 35 : FALSE, FALSE );
749 : /* Skip comment lines */
750 35 : int nTokens = CSLCount(papszFields);
751 105 : if (nTokens >= 2 &&
752 35 : papszFields[0][0] != '#' &&
753 35 : papszFields[0][0] != '/')
754 : {
755 : pasColorAssociation =
756 : (ColorAssociation*)CPLRealloc(pasColorAssociation,
757 35 : (nColorAssociation + 1) * sizeof(ColorAssociation));
758 35 : if (EQUAL(papszFields[0], "nv") && bSrcHasNoData)
759 0 : pasColorAssociation[nColorAssociation].dfVal = dfSrcNoDataValue;
760 35 : else if (strlen(papszFields[0]) > 1 && papszFields[0][strlen(papszFields[0])-1] == '%')
761 : {
762 0 : double dfPct = atof(papszFields[0]) / 100.;
763 0 : if (dfPct < 0.0 || dfPct > 1.0)
764 : {
765 : CPLError(CE_Failure, CPLE_AppDefined,
766 0 : "Wrong value for a percentage : %s", papszFields[0]);
767 0 : CSLDestroy(papszFields);
768 0 : VSIFClose(fpColorFile);
769 0 : CPLFree(pasColorAssociation);
770 0 : *pnColors = 0;
771 0 : return NULL;
772 : }
773 0 : pasColorAssociation[nColorAssociation].dfVal =
774 0 : GDALColorReliefGetAbsoluteValFromPct(hSrcBand, dfPct);
775 : }
776 : else
777 35 : pasColorAssociation[nColorAssociation].dfVal = atof(papszFields[0]);
778 :
779 35 : if (nTokens >= 4)
780 : {
781 35 : pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]);
782 35 : pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]);
783 35 : pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]);
784 35 : pasColorAssociation[nColorAssociation].nA =
785 35 : (CSLCount(papszFields) >= 5 ) ? atoi(papszFields[4]) : 255;
786 : }
787 : else
788 : {
789 : int nR, nG, nB;
790 0 : if (!GDALColorReliefFindNamedColor(papszFields[1], &nR, &nG, &nB))
791 : {
792 : CPLError(CE_Failure, CPLE_AppDefined,
793 0 : "Unknown color : %s", papszFields[1]);
794 0 : CSLDestroy(papszFields);
795 0 : VSIFClose(fpColorFile);
796 0 : CPLFree(pasColorAssociation);
797 0 : *pnColors = 0;
798 0 : return NULL;
799 : }
800 0 : pasColorAssociation[nColorAssociation].nR = nR;
801 0 : pasColorAssociation[nColorAssociation].nG = nG;
802 0 : pasColorAssociation[nColorAssociation].nB = nB;
803 0 : pasColorAssociation[nColorAssociation].nA =
804 0 : (CSLCount(papszFields) >= 3 ) ? atoi(papszFields[2]) : 255;
805 : }
806 :
807 35 : nColorAssociation ++;
808 : }
809 35 : CSLDestroy(papszFields);
810 : }
811 5 : VSIFClose(fpColorFile);
812 :
813 5 : if (nColorAssociation == 0)
814 : {
815 : CPLError(CE_Failure, CPLE_AppDefined,
816 0 : "No color association found in %s", pszColorFilename);
817 0 : *pnColors = 0;
818 0 : return NULL;
819 : }
820 :
821 : qsort(pasColorAssociation, nColorAssociation,
822 5 : sizeof(ColorAssociation), GDALColorReliefSortColors);
823 :
824 5 : *pnColors = nColorAssociation;
825 5 : return pasColorAssociation;
826 : }
827 :
828 : static
829 5 : GByte* GDALColorReliefPrecompute(GDALRasterBandH hSrcBand,
830 : ColorAssociation* pasColorAssociation,
831 : int nColorAssociation,
832 : ColorSelectionMode eColorSelectionMode,
833 : int* pnIndexOffset)
834 : {
835 5 : GDALDataType eDT = GDALGetRasterDataType(hSrcBand);
836 5 : GByte* pabyPrecomputed = NULL;
837 5 : int nIndexOffset = (eDT == GDT_Int16) ? 32768 : 0;
838 5 : *pnIndexOffset = nIndexOffset;
839 5 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
840 5 : int nYSize = GDALGetRasterBandXSize(hSrcBand);
841 5 : if (eDT == GDT_Byte ||
842 : ((eDT == GDT_Int16 || eDT == GDT_UInt16) && nXSize * nYSize > 65536))
843 : {
844 0 : int iMax = (eDT == GDT_Byte) ? 256: 65536;
845 0 : pabyPrecomputed = (GByte*) VSIMalloc(4 * iMax);
846 0 : if (pabyPrecomputed)
847 : {
848 : int i;
849 0 : for(i=0;i<iMax;i++)
850 : {
851 : int nR, nG, nB, nA;
852 : GDALColorReliefGetRGBA (pasColorAssociation,
853 : nColorAssociation,
854 : i - nIndexOffset,
855 : eColorSelectionMode,
856 0 : &nR, &nG, &nB, &nA);
857 0 : pabyPrecomputed[4 * i] = (GByte) nR;
858 0 : pabyPrecomputed[4 * i + 1] = (GByte) nG;
859 0 : pabyPrecomputed[4 * i + 2] = (GByte) nB;
860 0 : pabyPrecomputed[4 * i + 3] = (GByte) nA;
861 : }
862 : }
863 : }
864 5 : return pabyPrecomputed;
865 : }
866 :
867 : /************************************************************************/
868 : /* ==================================================================== */
869 : /* GDALColorReliefDataset */
870 : /* ==================================================================== */
871 : /************************************************************************/
872 :
873 : class GDALColorReliefRasterBand;
874 :
875 : class GDALColorReliefDataset : public GDALDataset
876 : {
877 : friend class GDALColorReliefRasterBand;
878 :
879 : GDALDatasetH hSrcDS;
880 : GDALRasterBandH hSrcBand;
881 : int nColorAssociation;
882 : ColorAssociation* pasColorAssociation;
883 : ColorSelectionMode eColorSelectionMode;
884 : GByte* pabyPrecomputed;
885 : int nIndexOffset;
886 : float* pafSourceBuf;
887 : int* panSourceBuf;
888 : int nCurBlockXOff;
889 : int nCurBlockYOff;
890 :
891 : public:
892 : GDALColorReliefDataset(GDALDatasetH hSrcDS,
893 : GDALRasterBandH hSrcBand,
894 : const char* pszColorFilename,
895 : ColorSelectionMode eColorSelectionMode,
896 : int bAlpha);
897 : ~GDALColorReliefDataset();
898 :
899 : CPLErr GetGeoTransform( double * padfGeoTransform );
900 : const char *GetProjectionRef();
901 : };
902 :
903 : /************************************************************************/
904 : /* ==================================================================== */
905 : /* GDALColorReliefRasterBand */
906 : /* ==================================================================== */
907 : /************************************************************************/
908 :
909 : class GDALColorReliefRasterBand : public GDALRasterBand
910 12 : {
911 : friend class GDALColorReliefDataset;
912 :
913 :
914 : public:
915 : GDALColorReliefRasterBand( GDALColorReliefDataset *, int );
916 :
917 : virtual CPLErr IReadBlock( int, int, void * );
918 : virtual GDALColorInterp GetColorInterpretation();
919 : };
920 :
921 2 : GDALColorReliefDataset::GDALColorReliefDataset(
922 : GDALDatasetH hSrcDS,
923 : GDALRasterBandH hSrcBand,
924 : const char* pszColorFilename,
925 : ColorSelectionMode eColorSelectionMode,
926 2 : int bAlpha)
927 : {
928 2 : this->hSrcDS = hSrcDS;
929 2 : this->hSrcBand = hSrcBand;
930 2 : nColorAssociation = 0;
931 : pasColorAssociation =
932 : GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
933 2 : &nColorAssociation);
934 2 : this->eColorSelectionMode = eColorSelectionMode;
935 :
936 2 : nRasterXSize = GDALGetRasterXSize(hSrcDS);
937 2 : nRasterYSize = GDALGetRasterYSize(hSrcDS);
938 :
939 : int nBlockXSize, nBlockYSize;
940 2 : GDALGetBlockSize( hSrcBand, &nBlockXSize, &nBlockYSize);
941 :
942 2 : nIndexOffset = 0;
943 : pabyPrecomputed =
944 : GDALColorReliefPrecompute(hSrcBand,
945 : pasColorAssociation,
946 : nColorAssociation,
947 : eColorSelectionMode,
948 2 : &nIndexOffset);
949 :
950 : int i;
951 8 : for(i=0;i<((bAlpha) ? 4 : 3);i++)
952 : {
953 6 : SetBand(i + 1, new GDALColorReliefRasterBand(this, i+1));
954 : }
955 :
956 2 : pafSourceBuf = NULL;
957 2 : panSourceBuf = NULL;
958 2 : if (pabyPrecomputed)
959 0 : panSourceBuf = (int *) CPLMalloc(sizeof(int)*nBlockXSize*nBlockYSize);
960 : else
961 2 : pafSourceBuf = (float *) CPLMalloc(sizeof(float)*nBlockXSize*nBlockYSize);
962 2 : nCurBlockXOff = -1;
963 2 : nCurBlockYOff = -1;
964 2 : }
965 :
966 4 : GDALColorReliefDataset::~GDALColorReliefDataset()
967 : {
968 2 : CPLFree(pasColorAssociation);
969 2 : CPLFree(pabyPrecomputed);
970 2 : CPLFree(panSourceBuf);
971 2 : CPLFree(pafSourceBuf);
972 4 : }
973 :
974 2 : CPLErr GDALColorReliefDataset::GetGeoTransform( double * padfGeoTransform )
975 : {
976 2 : return GDALGetGeoTransform(hSrcDS, padfGeoTransform);
977 : }
978 :
979 2 : const char *GDALColorReliefDataset::GetProjectionRef()
980 : {
981 2 : return GDALGetProjectionRef(hSrcDS);
982 : }
983 :
984 6 : GDALColorReliefRasterBand::GDALColorReliefRasterBand(
985 6 : GDALColorReliefDataset * poDS, int nBand)
986 : {
987 6 : this->poDS = poDS;
988 6 : this->nBand = nBand;
989 6 : eDataType = GDT_Byte;
990 6 : GDALGetBlockSize( poDS->hSrcBand, &nBlockXSize, &nBlockYSize);
991 6 : }
992 :
993 387 : CPLErr GDALColorReliefRasterBand::IReadBlock( int nBlockXOff,
994 : int nBlockYOff,
995 : void *pImage )
996 : {
997 387 : GDALColorReliefDataset * poGDS = (GDALColorReliefDataset *) poDS;
998 : int nReqXSize, nReqYSize;
999 :
1000 387 : if ((nBlockXOff + 1) * nBlockXSize >= nRasterXSize)
1001 27 : nReqXSize = nRasterXSize - nBlockXOff * nBlockXSize;
1002 : else
1003 360 : nReqXSize = nBlockXSize;
1004 :
1005 387 : if ((nBlockYOff + 1) * nBlockYSize >= nRasterYSize)
1006 366 : nReqYSize = nRasterYSize - nBlockYOff * nBlockYSize;
1007 : else
1008 21 : nReqYSize = nBlockYSize;
1009 :
1010 387 : int nCount = nReqXSize * nReqYSize;
1011 387 : if ( poGDS->nCurBlockXOff != nBlockXOff ||
1012 : poGDS->nCurBlockYOff != nBlockYOff )
1013 : {
1014 371 : poGDS->nCurBlockXOff = nBlockXOff;
1015 371 : poGDS->nCurBlockYOff = nBlockYOff;
1016 :
1017 : CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1018 : GF_Read,
1019 : nBlockXOff * nBlockXSize,
1020 : nBlockYOff * nBlockYSize,
1021 : nReqXSize, nReqYSize,
1022 : (poGDS->panSourceBuf) ?
1023 : (void*) poGDS->panSourceBuf :
1024 : (void* )poGDS->pafSourceBuf,
1025 : nReqXSize, nReqYSize,
1026 : (poGDS->panSourceBuf) ? GDT_Int32 : GDT_Float32,
1027 371 : 0, 0);
1028 371 : if (eErr != CE_None)
1029 : {
1030 0 : memset(pImage, 0, nCount);
1031 0 : return eErr;
1032 : }
1033 : }
1034 :
1035 : int j;
1036 387 : if (poGDS->panSourceBuf)
1037 : {
1038 0 : for ( j = 0; j < nCount; j++)
1039 : {
1040 0 : int nIndex = poGDS->panSourceBuf[j] + poGDS->nIndexOffset;
1041 0 : ((GByte*)pImage)[j] = poGDS->pabyPrecomputed[4*nIndex + nBand-1];
1042 : }
1043 : }
1044 : else
1045 : {
1046 : int anComponents[4];
1047 88233 : for ( j = 0; j < nCount; j++)
1048 : {
1049 : GDALColorReliefGetRGBA (poGDS->pasColorAssociation,
1050 : poGDS->nColorAssociation,
1051 87846 : poGDS->pafSourceBuf[j],
1052 : poGDS->eColorSelectionMode,
1053 : &anComponents[0],
1054 : &anComponents[1],
1055 : &anComponents[2],
1056 175692 : &anComponents[3]);
1057 87846 : ((GByte*)pImage)[j] = (GByte) anComponents[nBand-1];
1058 : }
1059 : }
1060 :
1061 387 : return CE_None;
1062 : }
1063 :
1064 12 : GDALColorInterp GDALColorReliefRasterBand::GetColorInterpretation()
1065 : {
1066 12 : return (GDALColorInterp)(GCI_RedBand + nBand - 1);
1067 : }
1068 :
1069 :
1070 3 : CPLErr GDALColorRelief (GDALRasterBandH hSrcBand,
1071 : GDALRasterBandH hDstBand1,
1072 : GDALRasterBandH hDstBand2,
1073 : GDALRasterBandH hDstBand3,
1074 : GDALRasterBandH hDstBand4,
1075 : const char* pszColorFilename,
1076 : ColorSelectionMode eColorSelectionMode,
1077 : GDALProgressFunc pfnProgress,
1078 : void * pProgressData)
1079 : {
1080 : CPLErr eErr;
1081 :
1082 3 : if (hSrcBand == NULL || hDstBand1 == NULL || hDstBand2 == NULL ||
1083 : hDstBand3 == NULL)
1084 0 : return CE_Failure;
1085 :
1086 3 : int nColorAssociation = 0;
1087 : ColorAssociation* pasColorAssociation =
1088 : GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
1089 3 : &nColorAssociation);
1090 3 : if (pasColorAssociation == NULL)
1091 0 : return CE_Failure;
1092 :
1093 3 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
1094 3 : int nYSize = GDALGetRasterBandYSize(hSrcBand);
1095 :
1096 3 : if (pfnProgress == NULL)
1097 0 : pfnProgress = GDALDummyProgress;
1098 :
1099 3 : int nR = 0, nG = 0, nB = 0, nA = 0;
1100 :
1101 : /* -------------------------------------------------------------------- */
1102 : /* Precompute the map from values to RGBA quadruplets */
1103 : /* for GDT_Byte, GDT_Int16 or GDT_UInt16 */
1104 : /* -------------------------------------------------------------------- */
1105 3 : int nIndexOffset = 0;
1106 : GByte* pabyPrecomputed =
1107 : GDALColorReliefPrecompute(hSrcBand,
1108 : pasColorAssociation,
1109 : nColorAssociation,
1110 : eColorSelectionMode,
1111 3 : &nIndexOffset);
1112 :
1113 : /* -------------------------------------------------------------------- */
1114 : /* Initialize progress counter. */
1115 : /* -------------------------------------------------------------------- */
1116 :
1117 3 : float* pafSourceBuf = NULL;
1118 3 : int* panSourceBuf = NULL;
1119 3 : if (pabyPrecomputed)
1120 0 : panSourceBuf = (int *) CPLMalloc(sizeof(int)*nXSize);
1121 : else
1122 3 : pafSourceBuf = (float *) CPLMalloc(sizeof(float)*nXSize);
1123 3 : GByte* pabyDestBuf1 = (GByte*) CPLMalloc( 4 * nXSize );
1124 3 : GByte* pabyDestBuf2 = pabyDestBuf1 + nXSize;
1125 3 : GByte* pabyDestBuf3 = pabyDestBuf2 + nXSize;
1126 3 : GByte* pabyDestBuf4 = pabyDestBuf3 + nXSize;
1127 : int i, j;
1128 :
1129 3 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
1130 : {
1131 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1132 0 : eErr = CE_Failure;
1133 0 : goto end;
1134 : }
1135 :
1136 366 : for ( i = 0; i < nYSize; i++)
1137 : {
1138 : /* Read source buffer */
1139 : eErr = GDALRasterIO( hSrcBand,
1140 : GF_Read,
1141 : 0, i,
1142 : nXSize, 1,
1143 : (panSourceBuf) ? (void*) panSourceBuf : (void* )pafSourceBuf,
1144 : nXSize, 1,
1145 : (panSourceBuf) ? GDT_Int32 : GDT_Float32,
1146 363 : 0, 0);
1147 363 : if (eErr != CE_None)
1148 0 : goto end;
1149 :
1150 363 : if (panSourceBuf)
1151 : {
1152 0 : for ( j = 0; j < nXSize; j++)
1153 : {
1154 0 : int nIndex = panSourceBuf[j] + nIndexOffset;
1155 0 : pabyDestBuf1[j] = pabyPrecomputed[4 * nIndex];
1156 0 : pabyDestBuf2[j] = pabyPrecomputed[4 * nIndex + 1];
1157 0 : pabyDestBuf3[j] = pabyPrecomputed[4 * nIndex + 2];
1158 0 : pabyDestBuf4[j] = pabyPrecomputed[4 * nIndex + 3];
1159 : }
1160 : }
1161 : else
1162 : {
1163 44286 : for ( j = 0; j < nXSize; j++)
1164 : {
1165 : GDALColorReliefGetRGBA (pasColorAssociation,
1166 : nColorAssociation,
1167 43923 : pafSourceBuf[j],
1168 : eColorSelectionMode,
1169 : &nR,
1170 : &nG,
1171 : &nB,
1172 43923 : &nA);
1173 43923 : pabyDestBuf1[j] = (GByte) nR;
1174 43923 : pabyDestBuf2[j] = (GByte) nG;
1175 43923 : pabyDestBuf3[j] = (GByte) nB;
1176 43923 : pabyDestBuf4[j] = (GByte) nA;
1177 : }
1178 : }
1179 :
1180 : /* -----------------------------------------
1181 : * Write Line to Raster
1182 : */
1183 : eErr = GDALRasterIO(hDstBand1,
1184 : GF_Write,
1185 : 0, i, nXSize,
1186 363 : 1, pabyDestBuf1, nXSize, 1, GDT_Byte, 0, 0);
1187 363 : if (eErr != CE_None)
1188 0 : goto end;
1189 :
1190 : eErr = GDALRasterIO(hDstBand2,
1191 : GF_Write,
1192 : 0, i, nXSize,
1193 363 : 1, pabyDestBuf2, nXSize, 1, GDT_Byte, 0, 0);
1194 363 : if (eErr != CE_None)
1195 0 : goto end;
1196 :
1197 : eErr = GDALRasterIO(hDstBand3,
1198 : GF_Write,
1199 : 0, i, nXSize,
1200 363 : 1, pabyDestBuf3, nXSize, 1, GDT_Byte, 0, 0);
1201 363 : if (eErr != CE_None)
1202 0 : goto end;
1203 :
1204 363 : if (hDstBand4)
1205 : {
1206 : eErr = GDALRasterIO(hDstBand4,
1207 : GF_Write,
1208 : 0, i, nXSize,
1209 0 : 1, pabyDestBuf4, nXSize, 1, GDT_Byte, 0, 0);
1210 0 : if (eErr != CE_None)
1211 0 : goto end;
1212 : }
1213 :
1214 363 : if( !pfnProgress( 1.0 * (i+1) / nYSize, NULL, pProgressData ) )
1215 : {
1216 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1217 0 : eErr = CE_Failure;
1218 0 : goto end;
1219 : }
1220 : }
1221 3 : pfnProgress( 1.0, NULL, pProgressData );
1222 3 : eErr = CE_None;
1223 :
1224 : end:
1225 3 : VSIFree(pabyPrecomputed);
1226 3 : CPLFree(pafSourceBuf);
1227 3 : CPLFree(panSourceBuf);
1228 3 : CPLFree(pabyDestBuf1);
1229 3 : CPLFree(pasColorAssociation);
1230 :
1231 3 : return eErr;
1232 : }
1233 :
1234 : /************************************************************************/
1235 : /* GDALTRI() */
1236 : /************************************************************************/
1237 :
1238 0 : float GDALTRIAlg (float* afWin, float fDstNoDataValue, void* pData)
1239 : {
1240 : // Terrain Ruggedness is average difference in height
1241 0 : return (fabs(afWin[0]-afWin[4]) +
1242 0 : fabs(afWin[1]-afWin[4]) +
1243 0 : fabs(afWin[2]-afWin[4]) +
1244 0 : fabs(afWin[3]-afWin[4]) +
1245 0 : fabs(afWin[5]-afWin[4]) +
1246 0 : fabs(afWin[6]-afWin[4]) +
1247 0 : fabs(afWin[7]-afWin[4]) +
1248 0 : fabs(afWin[8]-afWin[4]))/8;
1249 : }
1250 :
1251 0 : CPLErr GDALTRI ( GDALRasterBandH hSrcBand,
1252 : GDALRasterBandH hDstBand,
1253 : GDALProgressFunc pfnProgress,
1254 : void * pProgressData)
1255 : {
1256 : return GDALGeneric3x3Processing(hSrcBand, hDstBand,
1257 : GDALTRIAlg, NULL,
1258 0 : pfnProgress, pProgressData);
1259 : }
1260 :
1261 : /************************************************************************/
1262 : /* GDALTPI() */
1263 : /************************************************************************/
1264 :
1265 0 : float GDALTPIAlg (float* afWin, float fDstNoDataValue, void* pData)
1266 : {
1267 : // Terrain Position is the difference between
1268 : // The central cell and the mean of the surrounding cells
1269 0 : return afWin[4] -
1270 0 : ((afWin[0]+
1271 0 : afWin[1]+
1272 0 : afWin[2]+
1273 0 : afWin[3]+
1274 0 : afWin[5]+
1275 0 : afWin[6]+
1276 0 : afWin[7]+
1277 0 : afWin[8])/8);
1278 : }
1279 :
1280 : /************************************************************************/
1281 : /* GDALRoughness() */
1282 : /************************************************************************/
1283 :
1284 0 : float GDALRoughnessAlg (float* afWin, float fDstNoDataValue, void* pData)
1285 : {
1286 : // Roughness is the largest difference
1287 : // between any two cells
1288 :
1289 0 : float pafRoughnessMin = afWin[0];
1290 0 : float pafRoughnessMax = afWin[0];
1291 :
1292 0 : for ( int k = 1; k < 9; k++)
1293 : {
1294 0 : if (afWin[k] > pafRoughnessMax)
1295 : {
1296 0 : pafRoughnessMax=afWin[k];
1297 : }
1298 0 : if (afWin[k] < pafRoughnessMin)
1299 : {
1300 0 : pafRoughnessMin=afWin[k];
1301 : }
1302 : }
1303 0 : return pafRoughnessMax - pafRoughnessMin;
1304 : }
1305 :
1306 : /************************************************************************/
1307 : /* ==================================================================== */
1308 : /* GDALGeneric3x3Dataset */
1309 : /* ==================================================================== */
1310 : /************************************************************************/
1311 :
1312 : class GDALGeneric3x3RasterBand;
1313 :
1314 : class GDALGeneric3x3Dataset : public GDALDataset
1315 : {
1316 : friend class GDALGeneric3x3RasterBand;
1317 :
1318 : GDALGeneric3x3ProcessingAlg pfnAlg;
1319 : void* pAlgData;
1320 : GDALDatasetH hSrcDS;
1321 : GDALRasterBandH hSrcBand;
1322 : float* apafSourceBuf[3];
1323 : int bDstHasNoData;
1324 : double dfDstNoDataValue;
1325 : int nCurLine;
1326 :
1327 : public:
1328 : GDALGeneric3x3Dataset(GDALDatasetH hSrcDS,
1329 : GDALRasterBandH hSrcBand,
1330 : GDALDataType eDstDataType,
1331 : int bDstHasNoData,
1332 : double dfDstNoDataValue,
1333 : GDALGeneric3x3ProcessingAlg pfnAlg,
1334 : void* pAlgData);
1335 : ~GDALGeneric3x3Dataset();
1336 :
1337 : CPLErr GetGeoTransform( double * padfGeoTransform );
1338 : const char *GetProjectionRef();
1339 : };
1340 :
1341 : /************************************************************************/
1342 : /* ==================================================================== */
1343 : /* GDALGeneric3x3RasterBand */
1344 : /* ==================================================================== */
1345 : /************************************************************************/
1346 :
1347 : class GDALGeneric3x3RasterBand : public GDALRasterBand
1348 2 : {
1349 : friend class GDALGeneric3x3Dataset;
1350 :
1351 :
1352 : public:
1353 : GDALGeneric3x3RasterBand( GDALGeneric3x3Dataset *poDS,
1354 : GDALDataType eDstDataType );
1355 :
1356 : virtual CPLErr IReadBlock( int, int, void * );
1357 : virtual double GetNoDataValue( int* pbHasNoData );
1358 : };
1359 :
1360 1 : GDALGeneric3x3Dataset::GDALGeneric3x3Dataset(
1361 : GDALDatasetH hSrcDS,
1362 : GDALRasterBandH hSrcBand,
1363 : GDALDataType eDstDataType,
1364 : int bDstHasNoData,
1365 : double dfDstNoDataValue,
1366 : GDALGeneric3x3ProcessingAlg pfnAlg,
1367 1 : void* pAlgData)
1368 : {
1369 1 : this->hSrcDS = hSrcDS;
1370 1 : this->hSrcBand = hSrcBand;
1371 1 : this->pfnAlg = pfnAlg;
1372 1 : this->pAlgData = pAlgData;
1373 1 : this->bDstHasNoData = bDstHasNoData;
1374 1 : this->dfDstNoDataValue = dfDstNoDataValue;
1375 :
1376 : CPLAssert(eDstDataType == GDT_Byte || eDstDataType == GDT_Float32);
1377 :
1378 1 : nRasterXSize = GDALGetRasterXSize(hSrcDS);
1379 1 : nRasterYSize = GDALGetRasterYSize(hSrcDS);
1380 :
1381 1 : SetBand(1, new GDALGeneric3x3RasterBand(this, eDstDataType));
1382 :
1383 1 : apafSourceBuf[0] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1384 1 : apafSourceBuf[1] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1385 1 : apafSourceBuf[2] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1386 :
1387 1 : nCurLine = -1;
1388 1 : }
1389 :
1390 2 : GDALGeneric3x3Dataset::~GDALGeneric3x3Dataset()
1391 : {
1392 1 : CPLFree(apafSourceBuf[0]);
1393 1 : CPLFree(apafSourceBuf[1]);
1394 1 : CPLFree(apafSourceBuf[2]);
1395 2 : }
1396 :
1397 1 : CPLErr GDALGeneric3x3Dataset::GetGeoTransform( double * padfGeoTransform )
1398 : {
1399 1 : return GDALGetGeoTransform(hSrcDS, padfGeoTransform);
1400 : }
1401 :
1402 1 : const char *GDALGeneric3x3Dataset::GetProjectionRef()
1403 : {
1404 1 : return GDALGetProjectionRef(hSrcDS);
1405 : }
1406 :
1407 1 : GDALGeneric3x3RasterBand::GDALGeneric3x3RasterBand(GDALGeneric3x3Dataset *poDS,
1408 1 : GDALDataType eDstDataType)
1409 : {
1410 1 : this->poDS = poDS;
1411 1 : this->nBand = 1;
1412 1 : eDataType = eDstDataType;
1413 1 : nBlockXSize = poDS->GetRasterXSize();
1414 1 : nBlockYSize = 1;
1415 1 : }
1416 :
1417 121 : CPLErr GDALGeneric3x3RasterBand::IReadBlock( int nBlockXOff,
1418 : int nBlockYOff,
1419 : void *pImage )
1420 : {
1421 121 : GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
1422 :
1423 121 : if ( nBlockYOff == 0 || nBlockYOff == nRasterYSize - 1)
1424 : {
1425 : int j;
1426 2 : if (eDataType == GDT_Byte)
1427 : {
1428 244 : for(j=0;j<nBlockXSize;j++)
1429 242 : ((GByte*)pImage)[j] = (GByte) poGDS->dfDstNoDataValue;
1430 : }
1431 : else
1432 : {
1433 0 : for(j=0;j<nBlockXSize;j++)
1434 0 : ((float*)pImage)[j] = (float) poGDS->dfDstNoDataValue;
1435 : }
1436 :
1437 2 : return CE_None;
1438 : }
1439 :
1440 119 : if ( poGDS->nCurLine != nBlockYOff )
1441 : {
1442 119 : if (nBlockYOff != 1 &&
1443 : poGDS->nCurLine == nBlockYOff + 1)
1444 : {
1445 0 : float* pafTmp = poGDS->apafSourceBuf[0];
1446 0 : poGDS->apafSourceBuf[0] = poGDS->apafSourceBuf[1];
1447 0 : poGDS->apafSourceBuf[1] = poGDS->apafSourceBuf[2];
1448 0 : poGDS->apafSourceBuf[2] = pafTmp;
1449 :
1450 0 : poGDS->nCurLine = nBlockYOff;
1451 :
1452 : CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1453 : GF_Read,
1454 : 0, nBlockYOff + 1, nBlockXSize, 1,
1455 0 : poGDS->apafSourceBuf[2],
1456 : nBlockXSize, 1,
1457 : GDT_Float32,
1458 0 : 0, 0);
1459 :
1460 0 : if (eErr != CE_None)
1461 : {
1462 0 : memset(pImage, 0, nBlockXSize * sizeof(float));
1463 0 : return eErr;
1464 : }
1465 : }
1466 : else
1467 : {
1468 119 : poGDS->nCurLine = nBlockYOff;
1469 : int i;
1470 476 : for(i=0;i<3;i++)
1471 : {
1472 : CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1473 : GF_Read,
1474 : 0, nBlockYOff + i - 1, nBlockXSize, 1,
1475 357 : poGDS->apafSourceBuf[i],
1476 : nBlockXSize, 1,
1477 : GDT_Float32,
1478 714 : 0, 0);
1479 357 : if (eErr != CE_None)
1480 : {
1481 0 : memset(pImage, 0, nBlockXSize * sizeof(float));
1482 0 : return eErr;
1483 : }
1484 : }
1485 : }
1486 : }
1487 :
1488 : int j;
1489 :
1490 119 : if (eDataType == GDT_Byte)
1491 : {
1492 119 : ((GByte*)pImage)[0] = (GByte) poGDS->dfDstNoDataValue;
1493 119 : if (nBlockXSize > 1)
1494 119 : ((GByte*)pImage)[nBlockXSize - 1] = (GByte) poGDS->dfDstNoDataValue;
1495 : }
1496 : else
1497 : {
1498 0 : ((float*)pImage)[0] = (float) poGDS->dfDstNoDataValue;
1499 0 : if (nBlockXSize > 1)
1500 0 : ((float*)pImage)[nBlockXSize - 1] = (float) poGDS->dfDstNoDataValue;
1501 : }
1502 :
1503 : int bSrcHasNoData;
1504 : double dfSrcNoDataValue = GDALGetRasterNoDataValue(poGDS->hSrcBand,
1505 119 : &bSrcHasNoData);
1506 :
1507 14280 : for(j=1;j<nBlockXSize - 1;j++)
1508 : {
1509 : float afWin[9];
1510 14161 : afWin[0] = poGDS->apafSourceBuf[0][j-1];
1511 14161 : afWin[1] = poGDS->apafSourceBuf[0][j];
1512 14161 : afWin[2] = poGDS->apafSourceBuf[0][j+1];
1513 14161 : afWin[3] = poGDS->apafSourceBuf[1][j-1];
1514 14161 : afWin[4] = poGDS->apafSourceBuf[1][j];
1515 14161 : afWin[5] = poGDS->apafSourceBuf[1][j+1];
1516 14161 : afWin[6] = poGDS->apafSourceBuf[2][j-1];
1517 14161 : afWin[7] = poGDS->apafSourceBuf[2][j];
1518 14161 : afWin[8] = poGDS->apafSourceBuf[2][j+1];
1519 :
1520 : // Check if afWin has null value
1521 14161 : int bContainsNull = FALSE;
1522 14161 : if (bSrcHasNoData)
1523 : {
1524 141610 : for ( int n = 0; n <= 8; n++)
1525 : {
1526 127449 : if(afWin[n] == dfSrcNoDataValue)
1527 : {
1528 0 : bContainsNull = TRUE;
1529 0 : break;
1530 : }
1531 : }
1532 : }
1533 :
1534 14161 : if (bContainsNull)
1535 : {
1536 : // We have nulls so write nullValue and move on
1537 0 : if (eDataType == GDT_Byte)
1538 0 : ((GByte*)pImage)[j] = (GByte) poGDS->dfDstNoDataValue;
1539 : else
1540 0 : ((float*)pImage)[j] = (float) poGDS->dfDstNoDataValue;
1541 : } else {
1542 : // We have a valid 3x3 window.
1543 14161 : if (eDataType == GDT_Byte)
1544 14161 : ((GByte*)pImage)[j] = (GByte) (poGDS->pfnAlg(
1545 14161 : afWin, poGDS->dfDstNoDataValue, poGDS->pAlgData) + 0.5);
1546 : else
1547 0 : ((float*)pImage)[j] = (float) poGDS->pfnAlg(
1548 0 : afWin, poGDS->dfDstNoDataValue, poGDS->pAlgData);
1549 : }
1550 : }
1551 :
1552 119 : return CE_None;
1553 : }
1554 :
1555 4 : double GDALGeneric3x3RasterBand::GetNoDataValue( int* pbHasNoData )
1556 : {
1557 4 : GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
1558 4 : if (pbHasNoData)
1559 3 : *pbHasNoData = poGDS->bDstHasNoData;
1560 4 : return poGDS->dfDstNoDataValue;
1561 : }
1562 :
1563 : /************************************************************************/
1564 : /* main() */
1565 : /************************************************************************/
1566 :
1567 : typedef enum
1568 : {
1569 : HILL_SHADE,
1570 : SLOPE,
1571 : ASPECT,
1572 : COLOR_RELIEF,
1573 : TRI,
1574 : TPI,
1575 : ROUGHNESS
1576 : } Algorithm;
1577 :
1578 10 : int main( int argc, char ** argv )
1579 :
1580 : {
1581 : Algorithm eUtilityMode;
1582 10 : double z = 1.0;
1583 10 : double scale = 1.0;
1584 10 : double az = 315.0;
1585 10 : double alt = 45.0;
1586 : // 0 = 'percent' or 1 = 'degrees'
1587 10 : int slopeFormat = 1;
1588 10 : int bAddAlpha = FALSE;
1589 10 : int bZeroForFlat = FALSE;
1590 10 : int bAngleAsAzimuth = TRUE;
1591 10 : ColorSelectionMode eColorSelectionMode = COLOR_SELECTION_INTERPOLATE;
1592 :
1593 10 : int nBand = 1;
1594 : double adfGeoTransform[6];
1595 :
1596 10 : const char *pszSrcFilename = NULL;
1597 10 : const char *pszDstFilename = NULL;
1598 10 : const char *pszColorFilename = NULL;
1599 10 : const char *pszFormat = "GTiff";
1600 10 : char **papszCreateOptions = NULL;
1601 :
1602 10 : GDALDatasetH hSrcDataset = NULL;
1603 10 : GDALDatasetH hDstDataset = NULL;
1604 10 : GDALRasterBandH hSrcBand = NULL;
1605 10 : GDALRasterBandH hDstBand = NULL;
1606 10 : GDALDriverH hDriver = NULL;
1607 :
1608 10 : GDALProgressFunc pfnProgress = GDALTermProgress;
1609 :
1610 10 : int nXSize = 0;
1611 10 : int nYSize = 0;
1612 :
1613 : /* Check strict compilation and runtime library version as we use C++ API */
1614 10 : if (! GDAL_CHECK_VERSION(argv[0]))
1615 0 : exit(1);
1616 :
1617 10 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
1618 10 : if( argc < 2 )
1619 : {
1620 0 : fprintf(stderr, "Not enough arguments\n");
1621 0 : Usage();
1622 0 : exit( 1 );
1623 : }
1624 :
1625 10 : if( EQUAL(argv[1], "--utility_version") || EQUAL(argv[1], "--utility-version") )
1626 : {
1627 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
1628 1 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
1629 1 : return 0;
1630 : }
1631 11 : else if ( EQUAL(argv[1], "shade") || EQUAL(argv[1], "hillshade") )
1632 : {
1633 2 : eUtilityMode = HILL_SHADE;
1634 : }
1635 7 : else if ( EQUAL(argv[1], "slope") )
1636 : {
1637 1 : eUtilityMode = SLOPE;
1638 : }
1639 6 : else if ( EQUAL(argv[1], "aspect") )
1640 : {
1641 1 : eUtilityMode = ASPECT;
1642 : }
1643 5 : else if ( EQUAL(argv[1], "color-relief") )
1644 : {
1645 5 : eUtilityMode = COLOR_RELIEF;
1646 : }
1647 0 : else if ( EQUAL(argv[1], "TRI") )
1648 : {
1649 0 : eUtilityMode = TRI;
1650 : }
1651 0 : else if ( EQUAL(argv[1], "TPI") )
1652 : {
1653 0 : eUtilityMode = TPI;
1654 : }
1655 0 : else if ( EQUAL(argv[1], "roughness") )
1656 : {
1657 0 : eUtilityMode = ROUGHNESS;
1658 : }
1659 : else
1660 : {
1661 0 : fprintf(stderr, "Missing valid sub-utility mention\n");
1662 0 : Usage();
1663 0 : exit( 1 );
1664 : }
1665 :
1666 : /* -------------------------------------------------------------------- */
1667 : /* Parse arguments. */
1668 : /* -------------------------------------------------------------------- */
1669 41 : for(int i = 2; i < argc; i++ )
1670 : {
1671 48 : if( eUtilityMode == HILL_SHADE && i + 1 < argc &&
1672 14 : (EQUAL(argv[i], "--z") || EQUAL(argv[i], "-z")))
1673 : {
1674 2 : z = atof(argv[++i]);
1675 : }
1676 30 : else if ( eUtilityMode == SLOPE && EQUAL(argv[i], "-p"))
1677 : {
1678 0 : slopeFormat = 0;
1679 : }
1680 30 : else if ( eUtilityMode == ASPECT && EQUAL(argv[i], "-trigonometric"))
1681 : {
1682 0 : bAngleAsAzimuth = FALSE;
1683 : }
1684 30 : else if ( eUtilityMode == ASPECT && EQUAL(argv[i], "-zero_for_flat"))
1685 : {
1686 0 : bZeroForFlat = TRUE;
1687 : }
1688 30 : else if ( eUtilityMode == COLOR_RELIEF && EQUAL(argv[i], "-exact_color_entry"))
1689 : {
1690 0 : eColorSelectionMode = COLOR_SELECTION_EXACT_ENTRY;
1691 : }
1692 31 : else if ( eUtilityMode == COLOR_RELIEF && EQUAL(argv[i], "-nearest_color_entry"))
1693 : {
1694 1 : eColorSelectionMode = COLOR_SELECTION_NEAREST_ENTRY;
1695 : }
1696 106 : else if( i + 1 < argc &&
1697 20 : (EQUAL(argv[i], "--s") ||
1698 20 : EQUAL(argv[i], "-s") ||
1699 17 : EQUAL(argv[i], "--scale") ||
1700 17 : EQUAL(argv[i], "-scale"))
1701 : )
1702 : {
1703 3 : scale = atof(argv[++i]);
1704 : }
1705 38 : else if( eUtilityMode == HILL_SHADE && i + 1 < argc &&
1706 3 : (EQUAL(argv[i], "--az") ||
1707 3 : EQUAL(argv[i], "-az") ||
1708 3 : EQUAL(argv[i], "--azimuth") ||
1709 3 : EQUAL(argv[i], "-azimuth"))
1710 : )
1711 : {
1712 0 : az = atof(argv[++i]);
1713 : }
1714 38 : else if( eUtilityMode == HILL_SHADE && i + 1 < argc &&
1715 3 : (EQUAL(argv[i], "--alt") ||
1716 3 : EQUAL(argv[i], "-alt") ||
1717 3 : EQUAL(argv[i], "--alt") ||
1718 3 : EQUAL(argv[i], "-alt"))
1719 : )
1720 : {
1721 0 : alt = atof(argv[++i]);
1722 : }
1723 43 : else if( eUtilityMode == COLOR_RELIEF &&
1724 17 : EQUAL(argv[i], "-alpha"))
1725 : {
1726 0 : bAddAlpha = TRUE;
1727 : }
1728 60 : else if( i + 1 < argc &&
1729 17 : (EQUAL(argv[i], "--b") ||
1730 17 : EQUAL(argv[i], "-b"))
1731 : )
1732 : {
1733 0 : nBand = atoi(argv[++i]);
1734 : }
1735 26 : else if ( EQUAL(argv[i], "-q") || EQUAL(argv[i], "-quiet") )
1736 : {
1737 0 : pfnProgress = GDALDummyProgress;
1738 : }
1739 26 : else if( EQUAL(argv[i],"-co") && i < argc-1 )
1740 : {
1741 0 : papszCreateOptions = CSLAddString( papszCreateOptions, argv[++i] );
1742 : }
1743 29 : else if( EQUAL(argv[i],"-of") && i < argc-1 )
1744 : {
1745 3 : pszFormat = argv[++i];
1746 : }
1747 23 : else if( argv[i][0] == '-' )
1748 : {
1749 : fprintf( stderr, "Option %s incomplete, or not recognised.\n\n",
1750 0 : argv[i] );
1751 0 : Usage();
1752 0 : GDALDestroyDriverManager();
1753 0 : exit( 2 );
1754 : }
1755 23 : else if( pszSrcFilename == NULL )
1756 : {
1757 9 : pszSrcFilename = argv[i];
1758 : }
1759 19 : else if( eUtilityMode == COLOR_RELIEF && pszColorFilename == NULL )
1760 : {
1761 5 : pszColorFilename = argv[i];
1762 : }
1763 9 : else if( pszDstFilename == NULL )
1764 : {
1765 9 : pszDstFilename = argv[i];
1766 : }
1767 : else
1768 0 : Usage();
1769 : }
1770 :
1771 9 : if( pszSrcFilename == NULL )
1772 : {
1773 0 : fprintf( stderr, "Missing source.\n\n" );
1774 0 : Usage();
1775 : }
1776 9 : if ( eUtilityMode == COLOR_RELIEF && pszColorFilename == NULL )
1777 : {
1778 0 : fprintf( stderr, "Missing color file.\n\n" );
1779 0 : Usage();
1780 : }
1781 9 : if( pszDstFilename == NULL )
1782 : {
1783 0 : fprintf( stderr, "Missing destination.\n\n" );
1784 0 : Usage();
1785 : }
1786 :
1787 9 : GDALAllRegister();
1788 :
1789 : // Open Dataset and get raster band
1790 9 : hSrcDataset = GDALOpen( pszSrcFilename, GA_ReadOnly );
1791 :
1792 9 : if( hSrcDataset == NULL )
1793 : {
1794 : fprintf( stderr,
1795 : "GDALOpen failed - %d\n%s\n",
1796 0 : CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
1797 0 : GDALDestroyDriverManager();
1798 0 : exit( 1 );
1799 : }
1800 :
1801 9 : nXSize = GDALGetRasterXSize(hSrcDataset);
1802 9 : nYSize = GDALGetRasterYSize(hSrcDataset);
1803 :
1804 9 : if( nBand <= 0 || nBand > GDALGetRasterCount(hSrcDataset) )
1805 : {
1806 : fprintf( stderr,
1807 0 : "Unable to fetch band #%d\n", nBand );
1808 0 : GDALDestroyDriverManager();
1809 0 : exit( 1 );
1810 : }
1811 9 : hSrcBand = GDALGetRasterBand( hSrcDataset, nBand );
1812 :
1813 9 : GDALGetGeoTransform(hSrcDataset, adfGeoTransform);
1814 :
1815 9 : hDriver = GDALGetDriverByName(pszFormat);
1816 9 : if( hDriver == NULL
1817 : || (GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL &&
1818 : GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, NULL ) == NULL))
1819 : {
1820 : int iDr;
1821 :
1822 : fprintf( stderr, "Output driver `%s' not recognised or does not support\n",
1823 0 : pszFormat );
1824 : fprintf( stderr, "direct output file creation. The following format drivers are configured\n"
1825 0 : "and support direct output:\n" );
1826 :
1827 0 : for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
1828 : {
1829 0 : GDALDriverH hDriver = GDALGetDriver(iDr);
1830 :
1831 0 : if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL &&
1832 : GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, NULL ) == NULL )
1833 : {
1834 : printf( " %s: %s\n",
1835 : GDALGetDriverShortName( hDriver ),
1836 0 : GDALGetDriverLongName( hDriver ) );
1837 : }
1838 : }
1839 0 : GDALDestroyDriverManager();
1840 0 : exit( 1 );
1841 : }
1842 :
1843 9 : double dfDstNoDataValue = 0;
1844 9 : int bDstHasNoData = FALSE;
1845 9 : void* pData = NULL;
1846 9 : GDALGeneric3x3ProcessingAlg pfnAlg = NULL;
1847 :
1848 9 : if (eUtilityMode == HILL_SHADE)
1849 : {
1850 2 : dfDstNoDataValue = 0;
1851 2 : bDstHasNoData = TRUE;
1852 : pData = GDALCreateHillshadeData (adfGeoTransform,
1853 : z,
1854 : scale,
1855 : alt,
1856 2 : az);
1857 2 : pfnAlg = GDALHillshadeAlg;
1858 : }
1859 7 : else if (eUtilityMode == SLOPE)
1860 : {
1861 1 : dfDstNoDataValue = -9999;
1862 1 : bDstHasNoData = TRUE;
1863 :
1864 1 : pData = GDALCreateSlopeData(adfGeoTransform, scale, slopeFormat);
1865 1 : pfnAlg = GDALSlopeAlg;
1866 : }
1867 :
1868 6 : else if (eUtilityMode == ASPECT)
1869 : {
1870 1 : if (!bZeroForFlat)
1871 : {
1872 1 : dfDstNoDataValue = -9999;
1873 1 : bDstHasNoData = TRUE;
1874 : }
1875 :
1876 1 : pData = GDALCreateAspectData(bAngleAsAzimuth);
1877 1 : pfnAlg = GDALAspectAlg;
1878 : }
1879 5 : else if (eUtilityMode == TRI)
1880 : {
1881 0 : dfDstNoDataValue = -9999;
1882 0 : bDstHasNoData = TRUE;
1883 0 : pfnAlg = GDALTRIAlg;
1884 : }
1885 5 : else if (eUtilityMode == TPI)
1886 : {
1887 0 : dfDstNoDataValue = -9999;
1888 0 : bDstHasNoData = TRUE;
1889 0 : pfnAlg = GDALTPIAlg;
1890 : }
1891 5 : else if (eUtilityMode == ROUGHNESS)
1892 : {
1893 0 : dfDstNoDataValue = -9999;
1894 0 : bDstHasNoData = TRUE;
1895 0 : pfnAlg = GDALRoughnessAlg;
1896 : }
1897 :
1898 : GDALDataType eDstDataType = (eUtilityMode == HILL_SHADE ||
1899 : eUtilityMode == COLOR_RELIEF) ? GDT_Byte :
1900 9 : GDT_Float32;
1901 :
1902 9 : if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL &&
1903 : GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, NULL ) != NULL)
1904 : {
1905 : GDALDatasetH hIntermediateDataset;
1906 :
1907 3 : if (eUtilityMode == COLOR_RELIEF)
1908 : hIntermediateDataset = (GDALDatasetH)
1909 : new GDALColorReliefDataset (hSrcDataset,
1910 : hSrcBand,
1911 : pszColorFilename,
1912 : eColorSelectionMode,
1913 2 : bAddAlpha);
1914 : else
1915 : hIntermediateDataset = (GDALDatasetH)
1916 : new GDALGeneric3x3Dataset(hSrcDataset, hSrcBand,
1917 : eDstDataType,
1918 : bDstHasNoData,
1919 : dfDstNoDataValue,
1920 : pfnAlg,
1921 1 : pData);
1922 :
1923 : GDALDatasetH hOutDS = GDALCreateCopy(
1924 : hDriver, pszDstFilename, hIntermediateDataset,
1925 : TRUE, papszCreateOptions,
1926 3 : pfnProgress, NULL );
1927 :
1928 3 : if( hOutDS != NULL )
1929 3 : GDALClose( hOutDS );
1930 3 : GDALClose(hIntermediateDataset);
1931 3 : GDALClose(hSrcDataset);
1932 :
1933 3 : CPLFree(pData);
1934 :
1935 3 : GDALDestroyDriverManager();
1936 3 : CSLDestroy( argv );
1937 3 : CSLDestroy( papszCreateOptions );
1938 :
1939 3 : return 0;
1940 : }
1941 :
1942 : int nDstBands;
1943 6 : if (eUtilityMode == COLOR_RELIEF)
1944 3 : nDstBands = (bAddAlpha) ? 4 : 3;
1945 : else
1946 3 : nDstBands = 1;
1947 :
1948 : hDstDataset = GDALCreate( hDriver,
1949 : pszDstFilename,
1950 : nXSize,
1951 : nYSize,
1952 : nDstBands,
1953 : eDstDataType,
1954 6 : papszCreateOptions);
1955 :
1956 6 : if( hDstDataset == NULL )
1957 : {
1958 : fprintf( stderr,
1959 : "Unable to create dataset %s %d\n%s\n", pszDstFilename,
1960 0 : CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
1961 0 : GDALDestroyDriverManager();
1962 0 : exit( 1 );
1963 : }
1964 :
1965 6 : hDstBand = GDALGetRasterBand( hDstDataset, 1 );
1966 :
1967 6 : GDALSetGeoTransform(hDstDataset, adfGeoTransform);
1968 6 : GDALSetProjection(hDstDataset, GDALGetProjectionRef(hSrcDataset));
1969 :
1970 6 : if (eUtilityMode == COLOR_RELIEF)
1971 : {
1972 : GDALColorRelief (hSrcBand,
1973 : GDALGetRasterBand(hDstDataset, 1),
1974 : GDALGetRasterBand(hDstDataset, 2),
1975 : GDALGetRasterBand(hDstDataset, 3),
1976 : (bAddAlpha) ? GDALGetRasterBand(hDstDataset, 4) : NULL,
1977 : pszColorFilename,
1978 : eColorSelectionMode,
1979 3 : pfnProgress, NULL);
1980 : }
1981 : else
1982 : {
1983 3 : if (bDstHasNoData)
1984 3 : GDALSetRasterNoDataValue(hDstBand, dfDstNoDataValue);
1985 :
1986 : GDALGeneric3x3Processing(hSrcBand, hDstBand,
1987 : pfnAlg, pData,
1988 3 : pfnProgress, NULL);
1989 :
1990 : }
1991 :
1992 6 : GDALClose(hSrcDataset);
1993 6 : GDALClose(hDstDataset);
1994 6 : CPLFree(pData);
1995 :
1996 6 : GDALDestroyDriverManager();
1997 6 : CSLDestroy( argv );
1998 6 : CSLDestroy( papszCreateOptions );
1999 :
2000 6 : return 0;
2001 : }
2002 :
|
