← Back to branch summary

~ubuntu-branches/debian/sid/gdal/sid

« back to all changes in this revision

Viewing changes to apps/gdaldem.cpp

  • Committer: Package Import Robot
  • Author(s): Francesco Paolo Lovergine
  • Date: 2012-05-07 15:04:42 UTC
  • mfrom: (5.5.16 experimental)
  • Revision ID: package-import@ubuntu.com-20120507150442-2eks97loeh6rq005
Tags: 1.9.0-1
* Ready for sid, starting transition.
* All symfiles updated to latest builds.
* Added dh_numpy call in debian/rules to depend on numpy ABI.
* Policy bumped to 3.9.3, no changes required.

Show diffs side-by-side

added added

removed removed

Lines of Context:
apps/gdaldem.cpp
1
1
/******************************************************************************
2
 
 * $Id: gdaldem.cpp 19777 2010-05-27 19:08:13Z rouault $
 
2
 * $Id: gdaldem.cpp 22857 2011-08-02 18:17:45Z rouault $
3
3
 *
4
4
 * Project:  GDAL DEM Utilities
5
5
 * Purpose:  
89
89
#include "cpl_string.h"
90
90
#include "gdal.h"
91
91
#include "gdal_priv.h"
 
92
#include "commonutils.h"
92
93
 
93
 
CPL_CVSID("$Id: gdaldem.cpp 19777 2010-05-27 19:08:13Z rouault $");
 
94
CPL_CVSID("$Id: gdaldem.cpp 22857 2011-08-02 18:17:45Z rouault $");
94
95
 
95
96
#ifndef M_PI
96
97
# define M_PI  3.1415926535897932384626433832795
97
98
#endif
98
99
 
 
100
#define INTERPOL(a,b) ((bSrcHasNoData && (ARE_REAL_EQUAL(a, fSrcNoDataValue) || ARE_REAL_EQUAL(b, fSrcNoDataValue))) ? fSrcNoDataValue : 2 * (a) - (b))
 
101
 
99
102
/************************************************************************/
100
103
/*                               Usage()                                */
101
104
/************************************************************************/
108
111
            "     gdaldem hillshade input_dem output_hillshade \n"
109
112
            "                 [-z ZFactor (default=1)] [-s scale* (default=1)] \n"
110
113
            "                 [-az Azimuth (default=315)] [-alt Altitude (default=45)]\n"
111
 
            "                 [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
 
114
            "                 [-alg ZevenbergenThorne]\n"
 
115
            "                 [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
112
116
            "\n"
113
117
            " - To generates a slope map from any GDAL-supported elevation raster :\n\n"
114
118
            "     gdaldem slope input_dem output_slope_map \n"
115
119
            "                 [-p use percent slope (default=degrees)] [-s scale* (default=1)]\n"
116
 
            "                 [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
 
120
            "                 [-alg ZevenbergenThorne]\n"
 
121
            "                 [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
117
122
            "\n"
118
123
            " - To generate an aspect map from any GDAL-supported elevation raster\n"
119
124
            "   Outputs a 32-bit float tiff with pixel values from 0-360 indicating azimuth :\n\n"
120
125
            "     gdaldem aspect input_dem output_aspect_map \n"
121
126
            "                 [-trigonometric] [-zero_for_flat]\n"
122
 
            "                 [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
 
127
            "                 [-alg ZevenbergenThorne]\n"
 
128
            "                 [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
123
129
            "\n"
124
130
            " - To generate a color relief map from any GDAL-supported elevation raster\n"
125
131
            "     gdaldem color-relief input_dem color_text_file output_color_relief_map\n"
129
135
            "\n"
130
136
            " - To generate a Terrain Ruggedness Index (TRI) map from any GDAL-supported elevation raster\n"
131
137
            "     gdaldem TRI input_dem output_TRI_map\n"
132
 
            "                 [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
 
138
            "                 [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
133
139
            "\n"
134
140
            " - To generate a Topographic Position Index (TPI) map from any GDAL-supported elevation raster\n"
135
141
            "     gdaldem TPI input_dem output_TPI_map\n"
136
 
            "                 [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
 
142
            "                 [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
137
143
            "\n"
138
144
            " - To generate a roughness map from any GDAL-supported elevation raster\n"
139
145
            "     gdaldem roughness input_dem output_roughness_map\n"
140
 
            "                 [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
 
146
            "                 [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
141
147
            "\n"
142
148
            " Notes : \n"
143
149
            "   Scale is the ratio of vertical units to horizontal\n"
146
152
}
147
153
 
148
154
/************************************************************************/
 
155
/*                          ComputeVal()                                */
 
156
/************************************************************************/
 
157
 
 
158
typedef float (*GDALGeneric3x3ProcessingAlg) (float* pafWindow, float fDstNoDataValue, void* pData);
 
159
 
 
160
static float ComputeVal(int bSrcHasNoData, float fSrcNoDataValue,
 
161
                        float* afWin, float fDstNoDataValue,
 
162
                        GDALGeneric3x3ProcessingAlg pfnAlg,
 
163
                        void* pData,
 
164
                        int bComputeAtEdges)
 
165
{
 
166
    if (bSrcHasNoData && ARE_REAL_EQUAL(afWin[4], fSrcNoDataValue))
 
167
    {
 
168
        return fDstNoDataValue;
 
169
    }
 
170
    else if (bSrcHasNoData)
 
171
    {
 
172
        int k;
 
173
        for(k=0;k<9;k++)
 
174
        {
 
175
            if (ARE_REAL_EQUAL(afWin[k], fSrcNoDataValue))
 
176
            {
 
177
                if (bComputeAtEdges)
 
178
                    afWin[k] = afWin[4];
 
179
                else
 
180
                    return fDstNoDataValue;
 
181
            }
 
182
        }
 
183
    }
 
184
 
 
185
    return pfnAlg(afWin, fDstNoDataValue, pData);
 
186
}
 
187
 
 
188
/************************************************************************/
149
189
/*                  GDALGeneric3x3Processing()                          */
150
190
/************************************************************************/
151
191
 
152
 
typedef float (*GDALGeneric3x3ProcessingAlg) (float* pafWindow, float fDstNoDataValue, void* pData);
153
 
 
154
192
CPLErr GDALGeneric3x3Processing  ( GDALRasterBandH hSrcBand,
155
193
                                   GDALRasterBandH hDstBand,
156
194
                                   GDALGeneric3x3ProcessingAlg pfnAlg,
157
195
                                   void* pData,
 
196
                                   int bComputeAtEdges,
158
197
                                   GDALProgressFunc pfnProgress,
159
198
                                   void * pProgressData)
160
199
{
209
248
                        0, 0);
210
249
    }
211
250
    
212
 
    // Exclude the edges
213
 
    for (j = 0; j < nXSize; j++)
 
251
    if (bComputeAtEdges && nXSize >= 2 && nYSize >= 2)
