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