source: trunk/zoo-project/zoo-services/gdal/dem/service.c

Last change on this file was 917, checked in by djay, 10 months ago

Merge prototype-v0 branch in trunk

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

Search

Context Navigation

ZOO Sponsors

http://www.zoo-project.org/trac/chrome/site/img/geolabs-logo.pnghttp://www.zoo-project.org/trac/chrome/site/img/neogeo-logo.png http://www.zoo-project.org/trac/chrome/site/img/apptech-logo.png http://www.zoo-project.org/trac/chrome/site/img/3liz-logo.png http://www.zoo-project.org/trac/chrome/site/img/gateway-logo.png

Become a sponsor !

Knowledge partners

http://www.zoo-project.org/trac/chrome/site/img/ocu-logo.png http://www.zoo-project.org/trac/chrome/site/img/gucas-logo.png http://www.zoo-project.org/trac/chrome/site/img/polimi-logo.png http://www.zoo-project.org/trac/chrome/site/img/fem-logo.png http://www.zoo-project.org/trac/chrome/site/img/supsi-logo.png http://www.zoo-project.org/trac/chrome/site/img/cumtb-logo.png

Become a knowledge partner

Related links

http://zoo-project.org/img/ogclogo.png http://zoo-project.org/img/osgeologo.png