214
252
    {
215
 
        pafOutputBuf[j] = fDstNoDataValue;
 
253
        for (j = 0; j < nXSize; j++)
 
254
        {
 
255
            float afWin[9];
 
256
            int jmin = (j == 0) ? j : j - 1;
 
257
            int jmax = (j == nXSize - 1) ? j : j + 1;
 
258
 
 
259
            afWin[0] = INTERPOL(pafThreeLineWin[jmin], pafThreeLineWin[nXSize + jmin]);
 
260
            afWin[1] = INTERPOL(pafThreeLineWin[j],    pafThreeLineWin[nXSize + j]);
 
261
            afWin[2] = INTERPOL(pafThreeLineWin[jmax], pafThreeLineWin[nXSize + jmax]);
 
262
            afWin[3] = pafThreeLineWin[jmin];
 
263
            afWin[4] = pafThreeLineWin[j];
 
264
            afWin[5] = pafThreeLineWin[jmax];
 
265
            afWin[6] = pafThreeLineWin[nXSize + jmin];
 
266
            afWin[7] = pafThreeLineWin[nXSize + j];
 
267
            afWin[8] = pafThreeLineWin[nXSize + jmax];
 
268
 
 
269
            pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
 
270
                                         afWin, fDstNoDataValue,
 
271
                                         pfnAlg, pData, bComputeAtEdges);
 
272
        }
 
273
        GDALRasterIO(hDstBand, GF_Write,
 
274
                    0, 0, nXSize, 1,
 
275
                    pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
216
276
    }
217
 
    GDALRasterIO(hDstBand, GF_Write,
218
 
                 0, 0, nXSize, 1,
219
 
                 pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
220
 
 
221
 
    if (nYSize > 1)
 
277
    else
222
278
    {
 
279
        // Exclude the edges
 
280
        for (j = 0; j < nXSize; j++)
 
281
        {
 
282
            pafOutputBuf[j] = fDstNoDataValue;
 
283
        }
223
284
        GDALRasterIO(hDstBand, GF_Write,
224
 
                     0, nYSize - 1, nXSize, 1,
225
 
                     pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
 
285
                    0, 0, nXSize, 1,
 
286
                    pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
 
287
    
 
288
        if (nYSize > 1)
 
289
        {
 
290
            GDALRasterIO(hDstBand, GF_Write,
 
291
                        0, nYSize - 1, nXSize, 1,
 
292
                        pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
 
293
        }
226
294
    }
227
295
    
228
296
    int nLine1Off = 0*nXSize;
243
311
        if (eErr != CE_None)
244
312
            goto end;
245
313
 
246
 
        // Exclude the edges
247
 
        pafOutputBuf[0] = fDstNoDataValue;
248
 
        if (nXSize > 1)
249
 
            pafOutputBuf[nXSize - 1] = fDstNoDataValue;
250
 
        
 
314
        if (bComputeAtEdges && nXSize >= 2)
 
315
        {
 
316
            float afWin[9];
 
317
 
 
318
            j = 0;
 
319
            afWin[0] = INTERPOL(pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j+1]);
 
320
            afWin[1] = pafThreeLineWin[nLine1Off + j];
 
321
            afWin[2] = pafThreeLineWin[nLine1Off + j+1];
 
322
            afWin[3] = INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j+1]);
 
323
            afWin[4] = pafThreeLineWin[nLine2Off + j];
 
324
            afWin[5] = pafThreeLineWin[nLine2Off + j+1];
 
325
            afWin[6] = INTERPOL(pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j+1]);
 
326
            afWin[7] = pafThreeLineWin[nLine3Off + j];
 
327
            afWin[8] = pafThreeLineWin[nLine3Off + j+1];
 
328
 
 
329
            pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
 
330
                                         afWin, fDstNoDataValue,
 
331
                                         pfnAlg, pData, bComputeAtEdges);
 
332
            j = nXSize - 1;
 
333
 
 
334
            afWin[0] = pafThreeLineWin[nLine1Off + j-1];
 
335
            afWin[1] = pafThreeLineWin[nLine1Off + j];
 
336
            afWin[2] = INTERPOL(pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j-1]);
 
337
            afWin[3] = pafThreeLineWin[nLine2Off + j-1];
 
338
            afWin[4] = pafThreeLineWin[nLine2Off + j];
 
339
            afWin[5] = INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j-1]);
 
340
            afWin[6] = pafThreeLineWin[nLine3Off + j-1];
 
341
            afWin[7] = pafThreeLineWin[nLine3Off + j];
 
342
            afWin[8] = INTERPOL(pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j-1]);
 
343
 
 
344
            pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
 
345
                                         afWin, fDstNoDataValue,
 
346
                                         pfnAlg, pData, bComputeAtEdges);
 
347
        }
 
348
        else
 
349
        {
 
350
            // Exclude the edges
 
351
            pafOutputBuf[0] = fDstNoDataValue;
 
352
            if (nXSize > 1)
 
353
                pafOutputBuf[nXSize - 1] = fDstNoDataValue;
 
354
        }
 
355
 
251
356
        for (j = 1; j < nXSize - 1; j++)
252
357
        {
253
358
            float afWin[9];
261
366
            afWin[7] = pafThreeLineWin[nLine3Off + j];
262
367
            afWin[8] = pafThreeLineWin[nLine3Off + j+1];
263
368
 
264
 
            if (bSrcHasNoData && (
265
 
                   afWin[0] == fSrcNoDataValue ||
266
 
                   afWin[1] == fSrcNoDataValue ||
267
 
                   afWin[2] == fSrcNoDataValue ||
268
 
                   afWin[3] == fSrcNoDataValue ||
269
 
                   afWin[4] == fSrcNoDataValue ||
270
 
                   afWin[5] == fSrcNoDataValue ||
271
 
                   afWin[6] == fSrcNoDataValue ||
272
 
                   afWin[7] == fSrcNoDataValue ||
273
 
                   afWin[8] == fSrcNoDataValue))
274
 
            {
275
 
                // We have nulls so write nullValue and move on
276
 
                pafOutputBuf[j] = fDstNoDataValue;
277
 
            }
278
 
            else
279
 
            {
280
 
                // We have a valid 3x3 window.
281
 
                pafOutputBuf[j] = pfnAlg(afWin, fDstNoDataValue, pData);
282
 
            }
 
369
            pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
 
370
                                         afWin, fDstNoDataValue,
 
371
                                         pfnAlg, pData, bComputeAtEdges);
283
372
        }
284
373
 
