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