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 191060 : static float ComputeVal(int bSrcHasNoData, float fSrcNoDataValue,
161 : float* afWin, float fDstNoDataValue,
162 : GDALGeneric3x3ProcessingAlg pfnAlg,
163 : void* pData,
164 : int bComputeAtEdges)
165 : {
166 191060 : if (bSrcHasNoData && ARE_REAL_EQUAL(afWin[4], fSrcNoDataValue))
167 : {
168 0 : return fDstNoDataValue;
169 : }
170 191060 : else if (bSrcHasNoData)
171 : {
172 : int k;
173 1718520 : for(k=0;k<9;k++)
174 : {
175 1546668 : 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 191060 : return pfnAlg(afWin, fDstNoDataValue, pData);
186 : }
187 :
188 : /************************************************************************/
189 : /* GDALGeneric3x3Processing() */
190 : /************************************************************************/
191 :
192 10 : 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 10 : float fSrcNoDataValue = 0.0, fDstNoDataValue = 0.0;
207 :
208 10 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
209 10 : int nYSize = GDALGetRasterBandYSize(hSrcBand);
210 :
211 10 : if (pfnProgress == NULL)
212 0 : pfnProgress = GDALDummyProgress;
213 :
214 : /* -------------------------------------------------------------------- */
215 : /* Initialize progress counter. */
216 : /* -------------------------------------------------------------------- */
217 10 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
218 : {
219 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
220 0 : return CE_Failure;
221 : }
222 :
223 10 : pafOutputBuf = (float *) CPLMalloc(sizeof(float)*nXSize);
224 10 : pafThreeLineWin = (float *) CPLMalloc(3*sizeof(float)*(nXSize+1));
225 :
226 10 : fSrcNoDataValue = (float) GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
227 10 : fDstNoDataValue = (float) GDALGetRasterNoDataValue(hDstBand, &bDstHasNoData);
228 10 : 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 30 : 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 20 : 0, 0);
249 : }
250 :
251 12 : if (bComputeAtEdges && nXSize >= 2 && nYSize >= 2)
252 : {
253 244 : for (j = 0; j < nXSize; j++)
254 : {
255 : float afWin[9];
256 242 : int jmin = (j == 0) ? j : j - 1;
257 242 : int jmax = (j == nXSize - 1) ? j : j + 1;
258 :
259 242 : afWin[0] = INTERPOL(pafThreeLineWin[jmin], pafThreeLineWin[nXSize + jmin]);
260 242 : afWin[1] = INTERPOL(pafThreeLineWin[j], pafThreeLineWin[nXSize + j]);
261 242 : afWin[2] = INTERPOL(pafThreeLineWin[jmax], pafThreeLineWin[nXSize + jmax]);
262 242 : afWin[3] = pafThreeLineWin[jmin];
263 242 : afWin[4] = pafThreeLineWin[j];
264 242 : afWin[5] = pafThreeLineWin[jmax];
265 242 : afWin[6] = pafThreeLineWin[nXSize + jmin];
266 242 : afWin[7] = pafThreeLineWin[nXSize + j];
267 242 : afWin[8] = pafThreeLineWin[nXSize + jmax];
268 :
269 242 : pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
270 : afWin, fDstNoDataValue,
271 242 : pfnAlg, pData, bComputeAtEdges);
272 : }
273 : GDALRasterIO(hDstBand, GF_Write,
274 : 0, 0, nXSize, 1,
275 2 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
276 : }
277 : else
278 : {
279 : // Exclude the edges
280 934 : for (j = 0; j < nXSize; j++)
281 : {
282 926 : pafOutputBuf[j] = fDstNoDataValue;
283 : }
284 : GDALRasterIO(hDstBand, GF_Write,
285 : 0, 0, nXSize, 1,
286 8 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
287 :
288 8 : if (nYSize > 1)
289 : {
290 : GDALRasterIO(hDstBand, GF_Write,
291 : 0, nYSize - 1, nXSize, 1,
292 8 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
293 : }
294 : }
295 :
296 10 : int nLine1Off = 0*nXSize;
297 10 : int nLine2Off = 1*nXSize;
298 10 : int nLine3Off = 2*nXSize;
299 :
300 1158 : 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 1148 : 0, 0);
311 1148 : if (eErr != CE_None)
312 0 : goto end;
313 :
314 1386 : if (bComputeAtEdges && nXSize >= 2)
315 : {
316 : float afWin[9];
317 :
318 238 : j = 0;
319 238 : afWin[0] = INTERPOL(pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j+1]);
320 238 : afWin[1] = pafThreeLineWin[nLine1Off + j];
321 238 : afWin[2] = pafThreeLineWin[nLine1Off + j+1];
322 238 : afWin[3] = INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j+1]);
323 238 : afWin[4] = pafThreeLineWin[nLine2Off + j];
324 238 : afWin[5] = pafThreeLineWin[nLine2Off + j+1];
325 238 : afWin[6] = INTERPOL(pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j+1]);
326 238 : afWin[7] = pafThreeLineWin[nLine3Off + j];
327 238 : afWin[8] = pafThreeLineWin[nLine3Off + j+1];
328 :
329 238 : pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
330 : afWin, fDstNoDataValue,
331 238 : pfnAlg, pData, bComputeAtEdges);
332 238 : j = nXSize - 1;
333 :
334 238 : afWin[0] = pafThreeLineWin[nLine1Off + j-1];
335 238 : afWin[1] = pafThreeLineWin[nLine1Off + j];
336 238 : afWin[2] = INTERPOL(pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j-1]);
337 238 : afWin[3] = pafThreeLineWin[nLine2Off + j-1];
338 238 : afWin[4] = pafThreeLineWin[nLine2Off + j];
339 238 : afWin[5] = INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j-1]);
340 238 : afWin[6] = pafThreeLineWin[nLine3Off + j-1];
341 238 : afWin[7] = pafThreeLineWin[nLine3Off + j];
342 238 : afWin[8] = INTERPOL(pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j-1]);
343 :
344 238 : pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
345 : afWin, fDstNoDataValue,
346 238 : pfnAlg, pData, bComputeAtEdges);
347 : }
348 : else
349 : {
350 : // Exclude the edges
351 910 : pafOutputBuf[0] = fDstNoDataValue;
352 910 : if (nXSize > 1)
353 910 : pafOutputBuf[nXSize - 1] = fDstNoDataValue;
354 : }
355 :
356 133644 : for (j = 1; j < nXSize - 1; j++)
357 : {
358 : float afWin[9];
359 132496 : afWin[0] = pafThreeLineWin[nLine1Off + j-1];
360 132496 : afWin[1] = pafThreeLineWin[nLine1Off + j];
361 132496 : afWin[2] = pafThreeLineWin[nLine1Off + j+1];
362 132496 : afWin[3] = pafThreeLineWin[nLine2Off + j-1];
363 132496 : afWin[4] = pafThreeLineWin[nLine2Off + j];
364 132496 : afWin[5] = pafThreeLineWin[nLine2Off + j+1];
365 132496 : afWin[6] = pafThreeLineWin[nLine3Off + j-1];
366 132496 : afWin[7] = pafThreeLineWin[nLine3Off + j];
367 132496 : afWin[8] = pafThreeLineWin[nLine3Off + j+1];
368 :
369 132496 : pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
370 : afWin, fDstNoDataValue,
371 132496 : pfnAlg, pData, bComputeAtEdges);
372 : }
373 :
374 : /* -----------------------------------------
375 : * Write Line to Raster
376 : */
377 : eErr = GDALRasterIO(hDstBand, GF_Write, 0, i, nXSize, 1,
378 1148 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
379 1148 : if (eErr != CE_None)
380 0 : goto end;
381 :
382 1148 : 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 1148 : int nTemp = nLine1Off;
390 1148 : nLine1Off = nLine2Off;
391 1148 : nLine2Off = nLine3Off;
392 1148 : nLine3Off = nTemp;
393 : }
394 :
395 10 : if (bComputeAtEdges && nXSize >= 2 && nYSize >= 2)
396 : {
397 244 : for (j = 0; j < nXSize; j++)
398 : {
399 : float afWin[9];
400 242 : int jmin = (j == 0) ? j : j - 1;
401 242 : int jmax = (j == nXSize - 1) ? j : j + 1;
402 :
403 242 : afWin[0] = pafThreeLineWin[nLine1Off + jmin];
404 242 : afWin[1] = pafThreeLineWin[nLine1Off + j];
405 242 : afWin[2] = pafThreeLineWin[nLine1Off + jmax];
406 242 : afWin[3] = pafThreeLineWin[nLine2Off + jmin];
407 242 : afWin[4] = pafThreeLineWin[nLine2Off + j];
408 242 : afWin[5] = pafThreeLineWin[nLine2Off + jmax];
409 242 : afWin[6] = INTERPOL(pafThreeLineWin[nLine2Off + jmin], pafThreeLineWin[nLine1Off + jmin]);
410 242 : afWin[7] = INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine1Off + j]);
411 242 : afWin[8] = INTERPOL(pafThreeLineWin[nLine2Off + jmax], pafThreeLineWin[nLine1Off + jmax]);
412 :
413 242 : pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
414 : afWin, fDstNoDataValue,
415 242 : pfnAlg, pData, bComputeAtEdges);
416 : }
417 : GDALRasterIO(hDstBand, GF_Write,
418 : 0, i, nXSize, 1,
419 2 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
420 : }
421 :
422 10 : pfnProgress( 1.0, NULL, pProgressData );
423 10 : eErr = CE_None;
424 :
425 : end:
426 10 : CPLFree(pafOutputBuf);
427 10 : CPLFree(pafThreeLineWin);
428 :
429 10 : 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 134416 : float GDALHillshadeAlg (float* afWin, float fDstNoDataValue, void* pData)
466 : {
467 134416 : GDALHillshadeAlgData* psData = (GDALHillshadeAlgData*)pData;
468 : double x, y, aspect, xx_plus_yy, cang;
469 :
470 : // First Slope ...
471 403248 : x = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
472 403248 : (afWin[2] + afWin[5] + afWin[5] + afWin[8])) / psData->ewres;
473 :
474 403248 : y = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
475 403248 : (afWin[0] + afWin[1] + afWin[1] + afWin[2])) / psData->nsres;
476 :
477 134416 : xx_plus_yy = x * x + y * y;
478 :
479 : // ... then aspect...
480 134416 : 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 134416 : sqrt(1 + psData->square_z_scale_factor * xx_plus_yy);
487 :
488 134416 : if (cang <= 0.0)
489 832 : cang = 1.0;
490 : else
491 133584 : cang = 1.0 + (254.0 * cang);
492 :
493 134416 : 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 10 : void* GDALCreateHillshadeData(double* adfGeoTransform,
526 : double z,
527 : double scale,
528 : double alt,
529 : double az,
530 : int bZevenbergenThorne)
531 : {
532 : GDALHillshadeAlgData* pData =
533 10 : (GDALHillshadeAlgData*)CPLMalloc(sizeof(GDALHillshadeAlgData));
534 :
535 10 : const double degreesToRadians = M_PI / 180.0;
536 10 : pData->nsres = adfGeoTransform[5];
537 10 : pData->ewres = adfGeoTransform[1];
538 10 : pData->sin_altRadians = sin(alt * degreesToRadians);
539 10 : pData->azRadians = az * degreesToRadians;
540 10 : double z_scale_factor = z / (((bZevenbergenThorne) ? 2 : 8) * scale);
541 : pData->cos_altRadians_mul_z_scale_factor =
542 10 : cos(alt * degreesToRadians) * z_scale_factor;
543 10 : pData->square_z_scale_factor = z_scale_factor * z_scale_factor;
544 10 : 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 28322 : float GDALSlopeHornAlg (float* afWin, float fDstNoDataValue, void* pData)
560 : {
561 28322 : const double radiansToDegrees = 180.0 / M_PI;
562 28322 : GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
563 : double dx, dy, key;
564 :
565 84966 : dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
566 84966 : (afWin[2] + afWin[5] + afWin[5] + afWin[8]))/psData->ewres;
567 :
568 84966 : dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
569 84966 : (afWin[0] + afWin[1] + afWin[1] + afWin[2]))/psData->nsres;
570 :
571 28322 : key = (dx * dx + dy * dy);
572 :
573 28322 : if (psData->slopeFormat == 1)
574 28322 : 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 2 : void* GDALCreateSlopeData(double* adfGeoTransform,
598 : double scale,
599 : int slopeFormat)
600 : {
601 : GDALSlopeAlgData* pData =
602 2 : (GDALSlopeAlgData*)CPLMalloc(sizeof(GDALSlopeAlgData));
603 :
604 2 : pData->nsres = adfGeoTransform[5];
605 2 : pData->ewres = adfGeoTransform[1];
606 2 : pData->scale = scale;
607 2 : pData->slopeFormat = slopeFormat;
608 2 : return pData;
609 : }
610 :
611 : /************************************************************************/
612 : /* GDALAspect() */
613 : /************************************************************************/
614 :
615 : typedef struct
616 : {
617 : int bAngleAsAzimuth;
618 : } GDALAspectAlgData;
619 :
620 28322 : float GDALAspectAlg (float* afWin, float fDstNoDataValue, void* pData)
621 : {
622 28322 : const double degreesToRadians = M_PI / 180.0;
623 28322 : GDALAspectAlgData* psData = (GDALAspectAlgData*)pData;
624 : double dx, dy;
625 : float aspect;
626 :
627 84966 : dx = ((afWin[2] + afWin[5] + afWin[5] + afWin[8]) -
628 84966 : (afWin[0] + afWin[3] + afWin[3] + afWin[6]));
629 :
630 84966 : dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
631 84966 : (afWin[0] + afWin[1] + afWin[1] + afWin[2]));
632 :
633 28322 : aspect = (float) (atan2(dy,-dx) / degreesToRadians);
634 :
635 36688 : if (dx == 0 && dy == 0)
636 : {
637 : /* Flat area */
638 8366 : aspect = fDstNoDataValue;
639 : }
640 19956 : else if ( psData->bAngleAsAzimuth )
641 : {
642 19956 : if (aspect > 90.0)
643 2486 : aspect = 450.0f - aspect;
644 : else
645 17470 : aspect = 90.0f - aspect;
646 : }
647 : else
648 : {
649 0 : if (aspect < 0)
650 0 : aspect += 360.0;
651 : }
652 :
653 28322 : if (aspect == 360.0)
654 0 : aspect = 0.0;
655 :
656 28322 : 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 2 : void* GDALCreateAspectData(int bAngleAsAzimuth)
696 : {
697 : GDALAspectAlgData* pData =
698 2 : (GDALAspectAlgData*)CPLMalloc(sizeof(GDALAspectAlgData));
699 :
700 2 : pData->bAngleAsAzimuth = bAngleAsAzimuth;
701 2 : 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 176 : static int GDALColorReliefSortColors(const void* pA, const void* pB)
718 : {
719 176 : ColorAssociation* pC1 = (ColorAssociation*)pA;
720 176 : ColorAssociation* pC2 = (ColorAssociation*)pB;
721 : return (pC1->dfVal < pC2->dfVal) ? -1 :
722 176 : (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 292820 : 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 292820 : int lower = 0;
743 292820 : 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 641740 : while(TRUE)
749 : {
750 934560 : mid = (lower + upper) / 2;
751 934560 : if (upper - lower <= 1)
752 : {
753 292820 : if (dfVal < pasColorAssociation[lower].dfVal)
754 0 : i = lower;
755 292820 : else if (dfVal < pasColorAssociation[upper].dfVal)
756 199300 : i = upper;
757 : else
758 93520 : i = upper + 1;
759 : break;
760 : }
761 641740 : else if (pasColorAssociation[mid].dfVal >= dfVal)
762 : {
763 385260 : upper = mid;
764 : }
765 : else
766 : {
767 256480 : lower = mid;
768 : }
769 : }
770 :
771 292820 : 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 292820 : 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 292820 : 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 322102 : if (eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY &&
824 29282 : pasColorAssociation[i-1].dfVal != dfVal)
825 : {
826 : int index;
827 39860 : if (dfVal - pasColorAssociation[i-1].dfVal <
828 19930 : pasColorAssociation[i].dfVal - dfVal)
829 13432 : index = i -1;
830 : else
831 6498 : index = i;
832 :
833 19930 : *pnR = pasColorAssociation[index].nR;
834 19930 : *pnG = pasColorAssociation[index].nG;
835 19930 : *pnB = pasColorAssociation[index].nB;
836 19930 : *pnA = pasColorAssociation[index].nA;
837 19930 : return TRUE;
838 : }
839 :
840 272890 : double dfRatio = (dfVal - pasColorAssociation[i-1].dfVal) /
841 272890 : (pasColorAssociation[i].dfVal - pasColorAssociation[i-1].dfVal);
842 272890 : *pnR = (int)(0.45 + pasColorAssociation[i-1].nR + dfRatio *
843 272890 : (pasColorAssociation[i].nR - pasColorAssociation[i-1].nR));
844 272890 : if (*pnR < 0) *pnR = 0;
845 272890 : else if (*pnR > 255) *pnR = 255;
846 272890 : *pnG = (int)(0.45 + pasColorAssociation[i-1].nG + dfRatio *
847 272890 : (pasColorAssociation[i].nG - pasColorAssociation[i-1].nG));
848 272890 : if (*pnG < 0) *pnG = 0;
849 272890 : else if (*pnG > 255) *pnG = 255;
850 272890 : *pnB = (int)(0.45 + pasColorAssociation[i-1].nB + dfRatio *
851 272890 : (pasColorAssociation[i].nB - pasColorAssociation[i-1].nB));
852 272890 : if (*pnB < 0) *pnB = 0;
853 272890 : else if (*pnB > 255) *pnB = 255;
854 272890 : *pnA = (int)(0.45 + pasColorAssociation[i-1].nA + dfRatio *
855 272890 : (pasColorAssociation[i].nA - pasColorAssociation[i-1].nA));
856 272890 : if (*pnA < 0) *pnA = 0;
857 272890 : else if (*pnA > 255) *pnA = 255;
858 :
859 272890 : 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 16 : ColorAssociation* GDALColorReliefParseColorFile(GDALRasterBandH hSrcBand,
927 : const char* pszColorFilename,
928 : int* pnColors)
929 : {
930 16 : VSILFILE* fpColorFile = VSIFOpenL(pszColorFilename, "rt");
931 16 : 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 16 : ColorAssociation* pasColorAssociation = NULL;
939 16 : int nColorAssociation = 0;
940 :
941 16 : int bSrcHasNoData = FALSE;
942 16 : double dfSrcNoDataValue = GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
943 :
944 : const char* pszLine;
945 16 : int bIsGMT_CPT = FALSE;
946 148 : while ((pszLine = CPLReadLineL(fpColorFile)) != NULL)
947 : {
948 116 : if (pszLine[0] == '#' && strstr(pszLine, "COLOR_MODEL"))
949 : {
950 2 : 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 2 : bIsGMT_CPT = TRUE;
958 : }
959 :
960 : char** papszFields = CSLTokenizeStringComplex(pszLine, " ,\t:",
961 116 : FALSE, FALSE );
962 : /* Skip comment lines */
963 116 : int nTokens = CSLCount(papszFields);
964 226 : if (nTokens >= 1 && (papszFields[0][0] == '#' ||
965 110 : papszFields[0][0] == '/'))
966 : {
967 6 : CSLDestroy(papszFields);
968 6 : continue;
969 : }
970 :
971 116 : if (bIsGMT_CPT && nTokens == 8)
972 : {
973 : pasColorAssociation =
974 : (ColorAssociation*)CPLRealloc(pasColorAssociation,
975 6 : (nColorAssociation + 2) * sizeof(ColorAssociation));
976 :
977 6 : pasColorAssociation[nColorAssociation].dfVal = atof(papszFields[0]);
978 6 : pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]);
979 6 : pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]);
980 6 : pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]);
981 6 : pasColorAssociation[nColorAssociation].nA = 255;
982 6 : nColorAssociation++;
983 :
984 6 : pasColorAssociation[nColorAssociation].dfVal = atof(papszFields[4]);
985 6 : pasColorAssociation[nColorAssociation].nR = atoi(papszFields[5]);
986 6 : pasColorAssociation[nColorAssociation].nG = atoi(papszFields[6]);
987 6 : pasColorAssociation[nColorAssociation].nB = atoi(papszFields[7]);
988 6 : pasColorAssociation[nColorAssociation].nA = 255;
989 6 : nColorAssociation++;
990 : }
991 110 : 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 6 : if (EQUAL(papszFields[0], "N") && bSrcHasNoData)
996 : {
997 : pasColorAssociation =
998 : (ColorAssociation*)CPLRealloc(pasColorAssociation,
999 2 : (nColorAssociation + 1) * sizeof(ColorAssociation));
1000 :
1001 2 : pasColorAssociation[nColorAssociation].dfVal = dfSrcNoDataValue;
1002 2 : pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]);
1003 2 : pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]);
1004 2 : pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]);
1005 2 : pasColorAssociation[nColorAssociation].nA = 255;
1006 2 : nColorAssociation++;
1007 : }
1008 : }
1009 98 : else if (!bIsGMT_CPT && nTokens >= 2)
1010 : {
1011 : pasColorAssociation =
1012 : (ColorAssociation*)CPLRealloc(pasColorAssociation,
1013 98 : (nColorAssociation + 1) * sizeof(ColorAssociation));
1014 98 : if (EQUAL(papszFields[0], "nv") && bSrcHasNoData)
1015 0 : pasColorAssociation[nColorAssociation].dfVal = dfSrcNoDataValue;
1016 98 : 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 98 : pasColorAssociation[nColorAssociation].dfVal = atof(papszFields[0]);
1034 :
1035 98 : if (nTokens >= 4)
1036 : {
1037 98 : pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]);
1038 98 : pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]);
1039 98 : pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]);
1040 98 : pasColorAssociation[nColorAssociation].nA =
1041 98 : (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 98 : nColorAssociation ++;
1064 : }
1065 110 : CSLDestroy(papszFields);
1066 : }
1067 16 : VSIFCloseL(fpColorFile);
1068 :
1069 16 : 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 16 : sizeof(ColorAssociation), GDALColorReliefSortColors);
1079 :
1080 16 : *pnColors = nColorAssociation;
1081 16 : return pasColorAssociation;
1082 : }
1083 :
1084 : static
1085 12 : GByte* GDALColorReliefPrecompute(GDALRasterBandH hSrcBand,
1086 : ColorAssociation* pasColorAssociation,
1087 : int nColorAssociation,
1088 : ColorSelectionMode eColorSelectionMode,
1089 : int* pnIndexOffset)
1090 : {
1091 12 : GDALDataType eDT = GDALGetRasterDataType(hSrcBand);
1092 12 : GByte* pabyPrecomputed = NULL;
1093 12 : int nIndexOffset = (eDT == GDT_Int16) ? 32768 : 0;
1094 12 : *pnIndexOffset = nIndexOffset;
1095 12 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
1096 12 : int nYSize = GDALGetRasterBandXSize(hSrcBand);
1097 12 : 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 12 : 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 12 : {
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 4 : GDALColorReliefDataset::GDALColorReliefDataset(
1178 : GDALDatasetH hSrcDS,
1179 : GDALRasterBandH hSrcBand,
1180 : const char* pszColorFilename,
1181 : ColorSelectionMode eColorSelectionMode,
1182 4 : int bAlpha)
1183 : {
1184 4 : this->hSrcDS = hSrcDS;
1185 4 : this->hSrcBand = hSrcBand;
1186 4 : nColorAssociation = 0;
1187 : pasColorAssociation =
1188 : GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
1189 4 : &nColorAssociation);
1190 4 : this->eColorSelectionMode = eColorSelectionMode;
1191 :
1192 4 : nRasterXSize = GDALGetRasterXSize(hSrcDS);
1193 4 : nRasterYSize = GDALGetRasterYSize(hSrcDS);
1194 :
1195 : int nBlockXSize, nBlockYSize;
1196 4 : GDALGetBlockSize( hSrcBand, &nBlockXSize, &nBlockYSize);
1197 :
1198 4 : nIndexOffset = 0;
1199 : pabyPrecomputed =
1200 : GDALColorReliefPrecompute(hSrcBand,
1201 : pasColorAssociation,
1202 : nColorAssociation,
1203 : eColorSelectionMode,
1204 4 : &nIndexOffset);
1205 :
1206 : int i;
1207 16 : for(i=0;i<((bAlpha) ? 4 : 3);i++)
1208 : {
1209 12 : SetBand(i + 1, new GDALColorReliefRasterBand(this, i+1));
1210 : }
1211 :
1212 4 : pafSourceBuf = NULL;
1213 4 : panSourceBuf = NULL;
1214 4 : if (pabyPrecomputed)
1215 0 : panSourceBuf = (int *) CPLMalloc(sizeof(int)*nBlockXSize*nBlockYSize);
1216 : else
1217 4 : pafSourceBuf = (float *) CPLMalloc(sizeof(float)*nBlockXSize*nBlockYSize);
1218 4 : nCurBlockXOff = -1;
1219 4 : nCurBlockYOff = -1;
1220 4 : }
1221 :
1222 4 : GDALColorReliefDataset::~GDALColorReliefDataset()
1223 : {
1224 4 : CPLFree(pasColorAssociation);
1225 4 : CPLFree(pabyPrecomputed);
1226 4 : CPLFree(panSourceBuf);
1227 4 : CPLFree(pafSourceBuf);
1228 4 : }
1229 :
1230 4 : CPLErr GDALColorReliefDataset::GetGeoTransform( double * padfGeoTransform )
1231 : {
1232 4 : return GDALGetGeoTransform(hSrcDS, padfGeoTransform);
1233 : }
1234 :
1235 4 : const char *GDALColorReliefDataset::GetProjectionRef()
1236 : {
1237 4 : return GDALGetProjectionRef(hSrcDS);
1238 : }
1239 :
1240 12 : GDALColorReliefRasterBand::GDALColorReliefRasterBand(
1241 12 : GDALColorReliefDataset * poDS, int nBand)
1242 : {
1243 12 : this->poDS = poDS;
1244 12 : this->nBand = nBand;
1245 12 : eDataType = GDT_Byte;
1246 12 : GDALGetBlockSize( poDS->hSrcBand, &nBlockXSize, &nBlockYSize);
1247 12 : }
1248 :
1249 774 : CPLErr GDALColorReliefRasterBand::IReadBlock( int nBlockXOff,
1250 : int nBlockYOff,
1251 : void *pImage )
1252 : {
1253 774 : GDALColorReliefDataset * poGDS = (GDALColorReliefDataset *) poDS;
1254 : int nReqXSize, nReqYSize;
1255 :
1256 774 : if ((nBlockXOff + 1) * nBlockXSize >= nRasterXSize)
1257 54 : nReqXSize = nRasterXSize - nBlockXOff * nBlockXSize;
1258 : else
1259 720 : nReqXSize = nBlockXSize;
1260 :
1261 774 : if ((nBlockYOff + 1) * nBlockYSize >= nRasterYSize)
1262 732 : nReqYSize = nRasterYSize - nBlockYOff * nBlockYSize;
1263 : else
1264 42 : nReqYSize = nBlockYSize;
1265 :
1266 774 : int nCount = nReqXSize * nReqYSize;
1267 774 : if ( poGDS->nCurBlockXOff != nBlockXOff ||
1268 : poGDS->nCurBlockYOff != nBlockYOff )
1269 : {
1270 742 : poGDS->nCurBlockXOff = nBlockXOff;
1271 742 : 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 742 : 0, 0);
1284 742 : if (eErr != CE_None)
1285 : {
1286 0 : memset(pImage, 0, nCount);
1287 0 : return eErr;
1288 : }
1289 : }
1290 :
1291 : int j;
1292 774 : 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 176466 : for ( j = 0; j < nCount; j++)
1304 : {
1305 : GDALColorReliefGetRGBA (poGDS->pasColorAssociation,
1306 : poGDS->nColorAssociation,
1307 175692 : poGDS->pafSourceBuf[j],
1308 : poGDS->eColorSelectionMode,
1309 : &anComponents[0],
1310 : &anComponents[1],
1311 : &anComponents[2],
1312 351384 : &anComponents[3]);
1313 175692 : ((GByte*)pImage)[j] = (GByte) anComponents[nBand-1];
1314 : }
1315 : }
1316 :
1317 774 : return CE_None;
1318 : }
1319 :
1320 24 : GDALColorInterp GDALColorReliefRasterBand::GetColorInterpretation()
1321 : {
1322 24 : return (GDALColorInterp)(GCI_RedBand + nBand - 1);
1323 : }
1324 :
1325 :
1326 8 : 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 8 : if (hSrcBand == NULL || hDstBand1 == NULL || hDstBand2 == NULL ||
1339 : hDstBand3 == NULL)
1340 0 : return CE_Failure;
1341 :
1342 8 : int nColorAssociation = 0;
1343 : ColorAssociation* pasColorAssociation =
1344 : GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
1345 8 : &nColorAssociation);
1346 8 : if (pasColorAssociation == NULL)
1347 0 : return CE_Failure;
1348 :
1349 8 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
1350 8 : int nYSize = GDALGetRasterBandYSize(hSrcBand);
1351 :
1352 8 : if (pfnProgress == NULL)
1353 0 : pfnProgress = GDALDummyProgress;
1354 :
1355 8 : 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 8 : int nIndexOffset = 0;
1362 : GByte* pabyPrecomputed =
1363 : GDALColorReliefPrecompute(hSrcBand,
1364 : pasColorAssociation,
1365 : nColorAssociation,
1366 : eColorSelectionMode,
1367 8 : &nIndexOffset);
1368 :
1369 : /* -------------------------------------------------------------------- */
1370 : /* Initialize progress counter. */
1371 : /* -------------------------------------------------------------------- */
1372 :
1373 8 : float* pafSourceBuf = NULL;
1374 8 : int* panSourceBuf = NULL;
1375 8 : if (pabyPrecomputed)
1376 0 : panSourceBuf = (int *) CPLMalloc(sizeof(int)*nXSize);
1377 : else
1378 8 : pafSourceBuf = (float *) CPLMalloc(sizeof(float)*nXSize);
1379 8 : GByte* pabyDestBuf1 = (GByte*) CPLMalloc( 4 * nXSize );
1380 8 : GByte* pabyDestBuf2 = pabyDestBuf1 + nXSize;
1381 8 : GByte* pabyDestBuf3 = pabyDestBuf2 + nXSize;
1382 8 : GByte* pabyDestBuf4 = pabyDestBuf3 + nXSize;
1383 : int i, j;
1384 :
1385 8 : 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 976 : 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 968 : 0, 0);
1403 968 : if (eErr != CE_None)
1404 0 : goto end;
1405 :
1406 968 : 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 118096 : for ( j = 0; j < nXSize; j++)
1420 : {
1421 : GDALColorReliefGetRGBA (pasColorAssociation,
1422 : nColorAssociation,
1423 117128 : pafSourceBuf[j],
1424 : eColorSelectionMode,
1425 : &nR,
1426 : &nG,
1427 : &nB,
1428 117128 : &nA);
1429 117128 : pabyDestBuf1[j] = (GByte) nR;
1430 117128 : pabyDestBuf2[j] = (GByte) nG;
1431 117128 : pabyDestBuf3[j] = (GByte) nB;
1432 117128 : 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 968 : 1, pabyDestBuf1, nXSize, 1, GDT_Byte, 0, 0);
1443 968 : if (eErr != CE_None)
1444 0 : goto end;
1445 :
1446 : eErr = GDALRasterIO(hDstBand2,
1447 : GF_Write,
1448 : 0, i, nXSize,
1449 968 : 1, pabyDestBuf2, nXSize, 1, GDT_Byte, 0, 0);
1450 968 : if (eErr != CE_None)
1451 0 : goto end;
1452 :
1453 : eErr = GDALRasterIO(hDstBand3,
1454 : GF_Write,
1455 : 0, i, nXSize,
1456 968 : 1, pabyDestBuf3, nXSize, 1, GDT_Byte, 0, 0);
1457 968 : if (eErr != CE_None)
1458 0 : goto end;
1459 :
1460 968 : 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 968 : 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 8 : pfnProgress( 1.0, NULL, pProgressData );
1478 8 : eErr = CE_None;
1479 :
1480 : end:
1481 8 : VSIFree(pabyPrecomputed);
1482 8 : CPLFree(pafSourceBuf);
1483 8 : CPLFree(panSourceBuf);
1484 8 : CPLFree(pabyDestBuf1);
1485 8 : CPLFree(pasColorAssociation);
1486 :
1487 8 : return eErr;
1488 : }
1489 :
1490 : /************************************************************************/
1491 : /* GDALGenerateVRTColorRelief() */
1492 : /************************************************************************/
1493 :
1494 4 : CPLErr GDALGenerateVRTColorRelief(const char* pszDstFilename,
1495 : GDALDatasetH hSrcDataset,
1496 : GDALRasterBandH hSrcBand,
1497 : const char* pszColorFilename,
1498 : ColorSelectionMode eColorSelectionMode,
1499 : int bAddAlpha)
1500 : {
1501 :
1502 4 : int nColorAssociation = 0;
1503 : ColorAssociation* pasColorAssociation =
1504 : GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
1505 4 : &nColorAssociation);
1506 4 : if (pasColorAssociation == NULL)
1507 0 : return CE_Failure;
1508 :
1509 4 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
1510 4 : int nYSize = GDALGetRasterBandYSize(hSrcBand);
1511 :
1512 4 : VSILFILE* fp = VSIFOpenL(pszDstFilename, "wt");
1513 4 : if (fp == NULL)
1514 : {
1515 0 : CPLFree(pasColorAssociation);
1516 0 : return CE_Failure;
1517 : }
1518 :
1519 4 : VSIFPrintfL(fp, "<VRTDataset rasterXSize=\"%d\" rasterYSize=\"%d\">\n", nXSize, nYSize);
1520 4 : const char* pszProjectionRef = GDALGetProjectionRef(hSrcDataset);
1521 4 : if (pszProjectionRef && pszProjectionRef[0] != '\0')
1522 : {
1523 4 : char* pszEscapedString = CPLEscapeString(pszProjectionRef, -1, CPLES_XML);
1524 4 : VSIFPrintfL(fp, " <SRS>%s</SRS>\n", pszEscapedString);
1525 4 : VSIFree(pszEscapedString);
1526 : }
1527 : double adfGT[6];
1528 4 : if (GDALGetGeoTransform(hSrcDataset, adfGT) == CE_None)
1529 : {
1530 : VSIFPrintfL(fp, " <GeoTransform> %.16g, %.16g, %.16g, "
1531 : "%.16g, %.16g, %.16g</GeoTransform>\n",
1532 4 : adfGT[0], adfGT[1], adfGT[2], adfGT[3], adfGT[4], adfGT[5]);
1533 : }
1534 4 : int nBands = 3 + (bAddAlpha ? 1 : 0);
1535 : int iBand;
1536 :
1537 : int nBlockXSize, nBlockYSize;
1538 4 : GDALGetBlockSize(hSrcBand, &nBlockXSize, &nBlockYSize);
1539 :
1540 : int bRelativeToVRT;
1541 4 : CPLString osPath = CPLGetPath(pszDstFilename);
1542 : char* pszSourceFilename = CPLStrdup(
1543 : CPLExtractRelativePath( osPath.c_str(), GDALGetDescription(hSrcDataset),
1544 4 : &bRelativeToVRT ));
1545 :
1546 16 : for(iBand = 0; iBand < nBands; iBand++)
1547 : {
1548 12 : VSIFPrintfL(fp, " <VRTRasterBand dataType=\"Byte\" band=\"%d\">\n", iBand + 1);
1549 : VSIFPrintfL(fp, " <ColorInterp>%s</ColorInterp>\n",
1550 12 : GDALGetColorInterpretationName((GDALColorInterp)(GCI_RedBand + iBand)));
1551 12 : VSIFPrintfL(fp, " <ComplexSource>\n");
1552 : VSIFPrintfL(fp, " <SourceFilename relativeToVRT=\"%d\">%s</SourceFilename>\n",
1553 12 : bRelativeToVRT, pszSourceFilename);
1554 12 : 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 12 : nBlockXSize, nBlockYSize);
1561 : VSIFPrintfL(fp, " <SrcRect xOff=\"0\" yOff=\"0\" xSize=\"%d\" ySize=\"%d\"/>\n",
1562 12 : nXSize, nYSize);
1563 : VSIFPrintfL(fp, " <DstRect xOff=\"0\" yOff=\"0\" xSize=\"%d\" ySize=\"%d\"/>\n",
1564 12 : nXSize, nYSize);
1565 :
1566 12 : VSIFPrintfL(fp, " <LUT>");
1567 : int iColor;
1568 : #define EPSILON 1e-8
1569 96 : for(iColor=0;iColor<nColorAssociation;iColor++)
1570 : {
1571 84 : if (eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY)
1572 : {
1573 42 : if (iColor > 1)
1574 30 : VSIFPrintfL(fp, ",");
1575 : }
1576 42 : else if (iColor > 0)
1577 36 : VSIFPrintfL(fp, ",");
1578 :
1579 84 : double dfVal = pasColorAssociation[iColor].dfVal;
1580 :
1581 84 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
1582 : {
1583 0 : VSIFPrintfL(fp, "%.12g:0,", dfVal - EPSILON);
1584 : }
1585 84 : else if (iColor > 0 &&
1586 : eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY)
1587 : {
1588 36 : double dfMidVal = (dfVal + pasColorAssociation[iColor-1].dfVal) / 2;
1589 : VSIFPrintfL(fp, "%.12g:%d", dfMidVal - EPSILON,
1590 12 : (iBand == 0) ? pasColorAssociation[iColor-1].nR :
1591 12 : (iBand == 1) ? pasColorAssociation[iColor-1].nG :
1592 12 : (iBand == 2) ? pasColorAssociation[iColor-1].nB :
1593 72 : pasColorAssociation[iColor-1].nA);
1594 : VSIFPrintfL(fp, ",%.12g:%d", dfMidVal ,
1595 12 : (iBand == 0) ? pasColorAssociation[iColor].nR :
1596 12 : (iBand == 1) ? pasColorAssociation[iColor].nG :
1597 12 : (iBand == 2) ? pasColorAssociation[iColor].nB :
1598 72 : pasColorAssociation[iColor].nA);
1599 :
1600 : }
1601 :
1602 84 : if (eColorSelectionMode != COLOR_SELECTION_NEAREST_ENTRY)
1603 : {
1604 42 : if (dfVal != (double)(int)dfVal)
1605 0 : VSIFPrintfL(fp, "%.12g", dfVal);
1606 : else
1607 42 : VSIFPrintfL(fp, "%d", (int)dfVal);
1608 : VSIFPrintfL(fp, ":%d",
1609 14 : (iBand == 0) ? pasColorAssociation[iColor].nR :
1610 14 : (iBand == 1) ? pasColorAssociation[iColor].nG :
1611 14 : (iBand == 2) ? pasColorAssociation[iColor].nB :
1612 84 : pasColorAssociation[iColor].nA);
1613 : }
1614 :
1615 84 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
1616 : {
1617 0 : VSIFPrintfL(fp, ",%.12g:0", dfVal + EPSILON);
1618 : }
1619 :
1620 : }
1621 12 : VSIFPrintfL(fp, "</LUT>\n");
1622 :
1623 12 : VSIFPrintfL(fp, " </ComplexSource>\n");
1624 12 : VSIFPrintfL(fp, " </VRTRasterBand>\n");
1625 : }
1626 :
1627 4 : CPLFree(pszSourceFilename);
1628 :
1629 4 : VSIFPrintfL(fp, "</VRTDataset>\n");
1630 :
1631 4 : VSIFCloseL(fp);
1632 :
1633 4 : CPLFree(pasColorAssociation);
1634 :
1635 4 : 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 4 : {
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 4 : GDALGeneric3x3Dataset::GDALGeneric3x3Dataset(
1762 : GDALDatasetH hSrcDS,
1763 : GDALRasterBandH hSrcBand,
1764 : GDALDataType eDstDataType,
1765 : int bDstHasNoData,
1766 : double dfDstNoDataValue,
1767 : GDALGeneric3x3ProcessingAlg pfnAlg,
1768 : void* pAlgData,
1769 4 : int bComputeAtEdges)
1770 : {
1771 4 : this->hSrcDS = hSrcDS;
1772 4 : this->hSrcBand = hSrcBand;
1773 4 : this->pfnAlg = pfnAlg;
1774 4 : this->pAlgData = pAlgData;
1775 4 : this->bDstHasNoData = bDstHasNoData;
1776 4 : this->dfDstNoDataValue = dfDstNoDataValue;
1777 4 : this->bComputeAtEdges = bComputeAtEdges;
1778 :
1779 4 : CPLAssert(eDstDataType == GDT_Byte || eDstDataType == GDT_Float32);
1780 :
1781 4 : nRasterXSize = GDALGetRasterXSize(hSrcDS);
1782 4 : nRasterYSize = GDALGetRasterYSize(hSrcDS);
1783 :
1784 4 : SetBand(1, new GDALGeneric3x3RasterBand(this, eDstDataType));
1785 :
1786 4 : apafSourceBuf[0] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1787 4 : apafSourceBuf[1] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1788 4 : apafSourceBuf[2] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1789 :
1790 4 : nCurLine = -1;
1791 4 : }
1792 :
1793 4 : GDALGeneric3x3Dataset::~GDALGeneric3x3Dataset()
1794 : {
1795 4 : CPLFree(apafSourceBuf[0]);
1796 4 : CPLFree(apafSourceBuf[1]);
1797 4 : CPLFree(apafSourceBuf[2]);
1798 4 : }
1799 :
1800 4 : CPLErr GDALGeneric3x3Dataset::GetGeoTransform( double * padfGeoTransform )
1801 : {
1802 4 : return GDALGetGeoTransform(hSrcDS, padfGeoTransform);
1803 : }
1804 :
1805 4 : const char *GDALGeneric3x3Dataset::GetProjectionRef()
1806 : {
1807 4 : return GDALGetProjectionRef(hSrcDS);
1808 : }
1809 :
1810 4 : GDALGeneric3x3RasterBand::GDALGeneric3x3RasterBand(GDALGeneric3x3Dataset *poDS,
1811 4 : GDALDataType eDstDataType)
1812 : {
1813 4 : this->poDS = poDS;
1814 4 : this->nBand = 1;
1815 4 : eDataType = eDstDataType;
1816 4 : nBlockXSize = poDS->GetRasterXSize();
1817 4 : nBlockYSize = 1;
1818 :
1819 4 : bSrcHasNoData = FALSE;
1820 : fSrcNoDataValue = (float)GDALGetRasterNoDataValue(poDS->hSrcBand,
1821 4 : &bSrcHasNoData);
1822 4 : }
1823 :
1824 4 : void GDALGeneric3x3RasterBand::InitWidthNoData(void* pImage)
1825 : {
1826 : int j;
1827 4 : GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
1828 4 : if (eDataType == GDT_Byte)
1829 : {
1830 488 : for(j=0;j<nBlockXSize;j++)
1831 484 : ((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 4 : }
1839 :
1840 484 : CPLErr GDALGeneric3x3RasterBand::IReadBlock( int nBlockXOff,
1841 : int nBlockYOff,
1842 : void *pImage )
1843 : {
1844 : int i, j;
1845 : float fVal;
1846 484 : GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
1847 :
1848 722 : if (poGDS->bComputeAtEdges && nRasterXSize >= 2 && nRasterYSize >= 2)
1849 : {
1850 242 : if (nBlockYOff == 0)
1851 : {
1852 6 : for(i=0;i<2;i++)
1853 : {
1854 : CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1855 : GF_Read,
1856 : 0, i, nBlockXSize, 1,
1857 4 : poGDS->apafSourceBuf[i+1],
1858 : nBlockXSize, 1,
1859 : GDT_Float32,
1860 8 : 0, 0);
1861 4 : if (eErr != CE_None)
1862 : {
1863 0 : InitWidthNoData(pImage);
1864 0 : return eErr;
1865 : }
1866 : }
1867 2 : poGDS->nCurLine = 0;
1868 :
1869 244 : for (j = 0; j < nRasterXSize; j++)
1870 : {
1871 : float afWin[9];
1872 242 : int jmin = (j == 0) ? j : j - 1;
1873 242 : int jmax = (j == nRasterXSize - 1) ? j : j + 1;
1874 :
1875 242 : afWin[0] = INTERPOL(poGDS->apafSourceBuf[1][jmin], poGDS->apafSourceBuf[2][jmin]);
1876 242 : afWin[1] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[2][j]);
1877 242 : afWin[2] = INTERPOL(poGDS->apafSourceBuf[1][jmax], poGDS->apafSourceBuf[2][jmax]);
1878 242 : afWin[3] = poGDS->apafSourceBuf[1][jmin];
1879 242 : afWin[4] = poGDS->apafSourceBuf[1][j];
1880 242 : afWin[5] = poGDS->apafSourceBuf[1][jmax];
1881 242 : afWin[6] = poGDS->apafSourceBuf[2][jmin];
1882 242 : afWin[7] = poGDS->apafSourceBuf[2][j];
1883 242 : afWin[8] = poGDS->apafSourceBuf[2][jmax];
1884 :
1885 : fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
1886 : afWin, (float) poGDS->dfDstNoDataValue,
1887 : poGDS->pfnAlg,
1888 : poGDS->pAlgData,
1889 242 : poGDS->bComputeAtEdges);
1890 :
1891 242 : if (eDataType == GDT_Byte)
1892 242 : ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
1893 : else
1894 0 : ((float*)pImage)[j] = fVal;
1895 : }
1896 :
1897 2 : return CE_None;
1898 : }
1899 240 : else if (nBlockYOff == nRasterYSize - 1)
1900 : {
1901 2 : 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 244 : for (j = 0; j < nRasterXSize; j++)
1921 : {
1922 : float afWin[9];
1923 242 : int jmin = (j == 0) ? j : j - 1;
1924 242 : int jmax = (j == nRasterXSize - 1) ? j : j + 1;
1925 :
1926 242 : afWin[0] = poGDS->apafSourceBuf[1][jmin];
1927 242 : afWin[1] = poGDS->apafSourceBuf[1][j];
1928 242 : afWin[2] = poGDS->apafSourceBuf[1][jmax];
1929 242 : afWin[3] = poGDS->apafSourceBuf[2][jmin];
1930 242 : afWin[4] = poGDS->apafSourceBuf[2][j];
1931 242 : afWin[5] = poGDS->apafSourceBuf[2][jmax];
1932 242 : afWin[6] = INTERPOL(poGDS->apafSourceBuf[2][jmin], poGDS->apafSourceBuf[1][jmin]);
1933 242 : afWin[7] = INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[1][j]);
1934 242 : 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 242 : poGDS->bComputeAtEdges);
1941 :
1942 242 : if (eDataType == GDT_Byte)
1943 242 : ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
1944 : else
1945 0 : ((float*)pImage)[j] = fVal;
1946 : }
1947 :
1948 2 : return CE_None;
1949 : }
1950 : }
1951 242 : else if ( nBlockYOff == 0 || nBlockYOff == nRasterYSize - 1)
1952 : {
1953 4 : InitWidthNoData(pImage);
1954 4 : return CE_None;
1955 : }
1956 :
1957 476 : if ( poGDS->nCurLine != nBlockYOff )
1958 : {
1959 476 : if (poGDS->nCurLine + 1 == nBlockYOff)
1960 : {
1961 474 : float* pafTmp = poGDS->apafSourceBuf[0];
1962 474 : poGDS->apafSourceBuf[0] = poGDS->apafSourceBuf[1];
1963 474 : poGDS->apafSourceBuf[1] = poGDS->apafSourceBuf[2];
1964 474 : poGDS->apafSourceBuf[2] = pafTmp;
1965 :
1966 : CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1967 : GF_Read,
1968 : 0, nBlockYOff + 1, nBlockXSize, 1,
1969 474 : poGDS->apafSourceBuf[2],
1970 : nBlockXSize, 1,
1971 : GDT_Float32,
1972 948 : 0, 0);
1973 :
1974 474 : if (eErr != CE_None)
1975 : {
1976 0 : InitWidthNoData(pImage);
1977 0 : return eErr;
1978 : }
1979 : }
1980 : else
1981 : {
1982 8 : for(i=0;i<3;i++)
1983 : {
1984 : CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1985 : GF_Read,
1986 : 0, nBlockYOff + i - 1, nBlockXSize, 1,
1987 6 : poGDS->apafSourceBuf[i],
1988 : nBlockXSize, 1,
1989 : GDT_Float32,
1990 12 : 0, 0);
1991 6 : if (eErr != CE_None)
1992 : {
1993 0 : InitWidthNoData(pImage);
1994 0 : return eErr;
1995 : }
1996 : }
1997 : }
1998 :
1999 476 : poGDS->nCurLine = nBlockYOff;
2000 : }
2001 :
2002 714 : if (poGDS->bComputeAtEdges && nRasterXSize >= 2)
2003 : {
2004 : float afWin[9];
2005 :
2006 238 : j = 0;
2007 238 : afWin[0] = INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j+1]);
2008 238 : afWin[1] = poGDS->apafSourceBuf[0][j];
2009 238 : afWin[2] = poGDS->apafSourceBuf[0][j+1];
2010 238 : afWin[3] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j+1]);
2011 238 : afWin[4] = poGDS->apafSourceBuf[1][j];
2012 238 : afWin[5] = poGDS->apafSourceBuf[1][j+1];
2013 238 : afWin[6] = INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j+1]);
2014 238 : afWin[7] = poGDS->apafSourceBuf[2][j];
2015 238 : 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 238 : poGDS->bComputeAtEdges);
2022 :
2023 238 : if (eDataType == GDT_Byte)
2024 238 : ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
2025 : else
2026 0 : ((float*)pImage)[j] = fVal;
2027 :
2028 238 : j = nRasterXSize - 1;
2029 :
2030 238 : afWin[0] = poGDS->apafSourceBuf[0][j-1];
2031 238 : afWin[1] = poGDS->apafSourceBuf[0][j];
2032 238 : afWin[2] = INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j-1]);
2033 238 : afWin[3] = poGDS->apafSourceBuf[1][j-1];
2034 238 : afWin[4] = poGDS->apafSourceBuf[1][j];
2035 238 : afWin[5] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j-1]);
2036 238 : afWin[6] = poGDS->apafSourceBuf[2][j-1];
2037 238 : afWin[7] = poGDS->apafSourceBuf[2][j];
2038 238 : 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 238 : poGDS->bComputeAtEdges);
2045 :
2046 238 : if (eDataType == GDT_Byte)
2047 238 : ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
2048 : else
2049 0 : ((float*)pImage)[j] = fVal;
2050 : }
2051 : else
2052 : {
2053 238 : if (eDataType == GDT_Byte)
2054 : {
2055 238 : ((GByte*)pImage)[0] = (GByte) poGDS->dfDstNoDataValue;
2056 238 : if (nBlockXSize > 1)
2057 238 : ((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 57120 : for(j=1;j<nBlockXSize - 1;j++)
2069 : {
2070 : float afWin[9];
2071 56644 : afWin[0] = poGDS->apafSourceBuf[0][j-1];
2072 56644 : afWin[1] = poGDS->apafSourceBuf[0][j];
2073 56644 : afWin[2] = poGDS->apafSourceBuf[0][j+1];
2074 56644 : afWin[3] = poGDS->apafSourceBuf[1][j-1];
2075 56644 : afWin[4] = poGDS->apafSourceBuf[1][j];
2076 56644 : afWin[5] = poGDS->apafSourceBuf[1][j+1];
2077 56644 : afWin[6] = poGDS->apafSourceBuf[2][j-1];
2078 56644 : afWin[7] = poGDS->apafSourceBuf[2][j];
2079 56644 : 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 56644 : poGDS->bComputeAtEdges);
2086 :
2087 56644 : if (eDataType == GDT_Byte)
2088 56644 : ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
2089 : else
2090 0 : ((float*)pImage)[j] = fVal;
2091 :
2092 : }
2093 :
2094 476 : return CE_None;
2095 : }
2096 :
2097 16 : double GDALGeneric3x3RasterBand::GetNoDataValue( int* pbHasNoData )
2098 : {
2099 16 : GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
2100 16 : if (pbHasNoData)
2101 12 : *pbHasNoData = poGDS->bDstHasNoData;
2102 16 : 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 32 : int main( int argc, char ** argv )
2121 :
2122 : {
2123 32 : Algorithm eUtilityMode = HILL_SHADE;
2124 32 : double z = 1.0;
2125 32 : double scale = 1.0;
2126 32 : double az = 315.0;
2127 32 : double alt = 45.0;
2128 : // 0 = 'percent' or 1 = 'degrees'
2129 32 : int slopeFormat = 1;
2130 32 : int bAddAlpha = FALSE;
2131 32 : int bZeroForFlat = FALSE;
2132 32 : int bAngleAsAzimuth = TRUE;
2133 32 : ColorSelectionMode eColorSelectionMode = COLOR_SELECTION_INTERPOLATE;
2134 :
2135 32 : int nBand = 1;
2136 : double adfGeoTransform[6];
2137 :
2138 32 : const char *pszSrcFilename = NULL;
2139 32 : const char *pszDstFilename = NULL;
2140 32 : const char *pszColorFilename = NULL;
2141 32 : const char *pszFormat = "GTiff";
2142 32 : int bFormatExplicitelySet = FALSE;
2143 32 : char **papszCreateOptions = NULL;
2144 :
2145 32 : GDALDatasetH hSrcDataset = NULL;
2146 32 : GDALDatasetH hDstDataset = NULL;
2147 32 : GDALRasterBandH hSrcBand = NULL;
2148 32 : GDALRasterBandH hDstBand = NULL;
2149 32 : GDALDriverH hDriver = NULL;
2150 :
2151 32 : GDALProgressFunc pfnProgress = GDALTermProgress;
2152 :
2153 32 : int nXSize = 0;
2154 32 : int nYSize = 0;
2155 :
2156 32 : int bComputeAtEdges = FALSE;
2157 32 : int bZevenbergenThorne = FALSE;
2158 :
2159 32 : int bQuiet = FALSE;
2160 :
2161 : /* Check strict compilation and runtime library version as we use C++ API */
2162 32 : if (! GDAL_CHECK_VERSION(argv[0]))
2163 0 : exit(1);
2164 :
2165 32 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
2166 32 : if( argc < 2 )
2167 : {
2168 0 : fprintf(stderr, "Not enough arguments\n");
2169 0 : Usage();
2170 : }
2171 :
2172 32 : 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 2 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
2176 2 : return 0;
2177 : }
2178 40 : else if ( EQUAL(argv[1], "shade") || EQUAL(argv[1], "hillshade") )
2179 : {
2180 10 : eUtilityMode = HILL_SHADE;
2181 : }
2182 20 : else if ( EQUAL(argv[1], "slope") )
2183 : {
2184 2 : eUtilityMode = SLOPE;
2185 : }
2186 18 : else if ( EQUAL(argv[1], "aspect") )
2187 : {
2188 2 : eUtilityMode = ASPECT;
2189 : }
2190 16 : else if ( EQUAL(argv[1], "color-relief") )
2191 : {
2192 16 : 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 152 : for(int i = 2; i < argc; i++ )
2216 : {
2217 216 : if( eUtilityMode == HILL_SHADE && i + 1 < argc &&
2218 84 : (EQUAL(argv[i], "--z") || EQUAL(argv[i], "-z")))
2219 : {
2220 10 : z = atof(argv[++i]);
2221 : }
2222 112 : else if ( eUtilityMode == SLOPE && EQUAL(argv[i], "-p"))
2223 : {
2224 0 : slopeFormat = 0;
2225 : }
2226 164 : else if ( (eUtilityMode == HILL_SHADE || eUtilityMode == SLOPE ||
2227 52 : 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 112 : else if ( eUtilityMode == ASPECT && EQUAL(argv[i], "-trigonometric"))
2239 : {
2240 0 : bAngleAsAzimuth = FALSE;
2241 : }
2242 112 : else if ( eUtilityMode == ASPECT && EQUAL(argv[i], "-zero_for_flat"))
2243 : {
2244 0 : bZeroForFlat = TRUE;
2245 : }
2246 112 : else if ( eUtilityMode == COLOR_RELIEF && EQUAL(argv[i], "-exact_color_entry"))
2247 : {
2248 0 : eColorSelectionMode = COLOR_SELECTION_EXACT_ENTRY;
2249 : }
2250 116 : else if ( eUtilityMode == COLOR_RELIEF && EQUAL(argv[i], "-nearest_color_entry"))
2251 : {
2252 4 : eColorSelectionMode = COLOR_SELECTION_NEAREST_ENTRY;
2253 : }
2254 408 : else if( i + 1 < argc &&
2255 78 : (EQUAL(argv[i], "--s") ||
2256 78 : EQUAL(argv[i], "-s") ||
2257 66 : EQUAL(argv[i], "--scale") ||
2258 66 : EQUAL(argv[i], "-scale"))
2259 : )
2260 : {
2261 12 : scale = atof(argv[++i]);
2262 : }
2263 182 : else if( eUtilityMode == HILL_SHADE && i + 1 < argc &&
2264 22 : (EQUAL(argv[i], "--az") ||
2265 22 : EQUAL(argv[i], "-az") ||
2266 20 : EQUAL(argv[i], "--azimuth") ||
2267 20 : EQUAL(argv[i], "-azimuth"))
2268 : )
2269 : {
2270 2 : az = atof(argv[++i]);
2271 : }
2272 174 : else if( eUtilityMode == HILL_SHADE && i + 1 < argc &&
2273 20 : (EQUAL(argv[i], "--alt") ||
2274 20 : EQUAL(argv[i], "-alt") ||
2275 20 : EQUAL(argv[i], "--alt") ||
2276 20 : EQUAL(argv[i], "-alt"))
2277 : )
2278 : {
2279 0 : alt = atof(argv[++i]);
2280 : }
2281 150 : else if( eUtilityMode == COLOR_RELIEF &&
2282 56 : EQUAL(argv[i], "-alpha"))
2283 : {
2284 0 : bAddAlpha = TRUE;
2285 : }
2286 136 : else if( eUtilityMode != COLOR_RELIEF &&
2287 38 : EQUAL(argv[i], "-compute_edges"))
2288 : {
2289 4 : bComputeAtEdges = TRUE;
2290 : }
2291 210 : else if( i + 1 < argc &&
2292 60 : (EQUAL(argv[i], "--b") ||
2293 60 : EQUAL(argv[i], "-b"))
2294 : )
2295 : {
2296 0 : nBand = atoi(argv[++i]);
2297 : }
2298 90 : else if ( EQUAL(argv[i], "-q") || EQUAL(argv[i], "-quiet") )
2299 : {
2300 0 : pfnProgress = GDALDummyProgress;
2301 0 : bQuiet = TRUE;
2302 : }
2303 92 : else if( EQUAL(argv[i],"-co") && i < argc-1 )
2304 : {
2305 2 : papszCreateOptions = CSLAddString( papszCreateOptions, argv[++i] );
2306 : }
2307 100 : else if( EQUAL(argv[i],"-of") && i < argc-1 )
2308 : {
2309 12 : pszFormat = argv[++i];
2310 12 : bFormatExplicitelySet = TRUE;
2311 : }
2312 76 : 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 76 : else if( pszSrcFilename == NULL )
2319 : {
2320 30 : pszSrcFilename = argv[i];
2321 : }
2322 62 : else if( eUtilityMode == COLOR_RELIEF && pszColorFilename == NULL )
2323 : {
2324 16 : pszColorFilename = argv[i];
2325 : }
2326 30 : else if( pszDstFilename == NULL )
2327 : {
2328 30 : pszDstFilename = argv[i];
2329 : }
2330 : else
2331 0 : Usage();
2332 : }
2333 :
2334 30 : if( pszSrcFilename == NULL )
2335 : {
2336 0 : fprintf( stderr, "Missing source.\n\n" );
2337 0 : Usage();
2338 : }
2339 30 : if ( eUtilityMode == COLOR_RELIEF && pszColorFilename == NULL )
2340 : {
2341 0 : fprintf( stderr, "Missing color file.\n\n" );
2342 0 : Usage();
2343 : }
2344 30 : if( pszDstFilename == NULL )
2345 : {
2346 0 : fprintf( stderr, "Missing destination.\n\n" );
2347 0 : Usage();
2348 : }
2349 :
2350 30 : GDALAllRegister();
2351 :
2352 : // Open Dataset and get raster band
2353 30 : hSrcDataset = GDALOpen( pszSrcFilename, GA_ReadOnly );
2354 :
2355 30 : 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 30 : nXSize = GDALGetRasterXSize(hSrcDataset);
2365 30 : nYSize = GDALGetRasterYSize(hSrcDataset);
2366 :
2367 30 : 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 30 : hSrcBand = GDALGetRasterBand( hSrcDataset, nBand );
2375 :
2376 30 : GDALGetGeoTransform(hSrcDataset, adfGeoTransform);
2377 :
2378 30 : hDriver = GDALGetDriverByName(pszFormat);
2379 30 : 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 30 : if (!bQuiet && !bFormatExplicitelySet)
2407 18 : CheckExtensionConsistency(pszDstFilename, pszFormat);
2408 :
2409 30 : double dfDstNoDataValue = 0;
2410 30 : int bDstHasNoData = FALSE;
2411 30 : void* pData = NULL;
2412 30 : GDALGeneric3x3ProcessingAlg pfnAlg = NULL;
2413 :
2414 30 : if (eUtilityMode == HILL_SHADE)
2415 : {
2416 10 : dfDstNoDataValue = 0;
2417 10 : bDstHasNoData = TRUE;
2418 : pData = GDALCreateHillshadeData (adfGeoTransform,
2419 : z,
2420 : scale,
2421 : alt,
2422 : az,
2423 10 : bZevenbergenThorne);
2424 10 : if (bZevenbergenThorne)
2425 0 : pfnAlg = GDALHillshadeZevenbergenThorneAlg;
2426 : else
2427 10 : pfnAlg = GDALHillshadeAlg;
2428 : }
2429 20 : else if (eUtilityMode == SLOPE)
2430 : {
2431 2 : dfDstNoDataValue = -9999;
2432 2 : bDstHasNoData = TRUE;
2433 :
2434 2 : pData = GDALCreateSlopeData(adfGeoTransform, scale, slopeFormat);
2435 2 : if (bZevenbergenThorne)
2436 0 : pfnAlg = GDALSlopeZevenbergenThorneAlg;
2437 : else
2438 2 : pfnAlg = GDALSlopeHornAlg;
2439 : }
2440 :
2441 18 : else if (eUtilityMode == ASPECT)
2442 : {
2443 2 : if (!bZeroForFlat)
2444 : {
2445 2 : dfDstNoDataValue = -9999;
2446 2 : bDstHasNoData = TRUE;
2447 : }
2448 :
2449 2 : pData = GDALCreateAspectData(bAngleAsAzimuth);
2450 2 : if (bZevenbergenThorne)
2451 0 : pfnAlg = GDALAspectZevenbergenThorneAlg;
2452 : else
2453 2 : pfnAlg = GDALAspectAlg;
2454 : }
2455 16 : else if (eUtilityMode == TRI)
2456 : {
2457 0 : dfDstNoDataValue = -9999;
2458 0 : bDstHasNoData = TRUE;
2459 0 : pfnAlg = GDALTRIAlg;
2460 : }
2461 16 : else if (eUtilityMode == TPI)
2462 : {
2463 0 : dfDstNoDataValue = -9999;
2464 0 : bDstHasNoData = TRUE;
2465 0 : pfnAlg = GDALTPIAlg;
2466 : }
2467 16 : 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 30 : GDT_Float32;
2477 :
2478 30 : if( EQUAL(pszFormat, "VRT") )
2479 : {
2480 4 : if (eUtilityMode == COLOR_RELIEF)
2481 : {
2482 : GDALGenerateVRTColorRelief(pszDstFilename,
2483 : hSrcDataset,
2484 : hSrcBand,
2485 : pszColorFilename,
2486 : eColorSelectionMode,
2487 4 : bAddAlpha);
2488 4 : GDALClose(hSrcDataset);
2489 :
2490 4 : CPLFree(pData);
2491 :
2492 4 : GDALDestroyDriverManager();
2493 4 : CSLDestroy( argv );
2494 4 : CSLDestroy( papszCreateOptions );
2495 :
2496 4 : 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 26 : if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL &&
2508 : GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, NULL ) != NULL)
2509 : {
2510 : GDALDatasetH hIntermediateDataset;
2511 :
2512 8 : if (eUtilityMode == COLOR_RELIEF)
2513 : hIntermediateDataset = (GDALDatasetH)
2514 : new GDALColorReliefDataset (hSrcDataset,
2515 : hSrcBand,
2516 : pszColorFilename,
2517 : eColorSelectionMode,
2518 4 : bAddAlpha);
2519 : else
2520 : hIntermediateDataset = (GDALDatasetH)
2521 : new GDALGeneric3x3Dataset(hSrcDataset, hSrcBand,
2522 : eDstDataType,
2523 : bDstHasNoData,
2524 : dfDstNoDataValue,
2525 : pfnAlg,
2526 : pData,
2527 4 : bComputeAtEdges);
2528 :
2529 : GDALDatasetH hOutDS = GDALCreateCopy(
2530 : hDriver, pszDstFilename, hIntermediateDataset,
2531 : TRUE, papszCreateOptions,
2532 8 : pfnProgress, NULL );
2533 :
2534 8 : if( hOutDS != NULL )
2535 8 : GDALClose( hOutDS );
2536 8 : GDALClose(hIntermediateDataset);
2537 8 : GDALClose(hSrcDataset);
2538 :
2539 8 : CPLFree(pData);
2540 :
2541 8 : GDALDestroyDriverManager();
2542 8 : CSLDestroy( argv );
2543 8 : CSLDestroy( papszCreateOptions );
2544 :
2545 8 : return 0;
2546 : }
2547 :
2548 : int nDstBands;
2549 18 : if (eUtilityMode == COLOR_RELIEF)
2550 8 : nDstBands = (bAddAlpha) ? 4 : 3;
2551 : else
2552 10 : nDstBands = 1;
2553 :
2554 : hDstDataset = GDALCreate( hDriver,
2555 : pszDstFilename,
2556 : nXSize,
2557 : nYSize,
2558 : nDstBands,
2559 : eDstDataType,
2560 18 : papszCreateOptions);
2561 :
2562 18 : 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 18 : hDstBand = GDALGetRasterBand( hDstDataset, 1 );
2572 :
2573 18 : GDALSetGeoTransform(hDstDataset, adfGeoTransform);
2574 18 : GDALSetProjection(hDstDataset, GDALGetProjectionRef(hSrcDataset));
2575 :
2576 18 : 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 8 : pfnProgress, NULL);
2586 : }
2587 : else
2588 : {
2589 10 : if (bDstHasNoData)
2590 10 : GDALSetRasterNoDataValue(hDstBand, dfDstNoDataValue);
2591 :
2592 : GDALGeneric3x3Processing(hSrcBand, hDstBand,
2593 : pfnAlg, pData,
2594 : bComputeAtEdges,
2595 10 : pfnProgress, NULL);
2596 :
2597 : }
2598 :
2599 18 : GDALClose(hSrcDataset);
2600 18 : GDALClose(hDstDataset);
2601 18 : CPLFree(pData);
2602 :
2603 18 : GDALDestroyDriverManager();
2604 18 : CSLDestroy( argv );
2605 18 : CSLDestroy( papszCreateOptions );
2606 :
2607 18 : return 0;
2608 : }
2609 :
|