285
374
        /* -----------------------------------------
303
392
        nLine3Off = nTemp;
304
393
    }
305
394
 
 
395
    if (bComputeAtEdges && nXSize >= 2 && nYSize >= 2)
 
396
    {
 
397
        for (j = 0; j < nXSize; j++)
 
398
        {
 
399
            float afWin[9];
 
400
            int jmin = (j == 0) ? j : j - 1;
 
401
            int jmax = (j == nXSize - 1) ? j : j + 1;
 
402
 
 
403
            afWin[0] = pafThreeLineWin[nLine1Off + jmin];
 
404
            afWin[1] = pafThreeLineWin[nLine1Off + j];
 
405
            afWin[2] = pafThreeLineWin[nLine1Off + jmax];
 
406
            afWin[3] = pafThreeLineWin[nLine2Off + jmin];
 
407
            afWin[4] = pafThreeLineWin[nLine2Off + j];
 
408
            afWin[5] = pafThreeLineWin[nLine2Off + jmax];
 
409
            afWin[6] = INTERPOL(pafThreeLineWin[nLine2Off + jmin], pafThreeLineWin[nLine1Off + jmin]);
 
410
            afWin[7] = INTERPOL(pafThreeLineWin[nLine2Off + j],    pafThreeLineWin[nLine1Off + j]);
 
411
            afWin[8] = INTERPOL(pafThreeLineWin[nLine2Off + jmax], pafThreeLineWin[nLine1Off + jmax]);
 
412
 
 
413
            pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
 
414
                                         afWin, fDstNoDataValue,
 
415
                                         pfnAlg, pData, bComputeAtEdges);
 
416
        }
 
417
        GDALRasterIO(hDstBand, GF_Write,
 
418
                     0, i, nXSize, 1,
 
419
                     pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
 
420
    }
 
421
 
306
422
    pfnProgress( 1.0, NULL, pProgressData );
307
423
    eErr = CE_None;
308
424
 
374
490
    else
375
491
        cang = 1.0 + (254.0 * cang);
376
492
        
377
 
    return cang;
 
493
    return (float) cang;
 
494
}
 
495
 
 
496
float GDALHillshadeZevenbergenThorneAlg (float* afWin, float fDstNoDataValue, void* pData)
 
497
{
 
498
    GDALHillshadeAlgData* psData = (GDALHillshadeAlgData*)pData;
 
499
    double x, y, aspect, xx_plus_yy, cang;
 
500
    
 
501
    // First Slope ...
 
502
    x = (afWin[3] - afWin[5]) / psData->ewres;
 
503
 
 
504
    y = (afWin[7] - afWin[1]) / psData->nsres;
 
505
 
 
506
    xx_plus_yy = x * x + y * y;
 
507
 
 
508
    // ... then aspect...
 
509
    aspect = atan2(y,x);
 
510
 
 
511
    // ... then the shade value
 
512
    cang = (psData->sin_altRadians -
 
513
           psData->cos_altRadians_mul_z_scale_factor * sqrt(xx_plus_yy) *
 
514
           sin(aspect - psData->azRadians)) /
 
515
           sqrt(1 + psData->square_z_scale_factor * xx_plus_yy);
 
516
 
 
517
    if (cang <= 0.0) 
 
518
        cang = 1.0;
 
519
    else
 
520
        cang = 1.0 + (254.0 * cang);
 
521
        
 
522
    return (float) cang;
378
523
}
379
524
 
380
525
void*  GDALCreateHillshadeData(double* adfGeoTransform,
381
526
                               double z,
382
527
                               double scale,
383
528
                               double alt,
384
 
                               double az)
 
529
                               double az,
 
530
                               int bZevenbergenThorne)
385
531
{
386
532
    GDALHillshadeAlgData* pData =
387
533
        (GDALHillshadeAlgData*)CPLMalloc(sizeof(GDALHillshadeAlgData));
391
537
    pData->ewres = adfGeoTransform[1];
392
538
    pData->sin_altRadians = sin(alt * degreesToRadians);
393
539
    pData->azRadians = az * degreesToRadians;
394
 
    double z_scale_factor = z / (8 * scale);
 
540
    double z_scale_factor = z / (((bZevenbergenThorne) ? 2 : 8) * scale);
395
541
    pData->cos_altRadians_mul_z_scale_factor =
396
542
        cos(alt * degreesToRadians) * z_scale_factor;
397
543
    pData->square_z_scale_factor = z_scale_factor * z_scale_factor;
410
556
    int    slopeFormat;
411
557
} GDALSlopeAlgData;
412
558
 
413
 
float GDALSlopeAlg (float* afWin, float fDstNoDataValue, void* pData)
 
559
float GDALSlopeHornAlg (float* afWin, float fDstNoDataValue, void* pData)
414
560
{
415
561
    const double radiansToDegrees = 180.0 / M_PI;
416
562
    GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
425
571
    key = (dx * dx + dy * dy);
426
572
 
427
573
    if (psData->slopeFormat == 1) 
428
 
        return atan(sqrt(key) / (8*psData->scale)) * radiansToDegrees;
429
 
    else
430
 
        return 100*(sqrt(key) / (8*psData->scale));
 
574
        return (float) (atan(sqrt(key) / (8*psData->scale)) * radiansToDegrees);
 
575
    else
 
576
        return (float) (100*(sqrt(key) / (8*psData->scale)));
 
577
}
 
578
 
 
579
float GDALSlopeZevenbergenThorneAlg (float* afWin, float fDstNoDataValue, void* pData)
 
580
{
 
581
    const double radiansToDegrees = 180.0 / M_PI;
 
582
    GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
 
583
    double dx, dy, key;
 
584
    
 
585
    dx = (afWin[3] - afWin[5])/psData->ewres;
 
586
 
 
587
    dy = (afWin[7] - afWin[1])/psData->nsres;
 
588
 
 
589
    key = (dx * dx + dy * dy);
 
590
 
 
591
    if (psData->slopeFormat == 1) 
 
592
        return (float) (atan(sqrt(key) / (2*psData->scale)) * radiansToDegrees);
 
593
    else
 
594
        return (float) (100*(sqrt(key) / (2*psData->scale)));
431
595
}
432
596
 
433
597
void*  GDALCreateSlopeData(double* adfGeoTransform,
466
630
    dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) - 
467
631
          (afWin[0] + afWin[1] + afWin[1] + afWin[2]));
468
632
 
469
 
    aspect = atan2(dy,-dx) / degreesToRadians;
470
 
 
471
 
    if (dx == 0 && dy == 0)
472
 
    {
473
 
        /* Flat area */
474
 
        aspect = fDstNoDataValue;
475
 
    } 
476
 
    else if ( psData->bAngleAsAzimuth )
477
 
    {
478
 
        if (aspect > 90.0) 
479
 
            aspect = 450.0 - aspect;
480
 
        else
481
 
            aspect = 90.0 - aspect;
482
 
    }
483
 
    else
484
 
    {
485
 
        if (aspect < 0)
486
 
            aspect += 360.0;
487
 
    }
488
 
 
489
 
    if (aspect == 360.0) 
490
 
        aspect = 0.0;
491
 
 
492
 
    return aspect;
