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