493
 
}
494
 
 
 
633
    aspect = (float) (atan2(dy,-dx) / degreesToRadians);
 
634
 
 
635
    if (dx == 0 && dy == 0)
 
636
    {
 
637
        /* Flat area */
 
638
        aspect = fDstNoDataValue;
 
639
    } 
 
640
    else if ( psData->bAngleAsAzimuth )
 
641
    {
 
642
        if (aspect > 90.0) 
 
643
            aspect = 450.0f - aspect;
 
644
        else
 
645
            aspect = 90.0f - aspect;
 
646
    }
 
647
    else
 
648
    {
 
649
        if (aspect < 0)
 
650
            aspect += 360.0;
 
651
    }
 
652
 
 
653
    if (aspect == 360.0) 
 
654
        aspect = 0.0;
 
655
 
 
656
    return aspect;
 
657
}
 
658
 
 
659
float GDALAspectZevenbergenThorneAlg (float* afWin, float fDstNoDataValue, void* pData)
 
660
{
 
661
    const double degreesToRadians = M_PI / 180.0;
 
662
    GDALAspectAlgData* psData = (GDALAspectAlgData*)pData;
 
663
    double dx, dy;
 
664
    float aspect;
 
665
    
 
666
    dx = (afWin[5] - afWin[3]);
 
667
 
 
668
    dy = (afWin[7] - afWin[1]);
 
669
 
 
670
    aspect = (float) (atan2(dy,-dx) / degreesToRadians);
 
671
 
 
672
    if (dx == 0 && dy == 0)
 
673
    {
 
674
        /* Flat area */
 
675
        aspect = fDstNoDataValue;
 
676
    } 
 
677
    else if ( psData->bAngleAsAzimuth )
 
678
    {
 
679
        if (aspect > 90.0) 
 
680
            aspect = 450.0f - aspect;
 
681
        else
 
682
            aspect = 90.0f - aspect;
 
683
    }
 
684
    else
 
685
    {
 
686
        if (aspect < 0)
 
687
            aspect += 360.0;
 
688
    }
 
689
 
 
690
    if (aspect == 360.0) 
 
691
        aspect = 0.0;
 
692
 
 
693
    return aspect;
 
694
}
495
695
void*  GDALCreateAspectData(int bAngleAsAzimuth)
496
696
{
497
697
    GDALAspectAlgData* pData =
727
927
                                                const char* pszColorFilename,
728
928
                                                int* pnColors)
729
929
{
730
 
    FILE* fpColorFile = VSIFOpenL(pszColorFilename, "rt");
 
930
    VSILFILE* fpColorFile = VSIFOpenL(pszColorFilename, "rt");
731
931
    if (fpColorFile == NULL)
732
932
    {
733
933
        CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszColorFilename);
742
942
    double dfSrcNoDataValue = GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
743
943
 
744
944
    const char* pszLine;
 
945
    int bIsGMT_CPT = FALSE;
745
946
    while ((pszLine = CPLReadLineL(fpColorFile)) != NULL)
746
947
    {
 
948
        if (pszLine[0] == '#' && strstr(pszLine, "COLOR_MODEL"))
 
949
        {
 
950
            if (strstr(pszLine, "COLOR_MODEL = RGB") == NULL)
 
951
            {
 
952
                CPLError(CE_Failure, CPLE_AppDefined, "Only COLOR_MODEL = RGB is supported");
 
953
                CPLFree(pasColorAssociation);
 
954
                *pnColors = 0;
 
955
                return NULL;
 
956
            }
 
957
            bIsGMT_CPT = TRUE;
 
958
        }
 
959
 
747
960
        char** papszFields = CSLTokenizeStringComplex(pszLine, " ,\t:", 
748
961
                                                      FALSE, FALSE );
749
962
        /* Skip comment lines */
750
963
        int nTokens = CSLCount(papszFields);
751
 
        if (nTokens >= 2 &&
752
 
            papszFields[0][0] != '#' &&
753
 
            papszFields[0][0] != '/')
 
964
        if (nTokens >= 1 && (papszFields[0][0] == '#' ||
 
965
                             papszFields[0][0] == '/'))
 
966
        {
 
967
            CSLDestroy(papszFields);
 
968
            continue;
 
969
        }
 
970
 
 
971
        if (bIsGMT_CPT && nTokens == 8)
 
972
        {
 
973
            pasColorAssociation =
 
974
                    (ColorAssociation*)CPLRealloc(pasColorAssociation,
 
975
                           (nColorAssociation + 2) * sizeof(ColorAssociation));
 
976
 
 
977
            pasColorAssociation[nColorAssociation].dfVal = atof(papszFields[0]);
 
978
            pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]);
 
979
            pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]);
 
980
            pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]);
 
981
            pasColorAssociation[nColorAssociation].nA = 255;
 
982
            nColorAssociation++;
 
983
 
 
984
            pasColorAssociation[nColorAssociation].dfVal = atof(papszFields[4]);
 
985
            pasColorAssociation[nColorAssociation].nR = atoi(papszFields[5]);
 
986
            pasColorAssociation[nColorAssociation].nG = atoi(papszFields[6]);
 
987
            pasColorAssociation[nColorAssociation].nB = atoi(papszFields[7]);
 
988
            pasColorAssociation[nColorAssociation].nA = 255;
 
989
            nColorAssociation++;
 
990
        }
 
991
        else if (bIsGMT_CPT && nTokens == 4)
 
992
        {
 
993
            /* The first token might be B (background), F (foreground) or N (nodata) */
 
994
            /* Just interested in N */
 
995
            if (EQUAL(papszFields[0], "N") && bSrcHasNoData)
 
996
            {
 
997
                 pasColorAssociation =
 
998
                    (ColorAssociation*)CPLRealloc(pasColorAssociation,
 
999
                           (nColorAssociation + 1) * sizeof(ColorAssociation));
 
1000
 
 
1001
                pasColorAssociation[nColorAssociation].dfVal = dfSrcNoDataValue;
 
1002
                pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]);
 
1003
                pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]);
 
1004
                pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]);
 
1005
                pasColorAssociation[nColorAssociation].nA = 255;
 
1006
                nColorAssociation++;
 
1007
            }
 
1008
        }
 
1009
        else if (!bIsGMT_CPT && nTokens >= 2)
754
1010
        {
755
1011
            pasColorAssociation =
756
1012
                    (ColorAssociation*)CPLRealloc(pasColorAssociation,
1253
1509
    int nXSize = GDALGetRasterBandXSize(hSrcBand);
1254
1510
    int nYSize = GDALGetRasterBandYSize(hSrcBand);
1255
1511
 
1256
 
    FILE* fp = VSIFOpenL(pszDstFilename, "wt");
 
1512
    VSILFILE* fp = VSIFOpenL(pszDstFilename, "wt");
1257
1513
    if (fp == NULL)
1258
1514
    {
1259
1515
        CPLFree(pasColorAssociation);
1282
1538
    GDALGetBlockSize(hSrcBand, &nBlockXSize, &nBlockYSize);
1283
1539
    
1284
1540
    int bRelativeToVRT;
 
1541
    CPLString osPath = CPLGetPath(pszDstFilename);
1285
1542
    char* pszSourceFilename = CPLStrdup(
1286
 
        CPLExtractRelativePath( pszDstFilename, GDALGetDescription(hSrcDataset), 
 
1543
        CPLExtractRelativePath( osPath.c_str(), GDALGetDescription(hSrcDataset), 
1287
1544
                                &bRelativeToVRT ));
1288
1545
 
1289
1546
    for(iBand = 0; iBand < nBands; iBand++)
1290
1547
    {
1291
1548
        VSIFPrintfL(fp, "  <VRTRasterBand dataType=\"Byte\" band=\"%d\">\n", iBand + 1);
 
1549
        VSIFPrintfL(fp, "    <ColorInterp>%s</ColorInterp>\n",
 
1550
                    GDALGetColorInterpretationName((GDALColorInterp)(GCI_RedBand + iBand)));
1292
1551
        VSIFPrintfL(fp, "    <ComplexSource>\n");
1293
1552
        VSIFPrintfL(fp, "      <SourceFilename relativeToVRT=\"%d\">%s</SourceFilename>\n",
1294
1553
                        bRelativeToVRT, pszSourceFilename);
1460
1719
    int                bDstHasNoData;
1461
1720
    double             dfDstNoDataValue;
1462
1721
    int                nCurLine;
 
1722
    int                bComputeAtEdges;
1463
1723
 
1464
1724
  public:
1465
1725
                        GDALGeneric3x3Dataset(GDALDatasetH hSrcDS,
1468
1728
                                              int bDstHasNoData,
1469
1729
                                              double dfDstNoDataValue,
1470
1730
                                              GDALGeneric3x3ProcessingAlg pfnAlg,
1471
 
                                              void* pAlgData);
 
1731
                                              void* pAlgData,
 
1732
                                              int bComputeAtEdges);
1472
1733
                       ~GDALGeneric3x3Dataset();
1473
1734
 
1474
1735
    CPLErr      GetGeoTransform( double * padfGeoTransform );
1484
1745
class GDALGeneric3x3RasterBand : public GDALRasterBand
1485
1746
{
1486
1747
    friend class GDALGeneric3x3Dataset;
1487
 
 
 
1748
    int bSrcHasNoData;
 
1749
    float fSrcNoDataValue;
 
1750
    
 
1751
    void                    InitWidthNoData(void* pImage);
1488
1752
    
1489
1753
  public:
1490
1754
                 GDALGeneric3x3RasterBand( GDALGeneric3x3Dataset *poDS,
1501
1765
                                     int bDstHasNoData,
1502
1766
                                     double dfDstNoDataValue,
1503
1767
                                     GDALGeneric3x3ProcessingAlg pfnAlg,
1504
 
                                     void* pAlgData)
 
1768
                                     void* pAlgData,
 
1769
                                     int bComputeAtEdges)
1505
1770
{
1506
1771
    this->hSrcDS = hSrcDS;
1507
1772
    this->hSrcBand = hSrcBand;
1509
1774
    this->pAlgData = pAlgData;
1510
1775
    this->bDstHasNoData = bDstHasNoData;
1511
1776
    this->dfDstNoDataValue = dfDstNoDataValue;
 
1777
    this->bComputeAtEdges = bComputeAtEdges;
1512
1778
    
1513
1779
    CPLAssert(eDstDataType == GDT_Byte || eDstDataType == GDT_Float32);
1514
1780
 
1520
1786
    apafSourceBuf[0] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1521
1787
    apafSourceBuf[1] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1522
1788
    apafSourceBuf[2] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
1523
 
    
 
1789
 
1524
1790
    nCurLine = -1;
1525
1791
}
1526
1792
 
1549
1815
    eDataType = eDstDataType;
1550
1816
    nBlockXSize = poDS->GetRasterXSize();
1551
1817
    nBlockYSize = 1;
 
1818
 
 
1819
    bSrcHasNoData = FALSE;
 
1820
    fSrcNoDataValue = (float)GDALGetRasterNoDataValue(poDS->hSrcBand,
 
1821
                                                      &bSrcHasNoData);
 
1822
}
 
1823
 
 
1824
void   GDALGeneric3x3RasterBand::InitWidthNoData(void* pImage)
 
1825
{
 
1826
    int j;
 
1827
    GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
 
1828
    if (eDataType == GDT_Byte)
 
1829
    {
 
1830
        for(j=0;j<nBlockXSize;j++)
 
1831
            ((GByte*)pImage)[j] = (GByte) poGDS->dfDstNoDataValue;
 
1832
    }
 
1833
    else
 
1834
    {
 
1835
        for(j=0;j<nBlockXSize;j++)
 
1836
            ((float*)pImage)[j] = (float) poGDS->dfDstNoDataValue;
 
1837
    }
1552
1838
}
1553
1839
 
1554
1840
CPLErr GDALGeneric3x3RasterBand::IReadBlock( int nBlockXOff,
1555
1841
                                             int nBlockYOff,
1556
1842
                                             void *pImage )
1557
1843
{
 
1844
    int i, j;
 
1845
    float fVal;
1558
1846
    GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
1559
 
    
1560
 
    if ( nBlockYOff == 0 || nBlockYOff == nRasterYSize - 1)
1561
 
    {
1562
 
        int j;
1563
 
        if (eDataType == GDT_Byte)
1564
 
        {
1565
 
            for(j=0;j<nBlockXSize;j++)
1566
 
                ((GByte*)pImage)[j] = (GByte) poGDS->dfDstNoDataValue;
1567
 
        }
1568
 
        else
1569
 
        {
1570
 
            for(j=0;j<nBlockXSize;j++)
1571
 
                ((float*)pImage)[j] = (float) poGDS->dfDstNoDataValue;
1572
 
        }
1573
 
            
 
1847
 
 
1848
    if (poGDS->bComputeAtEdges && nRasterXSize >= 2 && nRasterYSize >= 2)
 
1849
    {
 
1850
        if (nBlockYOff == 0)
 
1851
        {
 
1852
            for(i=0;i<2;i++)
 
1853
            {
 
1854
                CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
 
1855
                                    GF_Read,
 
1856
                                    0, i, nBlockXSize, 1,
 
1857
                                    poGDS->apafSourceBuf[i+1],
 
1858
                                    nBlockXSize, 1,
 
1859
                                    GDT_Float32,
 
1860
                                    0, 0);
 
1861
                if (eErr != CE_None)
 
1862
                {
 
1863
                    InitWidthNoData(pImage);
 
1864
                    return eErr;
 
1865
                }
 
1866
            }
 
1867
            poGDS->nCurLine = 0;
 
1868
 
 
1869
            for (j = 0; j < nRasterXSize; j++)
 
1870
            {
 
1871
                float afWin[9];
 
1872
                int jmin = (j == 0) ? j : j - 1;
 
1873
                int jmax = (j == nRasterXSize - 1) ? j : j + 1;
 
1874
 
 
1875
                afWin[0] = INTERPOL(poGDS->apafSourceBuf[1][jmin], poGDS->apafSourceBuf[2][jmin]);
 
1876
                afWin[1] = INTERPOL(poGDS->apafSourceBuf[1][j],    poGDS->apafSourceBuf[2][j]);
 
1877
                afWin[2] = INTERPOL(poGDS->apafSourceBuf[1][jmax], poGDS->apafSourceBuf[2][jmax]);
 
1878
                afWin[3] = poGDS->apafSourceBuf[1][jmin];
 
1879
                afWin[4] = poGDS->apafSourceBuf[1][j];
 
1880
                afWin[5] = poGDS->apafSourceBuf[1][jmax];
 
1881
                afWin[6] = poGDS->apafSourceBuf[2][jmin];
 
1882
                afWin[7] = poGDS->apafSourceBuf[2][j];
 
1883
                afWin[8] = poGDS->apafSourceBuf[2][jmax];
 
1884
 
 
1885
                fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
 
1886
                                    afWin, (float) poGDS->dfDstNoDataValue,
 
1887
                                    poGDS->pfnAlg,
 
1888
                                    poGDS->pAlgData,
 
1889
                                    poGDS->bComputeAtEdges);
 
1890
 
 
1891
                if (eDataType == GDT_Byte)
 
1892
                    ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
 
1893
                else
 
1894
                    ((float*)pImage)[j] = fVal;
 
1895
            }
 
1896
 
 
1897
            return CE_None;
 
1898
        }
 
1899
        else if (nBlockYOff == nRasterYSize - 1)
 
1900
        {
 
1901
            if (poGDS->nCurLine != nRasterYSize - 2)
 
1902
            {
 
1903
                for(i=0;i<2;i++)
 
1904
                {
 
1905
                    CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
 
1906
                                        GF_Read,
 
1907
                                        0, nRasterYSize - 2 + i, nBlockXSize, 1,
 
1908
                                        poGDS->apafSourceBuf[i+1],
 
1909
                                        nBlockXSize, 1,
 
1910
                                        GDT_Float32,
 
1911
                                        0, 0);
 
1912
                    if (eErr != CE_None)
 
1913
                    {
 
1914
                        InitWidthNoData(pImage);
 
1915
                        return eErr;
 
1916
                    }
 
1917
                }
 
1918
            }
 
1919
 
 
1920
            for (j = 0; j < nRasterXSize; j++)
 
1921
            {
 
1922
                float afWin[9];
 
1923
                int jmin = (j == 0) ? j : j - 1;
 
1924
                int jmax = (j == nRasterXSize - 1) ? j : j + 1;
 
1925
 
 
1926
                afWin[0] = poGDS->apafSourceBuf[1][jmin];
 
1927
                afWin[1] = poGDS->apafSourceBuf[1][j];
 
1928
                afWin[2] = poGDS->apafSourceBuf[1][jmax];
 
1929
                afWin[3] = poGDS->apafSourceBuf[2][jmin];
 
1930
                afWin[4] = poGDS->apafSourceBuf[2][j];
 
1931
                afWin[5] = poGDS->apafSourceBuf[2][jmax];
 
1932
                afWin[6] = INTERPOL(poGDS->apafSourceBuf[2][jmin], poGDS->apafSourceBuf[1][jmin]);
 
1933
                afWin[7] = INTERPOL(poGDS->apafSourceBuf[2][j],    poGDS->apafSourceBuf[1][j]);
 
1934
                afWin[8] = INTERPOL(poGDS->apafSourceBuf[2][jmax], poGDS->apafSourceBuf[1][jmax]);
 
1935
 
 
1936
                fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
 
1937
                                    afWin, (float) poGDS->dfDstNoDataValue,
 
1938
                                    poGDS->pfnAlg,
 
1939
                                    poGDS->pAlgData,
 
1940
                                    poGDS->bComputeAtEdges);
 
1941
 
 
1942
                if (eDataType == GDT_Byte)
 
1943
                    ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
 
1944
                else
 
1945
                    ((float*)pImage)[j] = fVal;
 
1946
            }
 
1947
 
 
1948
            return CE_None;
 
1949
        }
 
1950
    }
 
1951
    else if ( nBlockYOff == 0 || nBlockYOff == nRasterYSize - 1)
 
1952
    {
 
1953
        InitWidthNoData(pImage);
1574
1954
        return CE_None;
1575
1955
    }
1576
1956
 
1577
1957
    if ( poGDS->nCurLine != nBlockYOff )
1578
1958
    {
1579
 
        if (nBlockYOff != 1 &&
1580
 
            poGDS->nCurLine == nBlockYOff + 1)
 
1959
        if (poGDS->nCurLine + 1 == nBlockYOff)
1581
1960
        {
1582
1961
            float* pafTmp =  poGDS->apafSourceBuf[0];
1583
1962
            poGDS->apafSourceBuf[0] = poGDS->apafSourceBuf[1];
1584
1963
            poGDS->apafSourceBuf[1] = poGDS->apafSourceBuf[2];
1585
1964
            poGDS->apafSourceBuf[2] = pafTmp;
1586
 
            
1587
 
            poGDS->nCurLine = nBlockYOff;
1588
 
            
 
1965
 
1589
1966
            CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1590
1967
                                    GF_Read,
1591
1968
                                    0, nBlockYOff + 1, nBlockXSize, 1,
1593
1970
                                    nBlockXSize, 1,
1594
1971
                                    GDT_Float32,
1595
1972
                                    0, 0);
1596
 
                                    
 
1973
 
1597
1974
            if (eErr != CE_None)
1598
1975
            {
1599
 
                memset(pImage, 0, nBlockXSize * sizeof(float));
 
1976
                InitWidthNoData(pImage);
1600
1977
                return eErr;
1601
1978
            }
1602
1979
        }
1603
1980
        else
1604
1981
        {
1605
 
            poGDS->nCurLine = nBlockYOff;
1606
 
            int i;
1607
1982
            for(i=0;i<3;i++)
1608
1983
            {
1609
1984
                CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
1615
1990
                                    0, 0);
1616
1991
                if (eErr != CE_None)
1617
1992
                {
1618
 
                    memset(pImage, 0, nBlockXSize * sizeof(float));
 
1993
                    InitWidthNoData(pImage);
1619
1994
                    return eErr;
1620
1995
                }
1621
1996
            }
1622
1997
        }
 
1998
 
 
1999
        poGDS->nCurLine = nBlockYOff;
1623
2000
    }
1624
2001
 
1625
 
    int j;
1626
 
    
1627
 
    if (eDataType == GDT_Byte)
 
2002
    if (poGDS->bComputeAtEdges && nRasterXSize >= 2)
1628
2003
    {
1629
 
        ((GByte*)pImage)[0] = (GByte) poGDS->dfDstNoDataValue;
1630
 
        if (nBlockXSize > 1)
1631
 
            ((GByte*)pImage)[nBlockXSize - 1] = (GByte) poGDS->dfDstNoDataValue;
 
2004
        float afWin[9];
 
2005
 
 
2006
        j = 0;
 
2007
        afWin[0] = INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j+1]);
 
2008
        afWin[1] = poGDS->apafSourceBuf[0][j];
 
2009
        afWin[2] = poGDS->apafSourceBuf[0][j+1];
 
2010
        afWin[3] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j+1]);
 
2011
        afWin[4] = poGDS->apafSourceBuf[1][j];
 
2012
        afWin[5] = poGDS->apafSourceBuf[1][j+1];
 
2013
        afWin[6] = INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j+1]);
 
2014
        afWin[7] = poGDS->apafSourceBuf[2][j];
 
2015
        afWin[8] = poGDS->apafSourceBuf[2][j+1];
 
2016
 
 
2017
        fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
 
2018
                                    afWin, (float) poGDS->dfDstNoDataValue,
 
2019
                                    poGDS->pfnAlg,
 
2020
                                    poGDS->pAlgData,
 
2021
                                    poGDS->bComputeAtEdges);
 
2022
 
 
2023
        if (eDataType == GDT_Byte)
 
2024
            ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
 
2025
        else
 
2026
            ((float*)pImage)[j] = fVal;
 
2027
 
 
2028
        j = nRasterXSize - 1;
 
2029
 
 
2030
        afWin[0] = poGDS->apafSourceBuf[0][j-1];
 
2031
        afWin[1] = poGDS->apafSourceBuf[0][j];
 
2032
        afWin[2] = INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j-1]);
 
2033
        afWin[3] = poGDS->apafSourceBuf[1][j-1];
 
2034
        afWin[4] = poGDS->apafSourceBuf[1][j];
 
2035
        afWin[5] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j-1]);
 
2036
        afWin[6] = poGDS->apafSourceBuf[2][j-1];
 
2037
        afWin[7] = poGDS->apafSourceBuf[2][j];
 
2038
        afWin[8] = INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j-1]);
 
2039
 
 
2040
        fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
 
2041
                                    afWin, (float) poGDS->dfDstNoDataValue,
 
2042
                                    poGDS->pfnAlg,
 
2043
                                    poGDS->pAlgData,
 
2044
                                    poGDS->bComputeAtEdges);
 
2045
 
 
2046
        if (eDataType == GDT_Byte)
 
2047
            ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
 
2048
        else
 
2049
            ((float*)pImage)[j] = fVal;
1632
2050
    }
1633
2051
    else
1634
2052
    {
1635
 
        ((float*)pImage)[0] = (float) poGDS->dfDstNoDataValue;
1636
 
        if (nBlockXSize > 1)
1637
 
            ((float*)pImage)[nBlockXSize - 1] = (float) poGDS->dfDstNoDataValue;
 
2053
        if (eDataType == GDT_Byte)
 
2054
        {
 
2055
            ((GByte*)pImage)[0] = (GByte) poGDS->dfDstNoDataValue;
 
2056
            if (nBlockXSize > 1)
 
2057
                ((GByte*)pImage)[nBlockXSize - 1] = (GByte) poGDS->dfDstNoDataValue;
 
2058
        }
 
2059
        else
 
2060
        {
 
2061
            ((float*)pImage)[0] = (float) poGDS->dfDstNoDataValue;
 
2062
            if (nBlockXSize > 1)
 
2063
                ((float*)pImage)[nBlockXSize - 1] = (float) poGDS->dfDstNoDataValue;
 
2064
        }
1638
2065
    }
1639
 
    
1640
 
    int bSrcHasNoData;
1641
 
    double dfSrcNoDataValue = GDALGetRasterNoDataValue(poGDS->hSrcBand,
1642
 
                                                       &bSrcHasNoData);
 
2066
 
1643
2067
 
1644
2068
    for(j=1;j<nBlockXSize - 1;j++)
1645
2069
    {
1654
2078
        afWin[7] = poGDS->apafSourceBuf[2][j];
1655
2079
        afWin[8] = poGDS->apafSourceBuf[2][j+1];
1656
2080
 
1657
 
        // Check if afWin has null value
1658
 
        int bContainsNull = FALSE;
1659
 
        if (bSrcHasNoData)
1660
 
        {
1661
 
            for ( int n = 0; n <= 8; n++)
1662
 
            {
1663
 
                if(afWin[n] == dfSrcNoDataValue)
1664
 
                {
1665
 
                    bContainsNull = TRUE;
1666
 
                    break;
1667
 
                }
1668
 
            }
1669
 
        }
1670
 
 
1671
 
        if (bContainsNull)
1672
 
        {
1673
 
            // We have nulls so write nullValue and move on
1674
 
            if (eDataType == GDT_Byte)
1675
 
                ((GByte*)pImage)[j] = (GByte) poGDS->dfDstNoDataValue;
1676
 
            else
1677
 
                ((float*)pImage)[j] = (float) poGDS->dfDstNoDataValue;
1678
 
        } else {
1679
 
            // We have a valid 3x3 window.
1680
 
            if (eDataType == GDT_Byte)
1681
 
                ((GByte*)pImage)[j] = (GByte) (poGDS->pfnAlg(
1682
 
                        afWin, poGDS->dfDstNoDataValue, poGDS->pAlgData) + 0.5);
1683
 
            else
1684
 
                ((float*)pImage)[j] = (float) poGDS->pfnAlg(
1685
 
                        afWin, poGDS->dfDstNoDataValue, poGDS->pAlgData);
1686
 
        }
 
2081
        fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
 
2082
                                afWin, (float) poGDS->dfDstNoDataValue,
 
2083
                                poGDS->pfnAlg,
 
2084
                                poGDS->pAlgData,
 
2085
                                poGDS->bComputeAtEdges);
 
2086
 
 
2087
        if (eDataType == GDT_Byte)
 
2088
            ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
 
2089
        else
 
2090
            ((float*)pImage)[j] = fVal;
 
2091
 
1687
2092
    }
1688
 
    
 
2093
 
1689
2094
    return CE_None;
1690
2095
}
1691
2096
 
1715
2120
int main( int argc, char ** argv )
1716
2121
 
1717
2122
{
1718
 
    Algorithm eUtilityMode;
 
2123
    Algorithm eUtilityMode = HILL_SHADE;
1719
2124
    double z = 1.0;
1720
2125
    double scale = 1.0;
1721
2126
    double az = 315.0;
1734
2139
    const char *pszDstFilename = NULL;
1735
2140
    const char *pszColorFilename = NULL;
1736
2141
    const char *pszFormat = "GTiff";
 
2142
    int bFormatExplicitelySet = FALSE;
1737
2143
    char **papszCreateOptions = NULL;
1738
2144
    
1739
2145
    GDALDatasetH hSrcDataset = NULL;
1747
2153
    int nXSize = 0;
1748
2154
    int nYSize = 0;
1749
2155
    
 
2156
    int bComputeAtEdges = FALSE;
 
2157
    int bZevenbergenThorne = FALSE;
 
2158
 
 
2159
    int bQuiet = FALSE;
 
2160
    
1750
2161
    /* Check strict compilation and runtime library version as we use C++ API */
1751
2162
    if (! GDAL_CHECK_VERSION(argv[0]))
1752
2163
        exit(1);
1756
2167
    {
1757
2168
        fprintf(stderr, "Not enough arguments\n");
1758
2169
        Usage();
1759
 
        exit( 1 );
1760
2170
    }
1761
2171
 
1762
2172
    if( EQUAL(argv[1], "--utility_version") || EQUAL(argv[1], "--utility-version") )
1797
2207
    {
1798
2208
        fprintf(stderr, "Missing valid sub-utility mention\n");
1799
2209
        Usage();
1800
 
        exit( 1 );
1801
2210
    }
1802
2211
 
1803
2212
/* -------------------------------------------------------------------- */
1814
2223
        {
1815
2224
            slopeFormat = 0;
1816
2225
        }
 
2226
        else if ( (eUtilityMode == HILL_SHADE || eUtilityMode == SLOPE ||
 
2227
                   eUtilityMode == ASPECT) && EQUAL(argv[i], "-alg") && i + 1 < argc)
 
2228
        {
 
2229
            i ++;
 
2230
            if (EQUAL(argv[i], "ZevenbergenThorne"))
 
2231
                bZevenbergenThorne = TRUE;
 
2232
            else if (!EQUAL(argv[i], "Horn"))
 
2233
            {
 
2234
                fprintf(stderr, "Wrong value for alg : %s\n", argv[i]);
 
2235
                Usage();
 
2236
            }
 
2237
        }
1817
2238
        else if ( eUtilityMode == ASPECT && EQUAL(argv[i], "-trigonometric"))
1818
2239
        {
1819
2240
            bAngleAsAzimuth = FALSE;
1862
2283
        {
1863
2284
            bAddAlpha = TRUE;
1864
2285
        }
 
2286
        else if( eUtilityMode != COLOR_RELIEF &&
 
2287
                 EQUAL(argv[i], "-compute_edges"))
 
2288
        {
 
2289
            bComputeAtEdges = TRUE;
 
2290
        }
1865
2291
        else if( i + 1 < argc &&
1866
2292
            (EQUAL(argv[i], "--b") || 
1867
2293
             EQUAL(argv[i], "-b"))
1872
2298
        else if ( EQUAL(argv[i], "-q") || EQUAL(argv[i], "-quiet") )
1873
2299
        {
1874
2300
            pfnProgress = GDALDummyProgress;
 
2301
            bQuiet = TRUE;
1875
2302
        }
1876
2303
        else if( EQUAL(argv[i],"-co") && i < argc-1 )
1877
2304
        {
1880
2307
        else if( EQUAL(argv[i],"-of") && i < argc-1 )
1881
2308
        {
1882
2309
            pszFormat = argv[++i];
 
2310
            bFormatExplicitelySet = TRUE;
1883
2311
        }
1884
2312
        else if( argv[i][0] == '-' )
1885
2313
        {
1886
2314
            fprintf( stderr, "Option %s incomplete, or not recognised.\n\n", 
1887
2315
                    argv[i] );
1888
2316
            Usage();
1889
 
            GDALDestroyDriverManager();
1890
 
            exit( 2 );
1891
2317
        }
1892
2318
        else if( pszSrcFilename == NULL )
1893
2319
        {
1977
2403
        exit( 1 );
1978
2404
    }
1979
2405
 
 
2406
    if (!bQuiet && !bFormatExplicitelySet)
 
2407
        CheckExtensionConsistency(pszDstFilename, pszFormat);
 
2408
 
1980
2409
    double dfDstNoDataValue = 0;
1981
2410
    int bDstHasNoData = FALSE;
1982
2411
    void* pData = NULL;
1990
2419
                                           z,
1991
2420
                                           scale,
1992
2421
                                           alt,
1993
 
                                           az);
1994
 
        pfnAlg = GDALHillshadeAlg;
 
2422
                                           az,
 
2423
                                           bZevenbergenThorne);
 
2424
        if (bZevenbergenThorne)
 
2425
            pfnAlg = GDALHillshadeZevenbergenThorneAlg;
 
2426
        else
 
2427
            pfnAlg = GDALHillshadeAlg;
1995
2428
    }
1996
2429
    else if (eUtilityMode == SLOPE)
1997
2430
    {
1999
2432
        bDstHasNoData = TRUE;
2000
2433
 
2001
2434
        pData = GDALCreateSlopeData(adfGeoTransform, scale, slopeFormat);
2002
 
        pfnAlg = GDALSlopeAlg;
 
2435
        if (bZevenbergenThorne)
 
2436
            pfnAlg = GDALSlopeZevenbergenThorneAlg;
 
2437
        else
 
2438
            pfnAlg = GDALSlopeHornAlg;
2003
2439
    }
2004
2440
 
2005
2441
    else if (eUtilityMode == ASPECT)
2011
2447
        }
2012
2448
 
2013
2449
        pData = GDALCreateAspectData(bAngleAsAzimuth);
2014
 
        pfnAlg = GDALAspectAlg;
 
2450
        if (bZevenbergenThorne)
 
2451
            pfnAlg = GDALAspectZevenbergenThorneAlg;
 
2452
        else
 
2453
            pfnAlg = GDALAspectAlg;
2015
2454
    }
2016
2455
    else if (eUtilityMode == TRI)
2017
2456
    {
2084
2523
                                          bDstHasNoData,
2085
2524
                                          dfDstNoDataValue,
2086
2525
                                          pfnAlg,
2087
 
                                          pData);
 
2526
                                          pData,
 
2527
                                          bComputeAtEdges);
2088
2528
 
2089
2529
        GDALDatasetH hOutDS = GDALCreateCopy(
2090
2530
                                 hDriver, pszDstFilename, hIntermediateDataset, 
2151
2591
        
2152
2592
        GDALGeneric3x3Processing(hSrcBand, hDstBand,
2153
2593
                                 pfnAlg, pData,
 
2594
                                 bComputeAtEdges,
2154
2595
                                 pfnProgress, NULL);
2155
2596
                                    
2156
2597
    }