Index: /trunk/zoo-project/zoo-services/gdal/dem/Makefile
===================================================================
--- /trunk/zoo-project/zoo-services/gdal/dem/Makefile (revision 412)
+++ /trunk/zoo-project/zoo-services/gdal/dem/Makefile (revision 412)
@@ -0,0 +1,10 @@
+ZRPATH=../../..
+include ${ZRPATH}/zoo-kernel/ZOOMakefile.opts
+CFLAGS=-DZOO_SERVICE ${ZOO_CFLAGS} ${XML2CFLAGS} ${GDAL_CFLAGS} ${PYTHONCFLAGS} -DLINUX_FREE_ISSUE #-DDEBUG
+CC=gcc
+
+cgi-env/dem_service.zo: service.c
+ g++ ${CFLAGS} -shared -fpic -o cgi-env/dem_service.zo ./service.c ${GDAL_LIBS} ${MACOS_LD_FLAGS}
+
+clean:
+ rm -f cgi-env/*.zo
Index: /trunk/zoo-project/zoo-services/gdal/dem/cgi-env/Gdal_Dem.zcfg
===================================================================
--- /trunk/zoo-project/zoo-services/gdal/dem/cgi-env/Gdal_Dem.zcfg (revision 412)
+++ /trunk/zoo-project/zoo-services/gdal/dem/cgi-env/Gdal_Dem.zcfg (revision 412)
@@ -0,0 +1,56 @@
+[Gdal_Dem]
+ Title = Tools to analyze and visualize DEMs.
+ Abstract = http://www.gdal.org/gdaldem.html
+ processVersion = 1
+ storeSupported = true
+ statusSupported = true
+ serviceType = C
+ serviceProvider = dem_service.zo
+
+ [InputDSN]
+ Title = The input data source name
+ Abstract = The input data source name to use.
+ minOccurs = 1
+ maxOccurs = 1
+
+ DataType = string
+
+
+ [OutputDSN]
+ Title = The output data source name
+ Abstract = The output data source name.
+ minOccurs = 1
+ maxOccurs = 1
+
+ DataType = string
+
+
+ [utility]
+ Title = The mode
+ Abstract = This service has 7 different modes : hillshade, slope, aspect, color-relief, TRI, TPI, roughness
+ minOccurs = 1
+ maxOccurs = 1
+
+ DataType = string
+
+
+ [s]
+ Title = The scale
+ Abstract = Scale is the ratio of vertical units to horizontal, for Feet:Latlong use 370400, for Meters:LatLong use 111120
+ minOccurs = 0
+ maxOccurs = 1
+
+ DataType = string
+
+
+
+
+ [Result]
+ Title = The resulting converted file
+ Abstract = The file name resulting of the convertion
+
+ DataType = string
+
+
+
+
Index: /trunk/zoo-project/zoo-services/gdal/dem/makefile.vc
===================================================================
--- /trunk/zoo-project/zoo-services/gdal/dem/makefile.vc (revision 412)
+++ /trunk/zoo-project/zoo-services/gdal/dem/makefile.vc (revision 412)
@@ -0,0 +1,23 @@
+TROOT= $(ZOOK_DIR)
+!INCLUDE $(ZOOK_DIR)/nmake.opt
+
+FILES=$(ZOOK_DIR)/service_internal.obj $(ZOOK_DIR)/ulinet.obj
+!IFDEF JS_DIR
+FILES=$(FILES) $(ZOOK_DIR)$(JS_FILE)
+!ENDIF
+!IFDEF JDK_DIR
+FILES=$(FILES) $(ZOOK_DIR)$(JAVA_FILE)
+!ENDIF
+!IFDEF PY_DIR
+FILES=$(FILES) $(ZOOK_DIR)$(PY_FILE)
+!ENDIF
+!IFDEF MS_DIR
+FILES=$(FILES) $(ZOOK_DIR)$(MS_FILE)
+!ENDIF
+
+cgi-env/dem_service.zo: service.c
+ $(CPP) /D_CRT_NONSTDC_NO_DEPRECATE -DWIN32 -DZOO_SERVICE -I$(GDAL_DIR)/frmts/vrt $(CFLAGS) /c service.c
+ link /dll /out:cgi-env/dem_service.zo $(FILES) ./service.obj $(LDFLAGS) /FORCE:MULTIPLE
+
+clean:
+ erase cgi-env\*.zo
Index: /trunk/zoo-project/zoo-services/gdal/dem/service.c
===================================================================
--- /trunk/zoo-project/zoo-services/gdal/dem/service.c (revision 412)
+++ /trunk/zoo-project/zoo-services/gdal/dem/service.c (revision 412)
@@ -0,0 +1,2899 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project: GDAL DEM Utilities
+ * Purpose:
+ * Authors: Matthew Perry, perrygeo at gmail.com
+ * Even Rouault, even dot rouault at mines dash paris dot org
+ * Howard Butler, hobu.inc at gmail.com
+ * Chris Yesson, chris dot yesson at ioz dot ac dot uk
+ *
+ ******************************************************************************
+ * Copyright (c) 2006, 2009 Matthew Perry
+ * Copyright (c) 2009 Even Rouault
+ * Portions derived from GRASS 4.1 (public domain) See
+ * http://trac.osgeo.org/gdal/ticket/2975 for more information regarding
+ * history of this code
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************
+ *
+ * Slope and aspect calculations based on original method for GRASS GIS 4.1
+ * by Michael Shapiro, U.S.Army Construction Engineering Research Laboratory
+ * Olga Waupotitsch, U.S.Army Construction Engineering Research Laboratory
+ * Marjorie Larson, U.S.Army Construction Engineering Research Laboratory
+ * as found in GRASS's r.slope.aspect module.
+ *
+ * Horn's formula is used to find the first order derivatives in x and y directions
+ * for slope and aspect calculations: Horn, B. K. P. (1981).
+ * "Hill Shading and the Reflectance Map", Proceedings of the IEEE, 69(1):14-47.
+ *
+ * Other reference :
+ * Burrough, P.A. and McDonell, R.A., 1998. Principles of Geographical Information
+ * Systems. p. 190.
+ *
+ * Shaded relief based on original method for GRASS GIS 4.1 by Jim Westervelt,
+ * U.S. Army Construction Engineering Research Laboratory
+ * as found in GRASS's r.shaded.relief (formerly shade.rel.sh) module.
+ * ref: "r.mapcalc: An Algebra for GIS and Image Processing",
+ * by Michael Shapiro and Jim Westervelt, U.S. Army Construction Engineering
+ * Research Laboratory (March/1991)
+ *
+ * Color table of named colors and lookup code derived from src/libes/gis/named_colr.c
+ * of GRASS 4.1
+ *
+ * TRI - Terrain Ruggedness Index is as descibed in Wilson et al (2007)
+ * this is based on the method of Valentine et al. (2004)
+ *
+ * TPI - Topographic Position Index follows the description in Wilson et al (2007), following Weiss (2001)
+ * The radius is fixed at 1 cell width/height
+ *
+ * Roughness - follows the definition in Wilson et al. (2007), which follows Dartnell (2000)
+ *
+ * References for TRI/TPI/Roughness:
+ * Dartnell, P. 2000. Applying Remote Sensing Techniques to map Seafloor
+ * Geology/Habitat Relationships. Masters Thesis, San Francisco State
+ * University, pp. 108.
+ * Valentine, P. C., S. J. Fuller, L. A. Scully. 2004. Terrain Ruggedness
+ * Analysis and Distribution of Boulder Ridges in the Stellwagen Bank National
+ * Marine Sanctuary Region (poster). Galway, Ireland: 5th International
+ * Symposium on Marine Geological and Biological Habitat Mapping (GeoHAB),
+ * May 2004.
+ * Weiss, A. D. 2001. Topographic Positions and Landforms Analysis (poster),
+ * ESRI International User Conference, July 2001. San Diego, CA: ESRI.
+ * Wilson, M. F. J.; O'Connell, B.; Brown, C.; Guinan, J. C. & Grehan, A. J.
+ * Multiscale terrain analysis of multibeam bathymetry data for habitat mapping
+ * on the continental slope Marine Geodesy, 2007, 30, 3-35
+ ****************************************************************************/
+
+#include
+#include
+
+#include "cpl_conv.h"
+#include "cpl_string.h"
+#include "gdal.h"
+#include "gdal_priv.h"
+#ifndef ZOO_SERVICE
+#include "commonutils.h"
+#else
+#include "service.h"
+#endif
+
+CPL_CVSID("$Id$");
+
+#ifdef WIN32
+#define strcasecmp _stricmp
+#define strncasecmp _strnicmp
+#endif
+
+#ifdef ZOO_SERVICE
+extern "C" {
+#endif
+
+#ifndef M_PI
+# define M_PI 3.1415926535897932384626433832795
+#endif
+
+#define INTERPOL(a,b) ((bSrcHasNoData && (ARE_REAL_EQUAL(a, fSrcNoDataValue) || ARE_REAL_EQUAL(b, fSrcNoDataValue))) ? fSrcNoDataValue : 2 * (a) - (b))
+
+#ifdef ZOO_SERVICE
+static int Usage(maps*& conf, const char* pszErrorMsg = NULL){
+ setMapInMaps(conf,"lenv","message","Missing source.");
+ return SERVICE_FAILED;
+}
+#else
+/************************************************************************/
+/* Usage() */
+/************************************************************************/
+
+static void Usage(const char* pszErrorMsg = NULL)
+
+{
+ printf( " Usage: \n"
+ " - To generate a shaded relief map from any GDAL-supported elevation raster : \n\n"
+ " gdaldem hillshade input_dem output_hillshade \n"
+ " [-z ZFactor (default=1)] [-s scale* (default=1)] \n"
+ " [-az Azimuth (default=315)] [-alt Altitude (default=45)]\n"
+ " [-alg ZevenbergenThorne] [-combined]\n"
+ " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
+ "\n"
+ " - To generates a slope map from any GDAL-supported elevation raster :\n\n"
+ " gdaldem slope input_dem output_slope_map \n"
+ " [-p use percent slope (default=degrees)] [-s scale* (default=1)]\n"
+ " [-alg ZevenbergenThorne]\n"
+ " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
+ "\n"
+ " - To generate an aspect map from any GDAL-supported elevation raster\n"
+ " Outputs a 32-bit float tiff with pixel values from 0-360 indicating azimuth :\n\n"
+ " gdaldem aspect input_dem output_aspect_map \n"
+ " [-trigonometric] [-zero_for_flat]\n"
+ " [-alg ZevenbergenThorne]\n"
+ " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
+ "\n"
+ " - To generate a color relief map from any GDAL-supported elevation raster\n"
+ " gdaldem color-relief input_dem color_text_file output_color_relief_map\n"
+ " [-alpha] [-exact_color_entry | -nearest_color_entry]\n"
+ " [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
+ " where color_text_file contains lines of the format \"elevation_value red green blue\"\n"
+ "\n"
+ " - To generate a Terrain Ruggedness Index (TRI) map from any GDAL-supported elevation raster\n"
+ " gdaldem TRI input_dem output_TRI_map\n"
+ " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
+ "\n"
+ " - To generate a Topographic Position Index (TPI) map from any GDAL-supported elevation raster\n"
+ " gdaldem TPI input_dem output_TPI_map\n"
+ " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
+ "\n"
+ " - To generate a roughness map from any GDAL-supported elevation raster\n"
+ " gdaldem roughness input_dem output_roughness_map\n"
+ " [-compute_edges] [-b Band (default=1)] [-of format] [-co \"NAME=VALUE\"]* [-q]\n"
+ "\n"
+ " Notes : \n"
+ " Scale is the ratio of vertical units to horizontal\n"
+ " for Feet:Latlong use scale=370400, for Meters:LatLong use scale=111120 \n\n");
+
+ if( pszErrorMsg != NULL )
+ fprintf(stderr, "\nFAILURE: %s\n", pszErrorMsg);
+
+ exit( 1 );
+}
+#endif
+/************************************************************************/
+/* ComputeVal() */
+/************************************************************************/
+
+typedef float (*GDALGeneric3x3ProcessingAlg) (float* pafWindow, float fDstNoDataValue, void* pData);
+
+static float ComputeVal(int bSrcHasNoData, float fSrcNoDataValue,
+ float* afWin, float fDstNoDataValue,
+ GDALGeneric3x3ProcessingAlg pfnAlg,
+ void* pData,
+ int bComputeAtEdges)
+{
+ if (bSrcHasNoData && ARE_REAL_EQUAL(afWin[4], fSrcNoDataValue))
+ {
+ return fDstNoDataValue;
+ }
+ else if (bSrcHasNoData)
+ {
+ int k;
+ for(k=0;k<9;k++)
+ {
+ if (ARE_REAL_EQUAL(afWin[k], fSrcNoDataValue))
+ {
+ if (bComputeAtEdges)
+ afWin[k] = afWin[4];
+ else
+ return fDstNoDataValue;
+ }
+ }
+ }
+
+ return pfnAlg(afWin, fDstNoDataValue, pData);
+}
+
+/************************************************************************/
+/* GDALGeneric3x3Processing() */
+/************************************************************************/
+
+CPLErr GDALGeneric3x3Processing ( GDALRasterBandH hSrcBand,
+ GDALRasterBandH hDstBand,
+ GDALGeneric3x3ProcessingAlg pfnAlg,
+ void* pData,
+ int bComputeAtEdges,
+ GDALProgressFunc pfnProgress,
+ void * pProgressData)
+{
+ CPLErr eErr;
+ float *pafThreeLineWin; /* 3 line rotating source buffer */
+ float *pafOutputBuf; /* 1 line destination buffer */
+ int i, j;
+
+ int bSrcHasNoData, bDstHasNoData;
+ float fSrcNoDataValue = 0.0, fDstNoDataValue = 0.0;
+
+ int nXSize = GDALGetRasterBandXSize(hSrcBand);
+ int nYSize = GDALGetRasterBandYSize(hSrcBand);
+
+ if (pfnProgress == NULL)
+ pfnProgress = GDALDummyProgress;
+
+/* -------------------------------------------------------------------- */
+/* Initialize progress counter. */
+/* -------------------------------------------------------------------- */
+ if( !pfnProgress( 0.0, NULL, pProgressData ) )
+ {
+ CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+ return CE_Failure;
+ }
+
+ pafOutputBuf = (float *) CPLMalloc(sizeof(float)*nXSize);
+ pafThreeLineWin = (float *) CPLMalloc(3*sizeof(float)*(nXSize+1));
+
+ fSrcNoDataValue = (float) GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
+ fDstNoDataValue = (float) GDALGetRasterNoDataValue(hDstBand, &bDstHasNoData);
+ if (!bDstHasNoData)
+ fDstNoDataValue = 0.0;
+
+ // Move a 3x3 pafWindow over each cell
+ // (where the cell in question is #4)
+ //
+ // 0 1 2
+ // 3 4 5
+ // 6 7 8
+
+ /* Preload the first 2 lines */
+ for ( i = 0; i < 2 && i < nYSize; i++)
+ {
+ GDALRasterIO( hSrcBand,
+ GF_Read,
+ 0, i,
+ nXSize, 1,
+ pafThreeLineWin + i * nXSize,
+ nXSize, 1,
+ GDT_Float32,
+ 0, 0);
+ }
+
+ if (bComputeAtEdges && nXSize >= 2 && nYSize >= 2)
+ {
+ for (j = 0; j < nXSize; j++)
+ {
+ float afWin[9];
+ int jmin = (j == 0) ? j : j - 1;
+ int jmax = (j == nXSize - 1) ? j : j + 1;
+
+ afWin[0] = INTERPOL(pafThreeLineWin[jmin], pafThreeLineWin[nXSize + jmin]);
+ afWin[1] = INTERPOL(pafThreeLineWin[j], pafThreeLineWin[nXSize + j]);
+ afWin[2] = INTERPOL(pafThreeLineWin[jmax], pafThreeLineWin[nXSize + jmax]);
+ afWin[3] = pafThreeLineWin[jmin];
+ afWin[4] = pafThreeLineWin[j];
+ afWin[5] = pafThreeLineWin[jmax];
+ afWin[6] = pafThreeLineWin[nXSize + jmin];
+ afWin[7] = pafThreeLineWin[nXSize + j];
+ afWin[8] = pafThreeLineWin[nXSize + jmax];
+
+ pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
+ afWin, fDstNoDataValue,
+ pfnAlg, pData, bComputeAtEdges);
+ }
+ GDALRasterIO(hDstBand, GF_Write,
+ 0, 0, nXSize, 1,
+ pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
+ }
+ else
+ {
+ // Exclude the edges
+ for (j = 0; j < nXSize; j++)
+ {
+ pafOutputBuf[j] = fDstNoDataValue;
+ }
+ GDALRasterIO(hDstBand, GF_Write,
+ 0, 0, nXSize, 1,
+ pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
+
+ if (nYSize > 1)
+ {
+ GDALRasterIO(hDstBand, GF_Write,
+ 0, nYSize - 1, nXSize, 1,
+ pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
+ }
+ }
+
+ int nLine1Off = 0*nXSize;
+ int nLine2Off = 1*nXSize;
+ int nLine3Off = 2*nXSize;
+
+ for ( i = 1; i < nYSize-1; i++)
+ {
+ /* Read third line of the line buffer */
+ eErr = GDALRasterIO( hSrcBand,
+ GF_Read,
+ 0, i+1,
+ nXSize, 1,
+ pafThreeLineWin + nLine3Off,
+ nXSize, 1,
+ GDT_Float32,
+ 0, 0);
+ if (eErr != CE_None)
+ goto end;
+
+ if (bComputeAtEdges && nXSize >= 2)
+ {
+ float afWin[9];
+
+ j = 0;
+ afWin[0] = INTERPOL(pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j+1]);
+ afWin[1] = pafThreeLineWin[nLine1Off + j];
+ afWin[2] = pafThreeLineWin[nLine1Off + j+1];
+ afWin[3] = INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j+1]);
+ afWin[4] = pafThreeLineWin[nLine2Off + j];
+ afWin[5] = pafThreeLineWin[nLine2Off + j+1];
+ afWin[6] = INTERPOL(pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j+1]);
+ afWin[7] = pafThreeLineWin[nLine3Off + j];
+ afWin[8] = pafThreeLineWin[nLine3Off + j+1];
+
+ pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
+ afWin, fDstNoDataValue,
+ pfnAlg, pData, bComputeAtEdges);
+ j = nXSize - 1;
+
+ afWin[0] = pafThreeLineWin[nLine1Off + j-1];
+ afWin[1] = pafThreeLineWin[nLine1Off + j];
+ afWin[2] = INTERPOL(pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j-1]);
+ afWin[3] = pafThreeLineWin[nLine2Off + j-1];
+ afWin[4] = pafThreeLineWin[nLine2Off + j];
+ afWin[5] = INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j-1]);
+ afWin[6] = pafThreeLineWin[nLine3Off + j-1];
+ afWin[7] = pafThreeLineWin[nLine3Off + j];
+ afWin[8] = INTERPOL(pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j-1]);
+
+ pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
+ afWin, fDstNoDataValue,
+ pfnAlg, pData, bComputeAtEdges);
+ }
+ else
+ {
+ // Exclude the edges
+ pafOutputBuf[0] = fDstNoDataValue;
+ if (nXSize > 1)
+ pafOutputBuf[nXSize - 1] = fDstNoDataValue;
+ }
+
+ for (j = 1; j < nXSize - 1; j++)
+ {
+ float afWin[9];
+ afWin[0] = pafThreeLineWin[nLine1Off + j-1];
+ afWin[1] = pafThreeLineWin[nLine1Off + j];
+ afWin[2] = pafThreeLineWin[nLine1Off + j+1];
+ afWin[3] = pafThreeLineWin[nLine2Off + j-1];
+ afWin[4] = pafThreeLineWin[nLine2Off + j];
+ afWin[5] = pafThreeLineWin[nLine2Off + j+1];
+ afWin[6] = pafThreeLineWin[nLine3Off + j-1];
+ afWin[7] = pafThreeLineWin[nLine3Off + j];
+ afWin[8] = pafThreeLineWin[nLine3Off + j+1];
+
+ pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
+ afWin, fDstNoDataValue,
+ pfnAlg, pData, bComputeAtEdges);
+ }
+
+ /* -----------------------------------------
+ * Write Line to Raster
+ */
+ eErr = GDALRasterIO(hDstBand, GF_Write, 0, i, nXSize, 1,
+ pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
+ if (eErr != CE_None)
+ goto end;
+
+ if( !pfnProgress( 1.0 * (i+1) / nYSize, NULL, pProgressData ) )
+ {
+ CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+ eErr = CE_Failure;
+ goto end;
+ }
+
+ int nTemp = nLine1Off;
+ nLine1Off = nLine2Off;
+ nLine2Off = nLine3Off;
+ nLine3Off = nTemp;
+ }
+
+ if (bComputeAtEdges && nXSize >= 2 && nYSize >= 2)
+ {
+ for (j = 0; j < nXSize; j++)
+ {
+ float afWin[9];
+ int jmin = (j == 0) ? j : j - 1;
+ int jmax = (j == nXSize - 1) ? j : j + 1;
+
+ afWin[0] = pafThreeLineWin[nLine1Off + jmin];
+ afWin[1] = pafThreeLineWin[nLine1Off + j];
+ afWin[2] = pafThreeLineWin[nLine1Off + jmax];
+ afWin[3] = pafThreeLineWin[nLine2Off + jmin];
+ afWin[4] = pafThreeLineWin[nLine2Off + j];
+ afWin[5] = pafThreeLineWin[nLine2Off + jmax];
+ afWin[6] = INTERPOL(pafThreeLineWin[nLine2Off + jmin], pafThreeLineWin[nLine1Off + jmin]);
+ afWin[7] = INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine1Off + j]);
+ afWin[8] = INTERPOL(pafThreeLineWin[nLine2Off + jmax], pafThreeLineWin[nLine1Off + jmax]);
+
+ pafOutputBuf[j] = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
+ afWin, fDstNoDataValue,
+ pfnAlg, pData, bComputeAtEdges);
+ }
+ GDALRasterIO(hDstBand, GF_Write,
+ 0, i, nXSize, 1,
+ pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
+ }
+
+ pfnProgress( 1.0, NULL, pProgressData );
+ eErr = CE_None;
+
+end:
+ CPLFree(pafOutputBuf);
+ CPLFree(pafThreeLineWin);
+
+ return eErr;
+}
+
+
+/************************************************************************/
+/* GDALHillshade() */
+/************************************************************************/
+
+typedef struct
+{
+ double nsres;
+ double ewres;
+ double sin_altRadians;
+ double cos_altRadians_mul_z_scale_factor;
+ double azRadians;
+ double square_z_scale_factor;
+ double square_M_PI_2;
+} GDALHillshadeAlgData;
+
+/* Unoptimized formulas are :
+ x = psData->z*((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
+ (afWin[2] + afWin[5] + afWin[5] + afWin[8])) /
+ (8.0 * psData->ewres * psData->scale);
+
+ y = psData->z*((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
+ (afWin[0] + afWin[1] + afWin[1] + afWin[2])) /
+ (8.0 * psData->nsres * psData->scale);
+
+ slope = M_PI / 2 - atan(sqrt(x*x + y*y));
+
+ aspect = atan2(y,x);
+
+ cang = sin(alt * degreesToRadians) * sin(slope) +
+ cos(alt * degreesToRadians) * cos(slope) *
+ cos(az * degreesToRadians - M_PI/2 - aspect);
+*/
+
+float GDALHillshadeAlg (float* afWin, float fDstNoDataValue, void* pData)
+{
+ GDALHillshadeAlgData* psData = (GDALHillshadeAlgData*)pData;
+ double x, y, aspect, xx_plus_yy, cang;
+
+ // First Slope ...
+ x = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
+ (afWin[2] + afWin[5] + afWin[5] + afWin[8])) / psData->ewres;
+
+ y = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
+ (afWin[0] + afWin[1] + afWin[1] + afWin[2])) / psData->nsres;
+
+ xx_plus_yy = x * x + y * y;
+
+ // ... then aspect...
+ aspect = atan2(y,x);
+
+ // ... then the shade value
+ cang = (psData->sin_altRadians -
+ psData->cos_altRadians_mul_z_scale_factor * sqrt(xx_plus_yy) *
+ sin(aspect - psData->azRadians)) /
+ sqrt(1 + psData->square_z_scale_factor * xx_plus_yy);
+
+ if (cang <= 0.0)
+ cang = 1.0;
+ else
+ cang = 1.0 + (254.0 * cang);
+
+ return (float) cang;
+}
+
+float GDALHillshadeCombinedAlg (float* afWin, float fDstNoDataValue, void* pData)
+{
+ GDALHillshadeAlgData* psData = (GDALHillshadeAlgData*)pData;
+ double x, y, aspect, xx_plus_yy, cang;
+
+ // First Slope ...
+ x = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
+ (afWin[2] + afWin[5] + afWin[5] + afWin[8])) / psData->ewres;
+
+ y = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
+ (afWin[0] + afWin[1] + afWin[1] + afWin[2])) / psData->nsres;
+
+ xx_plus_yy = x * x + y * y;
+
+ // ... then aspect...
+ aspect = atan2(y,x);
+ double slope = xx_plus_yy * psData->square_z_scale_factor;
+
+ // ... then the shade value
+ cang = acos((psData->sin_altRadians -
+ psData->cos_altRadians_mul_z_scale_factor * sqrt(xx_plus_yy) *
+ sin(aspect - psData->azRadians)) /
+ sqrt(1 + slope));
+
+ // combined shading
+ cang = 1 - cang * atan(sqrt(slope)) / psData->square_M_PI_2;
+
+ if (cang <= 0.0)
+ cang = 1.0;
+ else
+ cang = 1.0 + (254.0 * cang);
+
+ return (float) cang;
+}
+
+float GDALHillshadeZevenbergenThorneAlg (float* afWin, float fDstNoDataValue, void* pData)
+{
+ GDALHillshadeAlgData* psData = (GDALHillshadeAlgData*)pData;
+ double x, y, aspect, xx_plus_yy, cang;
+
+ // First Slope ...
+ x = (afWin[3] - afWin[5]) / psData->ewres;
+
+ y = (afWin[7] - afWin[1]) / psData->nsres;
+
+ xx_plus_yy = x * x + y * y;
+
+ // ... then aspect...
+ aspect = atan2(y,x);
+
+ // ... then the shade value
+ cang = (psData->sin_altRadians -
+ psData->cos_altRadians_mul_z_scale_factor * sqrt(xx_plus_yy) *
+ sin(aspect - psData->azRadians)) /
+ sqrt(1 + psData->square_z_scale_factor * xx_plus_yy);
+
+ if (cang <= 0.0)
+ cang = 1.0;
+ else
+ cang = 1.0 + (254.0 * cang);
+
+ return (float) cang;
+}
+
+float GDALHillshadeZevenbergenThorneCombinedAlg (float* afWin, float fDstNoDataValue, void* pData)
+{
+ GDALHillshadeAlgData* psData = (GDALHillshadeAlgData*)pData;
+ double x, y, aspect, xx_plus_yy, cang;
+
+ // First Slope ...
+ x = (afWin[3] - afWin[5]) / psData->ewres;
+
+ y = (afWin[7] - afWin[1]) / psData->nsres;
+
+ xx_plus_yy = x * x + y * y;
+
+ // ... then aspect...
+ aspect = atan2(y,x);
+ double slope = xx_plus_yy * psData->square_z_scale_factor;
+
+ // ... then the shade value
+ cang = acos((psData->sin_altRadians -
+ psData->cos_altRadians_mul_z_scale_factor * sqrt(xx_plus_yy) *
+ sin(aspect - psData->azRadians)) /
+ sqrt(1 + slope));
+
+ // combined shading
+ cang = 1 - cang * atan(sqrt(slope)) / psData->square_M_PI_2;
+
+ if (cang <= 0.0)
+ cang = 1.0;
+ else
+ cang = 1.0 + (254.0 * cang);
+
+ return (float) cang;
+}
+
+void* GDALCreateHillshadeData(double* adfGeoTransform,
+ double z,
+ double scale,
+ double alt,
+ double az,
+ int bZevenbergenThorne)
+{
+ GDALHillshadeAlgData* pData =
+ (GDALHillshadeAlgData*)CPLMalloc(sizeof(GDALHillshadeAlgData));
+
+ const double degreesToRadians = M_PI / 180.0;
+ pData->nsres = adfGeoTransform[5];
+ pData->ewres = adfGeoTransform[1];
+ pData->sin_altRadians = sin(alt * degreesToRadians);
+ pData->azRadians = az * degreesToRadians;
+ double z_scale_factor = z / (((bZevenbergenThorne) ? 2 : 8) * scale);
+ pData->cos_altRadians_mul_z_scale_factor =
+ cos(alt * degreesToRadians) * z_scale_factor;
+ pData->square_z_scale_factor = z_scale_factor * z_scale_factor;
+ pData->square_M_PI_2 = (M_PI*M_PI)/4;
+ return pData;
+}
+
+/************************************************************************/
+/* GDALSlope() */
+/************************************************************************/
+
+typedef struct
+{
+ double nsres;
+ double ewres;
+ double scale;
+ int slopeFormat;
+} GDALSlopeAlgData;
+
+float GDALSlopeHornAlg (float* afWin, float fDstNoDataValue, void* pData)
+{
+ const double radiansToDegrees = 180.0 / M_PI;
+ GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
+ double dx, dy, key;
+
+ dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
+ (afWin[2] + afWin[5] + afWin[5] + afWin[8]))/psData->ewres;
+
+ dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
+ (afWin[0] + afWin[1] + afWin[1] + afWin[2]))/psData->nsres;
+
+ key = (dx * dx + dy * dy);
+
+ if (psData->slopeFormat == 1)
+ return (float) (atan(sqrt(key) / (8*psData->scale)) * radiansToDegrees);
+ else
+ return (float) (100*(sqrt(key) / (8*psData->scale)));
+}
+
+float GDALSlopeZevenbergenThorneAlg (float* afWin, float fDstNoDataValue, void* pData)
+{
+ const double radiansToDegrees = 180.0 / M_PI;
+ GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
+ double dx, dy, key;
+
+ dx = (afWin[3] - afWin[5])/psData->ewres;
+
+ dy = (afWin[7] - afWin[1])/psData->nsres;
+
+ key = (dx * dx + dy * dy);
+
+ if (psData->slopeFormat == 1)
+ return (float) (atan(sqrt(key) / (2*psData->scale)) * radiansToDegrees);
+ else
+ return (float) (100*(sqrt(key) / (2*psData->scale)));
+}
+
+void* GDALCreateSlopeData(double* adfGeoTransform,
+ double scale,
+ int slopeFormat)
+{
+ GDALSlopeAlgData* pData =
+ (GDALSlopeAlgData*)CPLMalloc(sizeof(GDALSlopeAlgData));
+
+ pData->nsres = adfGeoTransform[5];
+ pData->ewres = adfGeoTransform[1];
+ pData->scale = scale;
+ pData->slopeFormat = slopeFormat;
+ return pData;
+}
+
+/************************************************************************/
+/* GDALAspect() */
+/************************************************************************/
+
+typedef struct
+{
+ int bAngleAsAzimuth;
+} GDALAspectAlgData;
+
+float GDALAspectAlg (float* afWin, float fDstNoDataValue, void* pData)
+{
+ const double degreesToRadians = M_PI / 180.0;
+ GDALAspectAlgData* psData = (GDALAspectAlgData*)pData;
+ double dx, dy;
+ float aspect;
+
+ dx = ((afWin[2] + afWin[5] + afWin[5] + afWin[8]) -
+ (afWin[0] + afWin[3] + afWin[3] + afWin[6]));
+
+ dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
+ (afWin[0] + afWin[1] + afWin[1] + afWin[2]));
+
+ aspect = (float) (atan2(dy,-dx) / degreesToRadians);
+
+ if (dx == 0 && dy == 0)
+ {
+ /* Flat area */
+ aspect = fDstNoDataValue;
+ }
+ else if ( psData->bAngleAsAzimuth )
+ {
+ if (aspect > 90.0)
+ aspect = 450.0f - aspect;
+ else
+ aspect = 90.0f - aspect;
+ }
+ else
+ {
+ if (aspect < 0)
+ aspect += 360.0;
+ }
+
+ if (aspect == 360.0)
+ aspect = 0.0;
+
+ return aspect;
+}
+
+float GDALAspectZevenbergenThorneAlg (float* afWin, float fDstNoDataValue, void* pData)
+{
+ const double degreesToRadians = M_PI / 180.0;
+ GDALAspectAlgData* psData = (GDALAspectAlgData*)pData;
+ double dx, dy;
+ float aspect;
+
+ dx = (afWin[5] - afWin[3]);
+
+ dy = (afWin[7] - afWin[1]);
+
+ aspect = (float) (atan2(dy,-dx) / degreesToRadians);
+
+ if (dx == 0 && dy == 0)
+ {
+ /* Flat area */
+ aspect = fDstNoDataValue;
+ }
+ else if ( psData->bAngleAsAzimuth )
+ {
+ if (aspect > 90.0)
+ aspect = 450.0f - aspect;
+ else
+ aspect = 90.0f - aspect;
+ }
+ else
+ {
+ if (aspect < 0)
+ aspect += 360.0;
+ }
+
+ if (aspect == 360.0)
+ aspect = 0.0;
+
+ return aspect;
+}
+void* GDALCreateAspectData(int bAngleAsAzimuth)
+{
+ GDALAspectAlgData* pData =
+ (GDALAspectAlgData*)CPLMalloc(sizeof(GDALAspectAlgData));
+
+ pData->bAngleAsAzimuth = bAngleAsAzimuth;
+ return pData;
+}
+
+/************************************************************************/
+/* GDALColorRelief() */
+/************************************************************************/
+
+typedef struct
+{
+ double dfVal;
+ int nR;
+ int nG;
+ int nB;
+ int nA;
+} ColorAssociation;
+
+static int GDALColorReliefSortColors(const void* pA, const void* pB)
+{
+ ColorAssociation* pC1 = (ColorAssociation*)pA;
+ ColorAssociation* pC2 = (ColorAssociation*)pB;
+ return (pC1->dfVal < pC2->dfVal) ? -1 :
+ (pC1->dfVal == pC2->dfVal) ? 0 : 1;
+}
+
+typedef enum
+{
+ COLOR_SELECTION_INTERPOLATE,
+ COLOR_SELECTION_NEAREST_ENTRY,
+ COLOR_SELECTION_EXACT_ENTRY
+} ColorSelectionMode;
+
+static int GDALColorReliefGetRGBA (ColorAssociation* pasColorAssociation,
+ int nColorAssociation,
+ double dfVal,
+ ColorSelectionMode eColorSelectionMode,
+ int* pnR,
+ int* pnG,
+ int* pnB,
+ int* pnA)
+{
+ int i;
+ int lower = 0;
+ int upper = nColorAssociation - 1;
+ int mid;
+
+ /* Find the index of the first element in the LUT input array that */
+ /* is not smaller than the dfVal value. */
+ while(TRUE)
+ {
+ mid = (lower + upper) / 2;
+ if (upper - lower <= 1)
+ {
+ if (dfVal < pasColorAssociation[lower].dfVal)
+ i = lower;
+ else if (dfVal < pasColorAssociation[upper].dfVal)
+ i = upper;
+ else
+ i = upper + 1;
+ break;
+ }
+ else if (pasColorAssociation[mid].dfVal >= dfVal)
+ {
+ upper = mid;
+ }
+ else
+ {
+ lower = mid;
+ }
+ }
+
+ if (i == 0)
+ {
+ if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
+ pasColorAssociation[0].dfVal != dfVal)
+ {
+ *pnR = 0;
+ *pnG = 0;
+ *pnB = 0;
+ *pnA = 0;
+ return FALSE;
+ }
+ else
+ {
+ *pnR = pasColorAssociation[0].nR;
+ *pnG = pasColorAssociation[0].nG;
+ *pnB = pasColorAssociation[0].nB;
+ *pnA = pasColorAssociation[0].nA;
+ return TRUE;
+ }
+ }
+ else if (i == nColorAssociation)
+ {
+ if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
+ pasColorAssociation[i-1].dfVal != dfVal)
+ {
+ *pnR = 0;
+ *pnG = 0;
+ *pnB = 0;
+ *pnA = 0;
+ return FALSE;
+ }
+ else
+ {
+ *pnR = pasColorAssociation[i-1].nR;
+ *pnG = pasColorAssociation[i-1].nG;
+ *pnB = pasColorAssociation[i-1].nB;
+ *pnA = pasColorAssociation[i-1].nA;
+ return TRUE;
+ }
+ }
+ else
+ {
+ if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
+ pasColorAssociation[i-1].dfVal != dfVal)
+ {
+ *pnR = 0;
+ *pnG = 0;
+ *pnB = 0;
+ *pnA = 0;
+ return FALSE;
+ }
+
+ if (eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY &&
+ pasColorAssociation[i-1].dfVal != dfVal)
+ {
+ int index;
+ if (dfVal - pasColorAssociation[i-1].dfVal <
+ pasColorAssociation[i].dfVal - dfVal)
+ index = i -1;
+ else
+ index = i;
+
+ *pnR = pasColorAssociation[index].nR;
+ *pnG = pasColorAssociation[index].nG;
+ *pnB = pasColorAssociation[index].nB;
+ *pnA = pasColorAssociation[index].nA;
+ return TRUE;
+ }
+
+ double dfRatio = (dfVal - pasColorAssociation[i-1].dfVal) /
+ (pasColorAssociation[i].dfVal - pasColorAssociation[i-1].dfVal);
+ *pnR = (int)(0.45 + pasColorAssociation[i-1].nR + dfRatio *
+ (pasColorAssociation[i].nR - pasColorAssociation[i-1].nR));
+ if (*pnR < 0) *pnR = 0;
+ else if (*pnR > 255) *pnR = 255;
+ *pnG = (int)(0.45 + pasColorAssociation[i-1].nG + dfRatio *
+ (pasColorAssociation[i].nG - pasColorAssociation[i-1].nG));
+ if (*pnG < 0) *pnG = 0;
+ else if (*pnG > 255) *pnG = 255;
+ *pnB = (int)(0.45 + pasColorAssociation[i-1].nB + dfRatio *
+ (pasColorAssociation[i].nB - pasColorAssociation[i-1].nB));
+ if (*pnB < 0) *pnB = 0;
+ else if (*pnB > 255) *pnB = 255;
+ *pnA = (int)(0.45 + pasColorAssociation[i-1].nA + dfRatio *
+ (pasColorAssociation[i].nA - pasColorAssociation[i-1].nA));
+ if (*pnA < 0) *pnA = 0;
+ else if (*pnA > 255) *pnA = 255;
+
+ return TRUE;
+ }
+}
+
+/* dfPct : percentage between 0 and 1 */
+static double GDALColorReliefGetAbsoluteValFromPct(GDALRasterBandH hSrcBand,
+ double dfPct)
+{
+ double dfMin, dfMax;
+ int bSuccessMin, bSuccessMax;
+ dfMin = GDALGetRasterMinimum(hSrcBand, &bSuccessMin);
+ dfMax = GDALGetRasterMaximum(hSrcBand, &bSuccessMax);
+ if (!bSuccessMin || !bSuccessMax)
+ {
+ double dfMean, dfStdDev;
+ fprintf(stderr, "Computing source raster statistics...\n");
+ GDALComputeRasterStatistics(hSrcBand, FALSE, &dfMin, &dfMax,
+ &dfMean, &dfStdDev, NULL, NULL);
+ }
+ return dfMin + dfPct * (dfMax - dfMin);
+}
+
+typedef struct
+{
+ const char *name;
+ float r, g, b;
+} NamedColor;
+
+static const NamedColor namedColors[] = {
+ { "white", 1.00, 1.00, 1.00 },
+ { "black", 0.00, 0.00, 0.00 },
+ { "red", 1.00, 0.00, 0.00 },
+ { "green", 0.00, 1.00, 0.00 },
+ { "blue", 0.00, 0.00, 1.00 },
+ { "yellow", 1.00, 1.00, 0.00 },
+ { "magenta",1.00, 0.00, 1.00 },
+ { "cyan", 0.00, 1.00, 1.00 },
+ { "aqua", 0.00, 0.75, 0.75 },
+ { "grey", 0.75, 0.75, 0.75 },
+ { "gray", 0.75, 0.75, 0.75 },
+ { "orange", 1.00, 0.50, 0.00 },
+ { "brown", 0.75, 0.50, 0.25 },
+ { "purple", 0.50, 0.00, 1.00 },
+ { "violet", 0.50, 0.00, 1.00 },
+ { "indigo", 0.00, 0.50, 1.00 },
+};
+
+static
+int GDALColorReliefFindNamedColor(const char *pszColorName, int *pnR, int *pnG, int *pnB)
+{
+ unsigned int i;
+
+ *pnR = *pnG = *pnB = 0;
+ for (i = 0; i < sizeof(namedColors) / sizeof(namedColors[0]); i++)
+ {
+ if (EQUAL(pszColorName, namedColors[i].name))
+ {
+ *pnR = (int)(255. * namedColors[i].r);
+ *pnG = (int)(255. * namedColors[i].g);
+ *pnB = (int)(255. * namedColors[i].b);
+ return TRUE;
+ }
+ }
+ return FALSE;
+}
+
+static
+ColorAssociation* GDALColorReliefParseColorFile(GDALRasterBandH hSrcBand,
+ const char* pszColorFilename,
+ int* pnColors)
+{
+ VSILFILE* fpColorFile = VSIFOpenL(pszColorFilename, "rt");
+ if (fpColorFile == NULL)
+ {
+ CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszColorFilename);
+ *pnColors = 0;
+ return NULL;
+ }
+
+ ColorAssociation* pasColorAssociation = NULL;
+ int nColorAssociation = 0;
+
+ int bSrcHasNoData = FALSE;
+ double dfSrcNoDataValue = GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
+
+ const char* pszLine;
+ int bIsGMT_CPT = FALSE;
+ while ((pszLine = CPLReadLineL(fpColorFile)) != NULL)
+ {
+ if (pszLine[0] == '#' && strstr(pszLine, "COLOR_MODEL"))
+ {
+ if (strstr(pszLine, "COLOR_MODEL = RGB") == NULL)
+ {
+ CPLError(CE_Failure, CPLE_AppDefined, "Only COLOR_MODEL = RGB is supported");
+ CPLFree(pasColorAssociation);
+ *pnColors = 0;
+ return NULL;
+ }
+ bIsGMT_CPT = TRUE;
+ }
+
+ char** papszFields = CSLTokenizeStringComplex(pszLine, " ,\t:",
+ FALSE, FALSE );
+ /* Skip comment lines */
+ int nTokens = CSLCount(papszFields);
+ if (nTokens >= 1 && (papszFields[0][0] == '#' ||
+ papszFields[0][0] == '/'))
+ {
+ CSLDestroy(papszFields);
+ continue;
+ }
+
+ if (bIsGMT_CPT && nTokens == 8)
+ {
+ pasColorAssociation =
+ (ColorAssociation*)CPLRealloc(pasColorAssociation,
+ (nColorAssociation + 2) * sizeof(ColorAssociation));
+
+ pasColorAssociation[nColorAssociation].dfVal = atof(papszFields[0]);
+ pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]);
+ pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]);
+ pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]);
+ pasColorAssociation[nColorAssociation].nA = 255;
+ nColorAssociation++;
+
+ pasColorAssociation[nColorAssociation].dfVal = atof(papszFields[4]);
+ pasColorAssociation[nColorAssociation].nR = atoi(papszFields[5]);
+ pasColorAssociation[nColorAssociation].nG = atoi(papszFields[6]);
+ pasColorAssociation[nColorAssociation].nB = atoi(papszFields[7]);
+ pasColorAssociation[nColorAssociation].nA = 255;
+ nColorAssociation++;
+ }
+ else if (bIsGMT_CPT && nTokens == 4)
+ {
+ /* The first token might be B (background), F (foreground) or N (nodata) */
+ /* Just interested in N */
+ if (EQUAL(papszFields[0], "N") && bSrcHasNoData)
+ {
+ pasColorAssociation =
+ (ColorAssociation*)CPLRealloc(pasColorAssociation,
+ (nColorAssociation + 1) * sizeof(ColorAssociation));
+
+ pasColorAssociation[nColorAssociation].dfVal = dfSrcNoDataValue;
+ pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]);
+ pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]);
+ pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]);
+ pasColorAssociation[nColorAssociation].nA = 255;
+ nColorAssociation++;
+ }
+ }
+ else if (!bIsGMT_CPT && nTokens >= 2)
+ {
+ pasColorAssociation =
+ (ColorAssociation*)CPLRealloc(pasColorAssociation,
+ (nColorAssociation + 1) * sizeof(ColorAssociation));
+ if (EQUAL(papszFields[0], "nv") && bSrcHasNoData)
+ pasColorAssociation[nColorAssociation].dfVal = dfSrcNoDataValue;
+ else if (strlen(papszFields[0]) > 1 && papszFields[0][strlen(papszFields[0])-1] == '%')
+ {
+ double dfPct = atof(papszFields[0]) / 100.;
+ if (dfPct < 0.0 || dfPct > 1.0)
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "Wrong value for a percentage : %s", papszFields[0]);
+ CSLDestroy(papszFields);
+ VSIFCloseL(fpColorFile);
+ CPLFree(pasColorAssociation);
+ *pnColors = 0;
+ return NULL;
+ }
+ pasColorAssociation[nColorAssociation].dfVal =
+ GDALColorReliefGetAbsoluteValFromPct(hSrcBand, dfPct);
+ }
+ else
+ pasColorAssociation[nColorAssociation].dfVal = atof(papszFields[0]);
+
+ if (nTokens >= 4)
+ {
+ pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]);
+ pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]);
+ pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]);
+ pasColorAssociation[nColorAssociation].nA =
+ (CSLCount(papszFields) >= 5 ) ? atoi(papszFields[4]) : 255;
+ }
+ else
+ {
+ int nR, nG, nB;
+ if (!GDALColorReliefFindNamedColor(papszFields[1], &nR, &nG, &nB))
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "Unknown color : %s", papszFields[1]);
+ CSLDestroy(papszFields);
+ VSIFCloseL(fpColorFile);
+ CPLFree(pasColorAssociation);
+ *pnColors = 0;
+ return NULL;
+ }
+ pasColorAssociation[nColorAssociation].nR = nR;
+ pasColorAssociation[nColorAssociation].nG = nG;
+ pasColorAssociation[nColorAssociation].nB = nB;
+ pasColorAssociation[nColorAssociation].nA =
+ (CSLCount(papszFields) >= 3 ) ? atoi(papszFields[2]) : 255;
+ }
+
+ nColorAssociation ++;
+ }
+ CSLDestroy(papszFields);
+ }
+ VSIFCloseL(fpColorFile);
+
+ if (nColorAssociation == 0)
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "No color association found in %s", pszColorFilename);
+ *pnColors = 0;
+ return NULL;
+ }
+
+ qsort(pasColorAssociation, nColorAssociation,
+ sizeof(ColorAssociation), GDALColorReliefSortColors);
+
+ *pnColors = nColorAssociation;
+ return pasColorAssociation;
+}
+
+static
+GByte* GDALColorReliefPrecompute(GDALRasterBandH hSrcBand,
+ ColorAssociation* pasColorAssociation,
+ int nColorAssociation,
+ ColorSelectionMode eColorSelectionMode,
+ int* pnIndexOffset)
+{
+ GDALDataType eDT = GDALGetRasterDataType(hSrcBand);
+ GByte* pabyPrecomputed = NULL;
+ int nIndexOffset = (eDT == GDT_Int16) ? 32768 : 0;
+ *pnIndexOffset = nIndexOffset;
+ int nXSize = GDALGetRasterBandXSize(hSrcBand);
+ int nYSize = GDALGetRasterBandXSize(hSrcBand);
+ if (eDT == GDT_Byte ||
+ ((eDT == GDT_Int16 || eDT == GDT_UInt16) && nXSize * nYSize > 65536))
+ {
+ int iMax = (eDT == GDT_Byte) ? 256: 65536;
+ pabyPrecomputed = (GByte*) VSIMalloc(4 * iMax);
+ if (pabyPrecomputed)
+ {
+ int i;
+ for(i=0;ihSrcDS = hSrcDS;
+ this->hSrcBand = hSrcBand;
+ nColorAssociation = 0;
+ pasColorAssociation =
+ GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
+ &nColorAssociation);
+ this->eColorSelectionMode = eColorSelectionMode;
+
+ nRasterXSize = GDALGetRasterXSize(hSrcDS);
+ nRasterYSize = GDALGetRasterYSize(hSrcDS);
+
+ int nBlockXSize, nBlockYSize;
+ GDALGetBlockSize( hSrcBand, &nBlockXSize, &nBlockYSize);
+
+ nIndexOffset = 0;
+ pabyPrecomputed =
+ GDALColorReliefPrecompute(hSrcBand,
+ pasColorAssociation,
+ nColorAssociation,
+ eColorSelectionMode,
+ &nIndexOffset);
+
+ int i;
+ for(i=0;i<((bAlpha) ? 4 : 3);i++)
+ {
+ SetBand(i + 1, new GDALColorReliefRasterBand(this, i+1));
+ }
+
+ pafSourceBuf = NULL;
+ panSourceBuf = NULL;
+ if (pabyPrecomputed)
+ panSourceBuf = (int *) CPLMalloc(sizeof(int)*nBlockXSize*nBlockYSize);
+ else
+ pafSourceBuf = (float *) CPLMalloc(sizeof(float)*nBlockXSize*nBlockYSize);
+ nCurBlockXOff = -1;
+ nCurBlockYOff = -1;
+}
+
+GDALColorReliefDataset::~GDALColorReliefDataset()
+{
+ CPLFree(pasColorAssociation);
+ CPLFree(pabyPrecomputed);
+ CPLFree(panSourceBuf);
+ CPLFree(pafSourceBuf);
+}
+
+CPLErr GDALColorReliefDataset::GetGeoTransform( double * padfGeoTransform )
+{
+ return GDALGetGeoTransform(hSrcDS, padfGeoTransform);
+}
+
+const char *GDALColorReliefDataset::GetProjectionRef()
+{
+ return GDALGetProjectionRef(hSrcDS);
+}
+
+GDALColorReliefRasterBand::GDALColorReliefRasterBand(
+ GDALColorReliefDataset * poDS, int nBand)
+{
+ this->poDS = poDS;
+ this->nBand = nBand;
+ eDataType = GDT_Byte;
+ GDALGetBlockSize( poDS->hSrcBand, &nBlockXSize, &nBlockYSize);
+}
+
+CPLErr GDALColorReliefRasterBand::IReadBlock( int nBlockXOff,
+ int nBlockYOff,
+ void *pImage )
+{
+ GDALColorReliefDataset * poGDS = (GDALColorReliefDataset *) poDS;
+ int nReqXSize, nReqYSize;
+
+ if ((nBlockXOff + 1) * nBlockXSize >= nRasterXSize)
+ nReqXSize = nRasterXSize - nBlockXOff * nBlockXSize;
+ else
+ nReqXSize = nBlockXSize;
+
+ if ((nBlockYOff + 1) * nBlockYSize >= nRasterYSize)
+ nReqYSize = nRasterYSize - nBlockYOff * nBlockYSize;
+ else
+ nReqYSize = nBlockYSize;
+
+ if ( poGDS->nCurBlockXOff != nBlockXOff ||
+ poGDS->nCurBlockYOff != nBlockYOff )
+ {
+ poGDS->nCurBlockXOff = nBlockXOff;
+ poGDS->nCurBlockYOff = nBlockYOff;
+
+ CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
+ GF_Read,
+ nBlockXOff * nBlockXSize,
+ nBlockYOff * nBlockYSize,
+ nReqXSize, nReqYSize,
+ (poGDS->panSourceBuf) ?
+ (void*) poGDS->panSourceBuf :
+ (void* )poGDS->pafSourceBuf,
+ nReqXSize, nReqYSize,
+ (poGDS->panSourceBuf) ? GDT_Int32 : GDT_Float32,
+ 0, 0);
+ if (eErr != CE_None)
+ {
+ memset(pImage, 0, nBlockXSize * nBlockYSize);
+ return eErr;
+ }
+ }
+
+ int x, y, j = 0;
+ if (poGDS->panSourceBuf)
+ {
+ for( y = 0; y < nReqYSize; y++ )
+ {
+ for( x = 0; x < nReqXSize; x++ )
+ {
+ int nIndex = poGDS->panSourceBuf[j] + poGDS->nIndexOffset;
+ ((GByte*)pImage)[y * nBlockXSize + x] = poGDS->pabyPrecomputed[4*nIndex + nBand-1];
+ j++;
+ }
+ }
+ }
+ else
+ {
+ int anComponents[4];
+ for( y = 0; y < nReqYSize; y++ )
+ {
+ for( x = 0; x < nReqXSize; x++ )
+ {
+ GDALColorReliefGetRGBA (poGDS->pasColorAssociation,
+ poGDS->nColorAssociation,
+ poGDS->pafSourceBuf[j],
+ poGDS->eColorSelectionMode,
+ &anComponents[0],
+ &anComponents[1],
+ &anComponents[2],
+ &anComponents[3]);
+ ((GByte*)pImage)[y * nBlockXSize + x] = (GByte) anComponents[nBand-1];
+ j++;
+ }
+ }
+ }
+
+ return CE_None;
+}
+
+GDALColorInterp GDALColorReliefRasterBand::GetColorInterpretation()
+{
+ return (GDALColorInterp)(GCI_RedBand + nBand - 1);
+}
+
+
+CPLErr GDALColorRelief (GDALRasterBandH hSrcBand,
+ GDALRasterBandH hDstBand1,
+ GDALRasterBandH hDstBand2,
+ GDALRasterBandH hDstBand3,
+ GDALRasterBandH hDstBand4,
+ const char* pszColorFilename,
+ ColorSelectionMode eColorSelectionMode,
+ GDALProgressFunc pfnProgress,
+ void * pProgressData)
+{
+ CPLErr eErr;
+
+ if (hSrcBand == NULL || hDstBand1 == NULL || hDstBand2 == NULL ||
+ hDstBand3 == NULL)
+ return CE_Failure;
+
+ int nColorAssociation = 0;
+ ColorAssociation* pasColorAssociation =
+ GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
+ &nColorAssociation);
+ if (pasColorAssociation == NULL)
+ return CE_Failure;
+
+ int nXSize = GDALGetRasterBandXSize(hSrcBand);
+ int nYSize = GDALGetRasterBandYSize(hSrcBand);
+
+ if (pfnProgress == NULL)
+ pfnProgress = GDALDummyProgress;
+
+ int nR = 0, nG = 0, nB = 0, nA = 0;
+
+/* -------------------------------------------------------------------- */
+/* Precompute the map from values to RGBA quadruplets */
+/* for GDT_Byte, GDT_Int16 or GDT_UInt16 */
+/* -------------------------------------------------------------------- */
+ int nIndexOffset = 0;
+ GByte* pabyPrecomputed =
+ GDALColorReliefPrecompute(hSrcBand,
+ pasColorAssociation,
+ nColorAssociation,
+ eColorSelectionMode,
+ &nIndexOffset);
+
+/* -------------------------------------------------------------------- */
+/* Initialize progress counter. */
+/* -------------------------------------------------------------------- */
+
+ float* pafSourceBuf = NULL;
+ int* panSourceBuf = NULL;
+ if (pabyPrecomputed)
+ panSourceBuf = (int *) CPLMalloc(sizeof(int)*nXSize);
+ else
+ pafSourceBuf = (float *) CPLMalloc(sizeof(float)*nXSize);
+ GByte* pabyDestBuf1 = (GByte*) CPLMalloc( 4 * nXSize );
+ GByte* pabyDestBuf2 = pabyDestBuf1 + nXSize;
+ GByte* pabyDestBuf3 = pabyDestBuf2 + nXSize;
+ GByte* pabyDestBuf4 = pabyDestBuf3 + nXSize;
+ int i, j;
+
+ if( !pfnProgress( 0.0, NULL, pProgressData ) )
+ {
+ CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+ eErr = CE_Failure;
+ goto end;
+ }
+
+ for ( i = 0; i < nYSize; i++)
+ {
+ /* Read source buffer */
+ eErr = GDALRasterIO( hSrcBand,
+ GF_Read,
+ 0, i,
+ nXSize, 1,
+ (panSourceBuf) ? (void*) panSourceBuf : (void* )pafSourceBuf,
+ nXSize, 1,
+ (panSourceBuf) ? GDT_Int32 : GDT_Float32,
+ 0, 0);
+ if (eErr != CE_None)
+ goto end;
+
+ if (panSourceBuf)
+ {
+ for ( j = 0; j < nXSize; j++)
+ {
+ int nIndex = panSourceBuf[j] + nIndexOffset;
+ pabyDestBuf1[j] = pabyPrecomputed[4 * nIndex];
+ pabyDestBuf2[j] = pabyPrecomputed[4 * nIndex + 1];
+ pabyDestBuf3[j] = pabyPrecomputed[4 * nIndex + 2];
+ pabyDestBuf4[j] = pabyPrecomputed[4 * nIndex + 3];
+ }
+ }
+ else
+ {
+ for ( j = 0; j < nXSize; j++)
+ {
+ GDALColorReliefGetRGBA (pasColorAssociation,
+ nColorAssociation,
+ pafSourceBuf[j],
+ eColorSelectionMode,
+ &nR,
+ &nG,
+ &nB,
+ &nA);
+ pabyDestBuf1[j] = (GByte) nR;
+ pabyDestBuf2[j] = (GByte) nG;
+ pabyDestBuf3[j] = (GByte) nB;
+ pabyDestBuf4[j] = (GByte) nA;
+ }
+ }
+
+ /* -----------------------------------------
+ * Write Line to Raster
+ */
+ eErr = GDALRasterIO(hDstBand1,
+ GF_Write,
+ 0, i, nXSize,
+ 1, pabyDestBuf1, nXSize, 1, GDT_Byte, 0, 0);
+ if (eErr != CE_None)
+ goto end;
+
+ eErr = GDALRasterIO(hDstBand2,
+ GF_Write,
+ 0, i, nXSize,
+ 1, pabyDestBuf2, nXSize, 1, GDT_Byte, 0, 0);
+ if (eErr != CE_None)
+ goto end;
+
+ eErr = GDALRasterIO(hDstBand3,
+ GF_Write,
+ 0, i, nXSize,
+ 1, pabyDestBuf3, nXSize, 1, GDT_Byte, 0, 0);
+ if (eErr != CE_None)
+ goto end;
+
+ if (hDstBand4)
+ {
+ eErr = GDALRasterIO(hDstBand4,
+ GF_Write,
+ 0, i, nXSize,
+ 1, pabyDestBuf4, nXSize, 1, GDT_Byte, 0, 0);
+ if (eErr != CE_None)
+ goto end;
+ }
+
+ if( !pfnProgress( 1.0 * (i+1) / nYSize, NULL, pProgressData ) )
+ {
+ CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+ eErr = CE_Failure;
+ goto end;
+ }
+ }
+ pfnProgress( 1.0, NULL, pProgressData );
+ eErr = CE_None;
+
+end:
+ VSIFree(pabyPrecomputed);
+ CPLFree(pafSourceBuf);
+ CPLFree(panSourceBuf);
+ CPLFree(pabyDestBuf1);
+ CPLFree(pasColorAssociation);
+
+ return eErr;
+}
+
+/************************************************************************/
+/* GDALGenerateVRTColorRelief() */
+/************************************************************************/
+
+CPLErr GDALGenerateVRTColorRelief(const char* pszDstFilename,
+ GDALDatasetH hSrcDataset,
+ GDALRasterBandH hSrcBand,
+ const char* pszColorFilename,
+ ColorSelectionMode eColorSelectionMode,
+ int bAddAlpha)
+{
+
+ int nColorAssociation = 0;
+ ColorAssociation* pasColorAssociation =
+ GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
+ &nColorAssociation);
+ if (pasColorAssociation == NULL)
+ return CE_Failure;
+
+ int nXSize = GDALGetRasterBandXSize(hSrcBand);
+ int nYSize = GDALGetRasterBandYSize(hSrcBand);
+
+ VSILFILE* fp = VSIFOpenL(pszDstFilename, "wt");
+ if (fp == NULL)
+ {
+ CPLFree(pasColorAssociation);
+ return CE_Failure;
+ }
+
+ VSIFPrintfL(fp, "\n", nXSize, nYSize);
+ const char* pszProjectionRef = GDALGetProjectionRef(hSrcDataset);
+ if (pszProjectionRef && pszProjectionRef[0] != '\0')
+ {
+ char* pszEscapedString = CPLEscapeString(pszProjectionRef, -1, CPLES_XML);
+ VSIFPrintfL(fp, " %s\n", pszEscapedString);
+ VSIFree(pszEscapedString);
+ }
+ double adfGT[6];
+ if (GDALGetGeoTransform(hSrcDataset, adfGT) == CE_None)
+ {
+ VSIFPrintfL(fp, " %.16g, %.16g, %.16g, "
+ "%.16g, %.16g, %.16g\n",
+ adfGT[0], adfGT[1], adfGT[2], adfGT[3], adfGT[4], adfGT[5]);
+ }
+ int nBands = 3 + (bAddAlpha ? 1 : 0);
+ int iBand;
+
+ int nBlockXSize, nBlockYSize;
+ GDALGetBlockSize(hSrcBand, &nBlockXSize, &nBlockYSize);
+
+ int bRelativeToVRT;
+ CPLString osPath = CPLGetPath(pszDstFilename);
+ char* pszSourceFilename = CPLStrdup(
+ CPLExtractRelativePath( osPath.c_str(), GDALGetDescription(hSrcDataset),
+ &bRelativeToVRT ));
+
+ for(iBand = 0; iBand < nBands; iBand++)
+ {
+ VSIFPrintfL(fp, " \n", iBand + 1);
+ VSIFPrintfL(fp, " %s\n",
+ GDALGetColorInterpretationName((GDALColorInterp)(GCI_RedBand + iBand)));
+ VSIFPrintfL(fp, " \n");
+ VSIFPrintfL(fp, " %s\n",
+ bRelativeToVRT, pszSourceFilename);
+ VSIFPrintfL(fp, " %d\n", GDALGetBandNumber(hSrcBand));
+ VSIFPrintfL(fp, " \n",
+ nXSize, nYSize,
+ GDALGetDataTypeName(GDALGetRasterDataType(hSrcBand)),
+ nBlockXSize, nBlockYSize);
+ VSIFPrintfL(fp, " \n",
+ nXSize, nYSize);
+ VSIFPrintfL(fp, " \n",
+ nXSize, nYSize);
+
+ VSIFPrintfL(fp, " ");
+ int iColor;
+#define EPSILON 1e-8
+ for(iColor=0;iColor 1)
+ VSIFPrintfL(fp, ",");
+ }
+ else if (iColor > 0)
+ VSIFPrintfL(fp, ",");
+
+ double dfVal = pasColorAssociation[iColor].dfVal;
+
+ if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
+ {
+ VSIFPrintfL(fp, "%.12g:0,", dfVal - EPSILON);
+ }
+ else if (iColor > 0 &&
+ eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY)
+ {
+ double dfMidVal = (dfVal + pasColorAssociation[iColor-1].dfVal) / 2;
+ VSIFPrintfL(fp, "%.12g:%d", dfMidVal - EPSILON,
+ (iBand == 0) ? pasColorAssociation[iColor-1].nR :
+ (iBand == 1) ? pasColorAssociation[iColor-1].nG :
+ (iBand == 2) ? pasColorAssociation[iColor-1].nB :
+ pasColorAssociation[iColor-1].nA);
+ VSIFPrintfL(fp, ",%.12g:%d", dfMidVal ,
+ (iBand == 0) ? pasColorAssociation[iColor].nR :
+ (iBand == 1) ? pasColorAssociation[iColor].nG :
+ (iBand == 2) ? pasColorAssociation[iColor].nB :
+ pasColorAssociation[iColor].nA);
+
+ }
+
+ if (eColorSelectionMode != COLOR_SELECTION_NEAREST_ENTRY)
+ {
+ if (dfVal != (double)(int)dfVal)
+ VSIFPrintfL(fp, "%.12g", dfVal);
+ else
+ VSIFPrintfL(fp, "%d", (int)dfVal);
+ VSIFPrintfL(fp, ":%d",
+ (iBand == 0) ? pasColorAssociation[iColor].nR :
+ (iBand == 1) ? pasColorAssociation[iColor].nG :
+ (iBand == 2) ? pasColorAssociation[iColor].nB :
+ pasColorAssociation[iColor].nA);
+ }
+
+ if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
+ {
+ VSIFPrintfL(fp, ",%.12g:0", dfVal + EPSILON);
+ }
+
+ }
+ VSIFPrintfL(fp, "\n");
+
+ VSIFPrintfL(fp, " \n");
+ VSIFPrintfL(fp, " \n");
+ }
+
+ CPLFree(pszSourceFilename);
+
+ VSIFPrintfL(fp, "\n");
+
+ VSIFCloseL(fp);
+
+ CPLFree(pasColorAssociation);
+
+ return CE_None;
+}
+
+
+/************************************************************************/
+/* GDALTRIAlg() */
+/************************************************************************/
+
+float GDALTRIAlg (float* afWin, float fDstNoDataValue, void* pData)
+{
+ // Terrain Ruggedness is average difference in height
+ return (fabs(afWin[0]-afWin[4]) +
+ fabs(afWin[1]-afWin[4]) +
+ fabs(afWin[2]-afWin[4]) +
+ fabs(afWin[3]-afWin[4]) +
+ fabs(afWin[5]-afWin[4]) +
+ fabs(afWin[6]-afWin[4]) +
+ fabs(afWin[7]-afWin[4]) +
+ fabs(afWin[8]-afWin[4]))/8;
+}
+
+
+/************************************************************************/
+/* GDALTPIAlg() */
+/************************************************************************/
+
+float GDALTPIAlg (float* afWin, float fDstNoDataValue, void* pData)
+{
+ // Terrain Position is the difference between
+ // The central cell and the mean of the surrounding cells
+ return afWin[4] -
+ ((afWin[0]+
+ afWin[1]+
+ afWin[2]+
+ afWin[3]+
+ afWin[5]+
+ afWin[6]+
+ afWin[7]+
+ afWin[8])/8);
+}
+
+/************************************************************************/
+/* GDALRoughnessAlg() */
+/************************************************************************/
+
+float GDALRoughnessAlg (float* afWin, float fDstNoDataValue, void* pData)
+{
+ // Roughness is the largest difference
+ // between any two cells
+
+ float pafRoughnessMin = afWin[0];
+ float pafRoughnessMax = afWin[0];
+
+ for ( int k = 1; k < 9; k++)
+ {
+ if (afWin[k] > pafRoughnessMax)
+ {
+ pafRoughnessMax=afWin[k];
+ }
+ if (afWin[k] < pafRoughnessMin)
+ {
+ pafRoughnessMin=afWin[k];
+ }
+ }
+ return pafRoughnessMax - pafRoughnessMin;
+}
+
+/************************************************************************/
+/* ==================================================================== */
+/* GDALGeneric3x3Dataset */
+/* ==================================================================== */
+/************************************************************************/
+
+class GDALGeneric3x3RasterBand;
+
+class GDALGeneric3x3Dataset : public GDALDataset
+{
+ friend class GDALGeneric3x3RasterBand;
+
+ GDALGeneric3x3ProcessingAlg pfnAlg;
+ void* pAlgData;
+ GDALDatasetH hSrcDS;
+ GDALRasterBandH hSrcBand;
+ float* apafSourceBuf[3];
+ int bDstHasNoData;
+ double dfDstNoDataValue;
+ int nCurLine;
+ int bComputeAtEdges;
+
+ public:
+ GDALGeneric3x3Dataset(GDALDatasetH hSrcDS,
+ GDALRasterBandH hSrcBand,
+ GDALDataType eDstDataType,
+ int bDstHasNoData,
+ double dfDstNoDataValue,
+ GDALGeneric3x3ProcessingAlg pfnAlg,
+ void* pAlgData,
+ int bComputeAtEdges);
+ ~GDALGeneric3x3Dataset();
+
+ CPLErr GetGeoTransform( double * padfGeoTransform );
+ const char *GetProjectionRef();
+};
+
+/************************************************************************/
+/* ==================================================================== */
+/* GDALGeneric3x3RasterBand */
+/* ==================================================================== */
+/************************************************************************/
+
+class GDALGeneric3x3RasterBand : public GDALRasterBand
+{
+ friend class GDALGeneric3x3Dataset;
+ int bSrcHasNoData;
+ float fSrcNoDataValue;
+
+ void InitWidthNoData(void* pImage);
+
+ public:
+ GDALGeneric3x3RasterBand( GDALGeneric3x3Dataset *poDS,
+ GDALDataType eDstDataType );
+
+ virtual CPLErr IReadBlock( int, int, void * );
+ virtual double GetNoDataValue( int* pbHasNoData );
+};
+
+GDALGeneric3x3Dataset::GDALGeneric3x3Dataset(
+ GDALDatasetH hSrcDS,
+ GDALRasterBandH hSrcBand,
+ GDALDataType eDstDataType,
+ int bDstHasNoData,
+ double dfDstNoDataValue,
+ GDALGeneric3x3ProcessingAlg pfnAlg,
+ void* pAlgData,
+ int bComputeAtEdges)
+{
+ this->hSrcDS = hSrcDS;
+ this->hSrcBand = hSrcBand;
+ this->pfnAlg = pfnAlg;
+ this->pAlgData = pAlgData;
+ this->bDstHasNoData = bDstHasNoData;
+ this->dfDstNoDataValue = dfDstNoDataValue;
+ this->bComputeAtEdges = bComputeAtEdges;
+
+ CPLAssert(eDstDataType == GDT_Byte || eDstDataType == GDT_Float32);
+
+ nRasterXSize = GDALGetRasterXSize(hSrcDS);
+ nRasterYSize = GDALGetRasterYSize(hSrcDS);
+
+ SetBand(1, new GDALGeneric3x3RasterBand(this, eDstDataType));
+
+ apafSourceBuf[0] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
+ apafSourceBuf[1] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
+ apafSourceBuf[2] = (float *) CPLMalloc(sizeof(float)*nRasterXSize);
+
+ nCurLine = -1;
+}
+
+GDALGeneric3x3Dataset::~GDALGeneric3x3Dataset()
+{
+ CPLFree(apafSourceBuf[0]);
+ CPLFree(apafSourceBuf[1]);
+ CPLFree(apafSourceBuf[2]);
+}
+
+CPLErr GDALGeneric3x3Dataset::GetGeoTransform( double * padfGeoTransform )
+{
+ return GDALGetGeoTransform(hSrcDS, padfGeoTransform);
+}
+
+const char *GDALGeneric3x3Dataset::GetProjectionRef()
+{
+ return GDALGetProjectionRef(hSrcDS);
+}
+
+GDALGeneric3x3RasterBand::GDALGeneric3x3RasterBand(GDALGeneric3x3Dataset *poDS,
+ GDALDataType eDstDataType)
+{
+ this->poDS = poDS;
+ this->nBand = 1;
+ eDataType = eDstDataType;
+ nBlockXSize = poDS->GetRasterXSize();
+ nBlockYSize = 1;
+
+ bSrcHasNoData = FALSE;
+ fSrcNoDataValue = (float)GDALGetRasterNoDataValue(poDS->hSrcBand,
+ &bSrcHasNoData);
+}
+
+void GDALGeneric3x3RasterBand::InitWidthNoData(void* pImage)
+{
+ int j;
+ GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
+ if (eDataType == GDT_Byte)
+ {
+ for(j=0;jdfDstNoDataValue;
+ }
+ else
+ {
+ for(j=0;jdfDstNoDataValue;
+ }
+}
+
+CPLErr GDALGeneric3x3RasterBand::IReadBlock( int nBlockXOff,
+ int nBlockYOff,
+ void *pImage )
+{
+ int i, j;
+ float fVal;
+ GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
+
+ if (poGDS->bComputeAtEdges && nRasterXSize >= 2 && nRasterYSize >= 2)
+ {
+ if (nBlockYOff == 0)
+ {
+ for(i=0;i<2;i++)
+ {
+ CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
+ GF_Read,
+ 0, i, nBlockXSize, 1,
+ poGDS->apafSourceBuf[i+1],
+ nBlockXSize, 1,
+ GDT_Float32,
+ 0, 0);
+ if (eErr != CE_None)
+ {
+ InitWidthNoData(pImage);
+ return eErr;
+ }
+ }
+ poGDS->nCurLine = 0;
+
+ for (j = 0; j < nRasterXSize; j++)
+ {
+ float afWin[9];
+ int jmin = (j == 0) ? j : j - 1;
+ int jmax = (j == nRasterXSize - 1) ? j : j + 1;
+
+ afWin[0] = INTERPOL(poGDS->apafSourceBuf[1][jmin], poGDS->apafSourceBuf[2][jmin]);
+ afWin[1] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[2][j]);
+ afWin[2] = INTERPOL(poGDS->apafSourceBuf[1][jmax], poGDS->apafSourceBuf[2][jmax]);
+ afWin[3] = poGDS->apafSourceBuf[1][jmin];
+ afWin[4] = poGDS->apafSourceBuf[1][j];
+ afWin[5] = poGDS->apafSourceBuf[1][jmax];
+ afWin[6] = poGDS->apafSourceBuf[2][jmin];
+ afWin[7] = poGDS->apafSourceBuf[2][j];
+ afWin[8] = poGDS->apafSourceBuf[2][jmax];
+
+ fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
+ afWin, (float) poGDS->dfDstNoDataValue,
+ poGDS->pfnAlg,
+ poGDS->pAlgData,
+ poGDS->bComputeAtEdges);
+
+ if (eDataType == GDT_Byte)
+ ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
+ else
+ ((float*)pImage)[j] = fVal;
+ }
+
+ return CE_None;
+ }
+ else if (nBlockYOff == nRasterYSize - 1)
+ {
+ if (poGDS->nCurLine != nRasterYSize - 2)
+ {
+ for(i=0;i<2;i++)
+ {
+ CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
+ GF_Read,
+ 0, nRasterYSize - 2 + i, nBlockXSize, 1,
+ poGDS->apafSourceBuf[i+1],
+ nBlockXSize, 1,
+ GDT_Float32,
+ 0, 0);
+ if (eErr != CE_None)
+ {
+ InitWidthNoData(pImage);
+ return eErr;
+ }
+ }
+ }
+
+ for (j = 0; j < nRasterXSize; j++)
+ {
+ float afWin[9];
+ int jmin = (j == 0) ? j : j - 1;
+ int jmax = (j == nRasterXSize - 1) ? j : j + 1;
+
+ afWin[0] = poGDS->apafSourceBuf[1][jmin];
+ afWin[1] = poGDS->apafSourceBuf[1][j];
+ afWin[2] = poGDS->apafSourceBuf[1][jmax];
+ afWin[3] = poGDS->apafSourceBuf[2][jmin];
+ afWin[4] = poGDS->apafSourceBuf[2][j];
+ afWin[5] = poGDS->apafSourceBuf[2][jmax];
+ afWin[6] = INTERPOL(poGDS->apafSourceBuf[2][jmin], poGDS->apafSourceBuf[1][jmin]);
+ afWin[7] = INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[1][j]);
+ afWin[8] = INTERPOL(poGDS->apafSourceBuf[2][jmax], poGDS->apafSourceBuf[1][jmax]);
+
+ fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
+ afWin, (float) poGDS->dfDstNoDataValue,
+ poGDS->pfnAlg,
+ poGDS->pAlgData,
+ poGDS->bComputeAtEdges);
+
+ if (eDataType == GDT_Byte)
+ ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
+ else
+ ((float*)pImage)[j] = fVal;
+ }
+
+ return CE_None;
+ }
+ }
+ else if ( nBlockYOff == 0 || nBlockYOff == nRasterYSize - 1)
+ {
+ InitWidthNoData(pImage);
+ return CE_None;
+ }
+
+ if ( poGDS->nCurLine != nBlockYOff )
+ {
+ if (poGDS->nCurLine + 1 == nBlockYOff)
+ {
+ float* pafTmp = poGDS->apafSourceBuf[0];
+ poGDS->apafSourceBuf[0] = poGDS->apafSourceBuf[1];
+ poGDS->apafSourceBuf[1] = poGDS->apafSourceBuf[2];
+ poGDS->apafSourceBuf[2] = pafTmp;
+
+ CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
+ GF_Read,
+ 0, nBlockYOff + 1, nBlockXSize, 1,
+ poGDS->apafSourceBuf[2],
+ nBlockXSize, 1,
+ GDT_Float32,
+ 0, 0);
+
+ if (eErr != CE_None)
+ {
+ InitWidthNoData(pImage);
+ return eErr;
+ }
+ }
+ else
+ {
+ for(i=0;i<3;i++)
+ {
+ CPLErr eErr = GDALRasterIO( poGDS->hSrcBand,
+ GF_Read,
+ 0, nBlockYOff + i - 1, nBlockXSize, 1,
+ poGDS->apafSourceBuf[i],
+ nBlockXSize, 1,
+ GDT_Float32,
+ 0, 0);
+ if (eErr != CE_None)
+ {
+ InitWidthNoData(pImage);
+ return eErr;
+ }
+ }
+ }
+
+ poGDS->nCurLine = nBlockYOff;
+ }
+
+ if (poGDS->bComputeAtEdges && nRasterXSize >= 2)
+ {
+ float afWin[9];
+
+ j = 0;
+ afWin[0] = INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j+1]);
+ afWin[1] = poGDS->apafSourceBuf[0][j];
+ afWin[2] = poGDS->apafSourceBuf[0][j+1];
+ afWin[3] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j+1]);
+ afWin[4] = poGDS->apafSourceBuf[1][j];
+ afWin[5] = poGDS->apafSourceBuf[1][j+1];
+ afWin[6] = INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j+1]);
+ afWin[7] = poGDS->apafSourceBuf[2][j];
+ afWin[8] = poGDS->apafSourceBuf[2][j+1];
+
+ fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
+ afWin, (float) poGDS->dfDstNoDataValue,
+ poGDS->pfnAlg,
+ poGDS->pAlgData,
+ poGDS->bComputeAtEdges);
+
+ if (eDataType == GDT_Byte)
+ ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
+ else
+ ((float*)pImage)[j] = fVal;
+
+ j = nRasterXSize - 1;
+
+ afWin[0] = poGDS->apafSourceBuf[0][j-1];
+ afWin[1] = poGDS->apafSourceBuf[0][j];
+ afWin[2] = INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j-1]);
+ afWin[3] = poGDS->apafSourceBuf[1][j-1];
+ afWin[4] = poGDS->apafSourceBuf[1][j];
+ afWin[5] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j-1]);
+ afWin[6] = poGDS->apafSourceBuf[2][j-1];
+ afWin[7] = poGDS->apafSourceBuf[2][j];
+ afWin[8] = INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j-1]);
+
+ fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
+ afWin, (float) poGDS->dfDstNoDataValue,
+ poGDS->pfnAlg,
+ poGDS->pAlgData,
+ poGDS->bComputeAtEdges);
+
+ if (eDataType == GDT_Byte)
+ ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
+ else
+ ((float*)pImage)[j] = fVal;
+ }
+ else
+ {
+ if (eDataType == GDT_Byte)
+ {
+ ((GByte*)pImage)[0] = (GByte) poGDS->dfDstNoDataValue;
+ if (nBlockXSize > 1)
+ ((GByte*)pImage)[nBlockXSize - 1] = (GByte) poGDS->dfDstNoDataValue;
+ }
+ else
+ {
+ ((float*)pImage)[0] = (float) poGDS->dfDstNoDataValue;
+ if (nBlockXSize > 1)
+ ((float*)pImage)[nBlockXSize - 1] = (float) poGDS->dfDstNoDataValue;
+ }
+ }
+
+
+ for(j=1;japafSourceBuf[0][j-1];
+ afWin[1] = poGDS->apafSourceBuf[0][j];
+ afWin[2] = poGDS->apafSourceBuf[0][j+1];
+ afWin[3] = poGDS->apafSourceBuf[1][j-1];
+ afWin[4] = poGDS->apafSourceBuf[1][j];
+ afWin[5] = poGDS->apafSourceBuf[1][j+1];
+ afWin[6] = poGDS->apafSourceBuf[2][j-1];
+ afWin[7] = poGDS->apafSourceBuf[2][j];
+ afWin[8] = poGDS->apafSourceBuf[2][j+1];
+
+ fVal = ComputeVal(bSrcHasNoData, fSrcNoDataValue,
+ afWin, (float) poGDS->dfDstNoDataValue,
+ poGDS->pfnAlg,
+ poGDS->pAlgData,
+ poGDS->bComputeAtEdges);
+
+ if (eDataType == GDT_Byte)
+ ((GByte*)pImage)[j] = (GByte) (fVal + 0.5);
+ else
+ ((float*)pImage)[j] = fVal;
+
+ }
+
+ return CE_None;
+}
+
+double GDALGeneric3x3RasterBand::GetNoDataValue( int* pbHasNoData )
+{
+ GDALGeneric3x3Dataset * poGDS = (GDALGeneric3x3Dataset *) poDS;
+ if (pbHasNoData)
+ *pbHasNoData = poGDS->bDstHasNoData;
+ return poGDS->dfDstNoDataValue;
+}
+
+/************************************************************************/
+/* main() */
+/************************************************************************/
+
+typedef enum
+{
+ HILL_SHADE,
+ SLOPE,
+ ASPECT,
+ COLOR_RELIEF,
+ TRI,
+ TPI,
+ ROUGHNESS
+} Algorithm;
+
+#define CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(nExtraArg) \
+ do { if (i + nExtraArg >= argc) \
+ Usage(CPLSPrintf("%s option requires %d argument(s)", argv[i], nExtraArg)); } while(0)
+
+#ifndef ZOO_SERVICE
+int main( int argc, char ** argv )
+#else
+#ifdef WIN32
+__declspec(dllexport)
+#endif
+ int Gdal_Dem(maps*& conf,maps*& inputs,maps*& outputs)
+#endif
+{
+ Algorithm eUtilityMode = HILL_SHADE;
+ double z = 1.0;
+ double scale = 1.0;
+ double az = 315.0;
+ double alt = 45.0;
+ // 0 = 'percent' or 1 = 'degrees'
+ int slopeFormat = 1;
+ int bAddAlpha = FALSE;
+ int bZeroForFlat = FALSE;
+ int bAngleAsAzimuth = TRUE;
+ ColorSelectionMode eColorSelectionMode = COLOR_SELECTION_INTERPOLATE;
+
+ int nBand = 1;
+ double adfGeoTransform[6];
+
+ const char *pszSrcFilename = NULL;
+ const char *pszDstFilename = NULL;
+ const char *pszColorFilename = NULL;
+ const char *pszFormat = "GTiff";
+ int bFormatExplicitelySet = FALSE;
+ char **papszCreateOptions = NULL;
+
+ GDALDatasetH hSrcDataset = NULL;
+ GDALDatasetH hDstDataset = NULL;
+ GDALRasterBandH hSrcBand = NULL;
+ GDALRasterBandH hDstBand = NULL;
+ GDALDriverH hDriver = NULL;
+
+ GDALProgressFunc pfnProgress = GDALTermProgress;
+
+ int nXSize = 0;
+ int nYSize = 0;
+
+ int bComputeAtEdges = FALSE;
+ int bZevenbergenThorne = FALSE;
+ int bCombined = FALSE;
+ int bQuiet = FALSE;
+
+#ifndef ZOO_SERVICE
+ /* Check strict compilation and runtime library version as we use C++ API */
+ if (! GDAL_CHECK_VERSION(argv[0]))
+ exit(1);
+
+ argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
+ if( argc < 2 )
+ {
+ Usage("Not enough arguments.");
+ }
+
+ if( EQUAL(argv[1], "--utility_version") || EQUAL(argv[1], "--utility-version") )
+ {
+ printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
+ argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
+ return 0;
+ }
+ else if( EQUAL(argv[1],"--help") )
+ Usage();
+ else if ( EQUAL(argv[1], "shade") || EQUAL(argv[1], "hillshade") )
+ {
+ eUtilityMode = HILL_SHADE;
+ }
+ else if ( EQUAL(argv[1], "slope") )
+ {
+ eUtilityMode = SLOPE;
+ }
+ else if ( EQUAL(argv[1], "aspect") )
+ {
+ eUtilityMode = ASPECT;
+ }
+ else if ( EQUAL(argv[1], "color-relief") )
+ {
+ eUtilityMode = COLOR_RELIEF;
+ }
+ else if ( EQUAL(argv[1], "TRI") )
+ {
+ eUtilityMode = TRI;
+ }
+ else if ( EQUAL(argv[1], "TPI") )
+ {
+ eUtilityMode = TPI;
+ }
+ else if ( EQUAL(argv[1], "roughness") )
+ {
+ eUtilityMode = ROUGHNESS;
+ }
+ else
+ {
+ Usage("Missing valid sub-utility mention.");
+ }
+
+/* -------------------------------------------------------------------- */
+/* Parse arguments. */
+/* -------------------------------------------------------------------- */
+ for(int i = 2; i < argc; i++ )
+ {
+ if( eUtilityMode == HILL_SHADE &&
+ (EQUAL(argv[i], "--z") || EQUAL(argv[i], "-z")))
+ {
+ CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
+ z = atof(argv[++i]);
+ }
+ else if ( eUtilityMode == SLOPE && EQUAL(argv[i], "-p"))
+ {
+ slopeFormat = 0;
+ }
+ else if ( (eUtilityMode == HILL_SHADE || eUtilityMode == SLOPE ||
+ eUtilityMode == ASPECT) && EQUAL(argv[i], "-alg") )
+ {
+ CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
+ i ++;
+ if (EQUAL(argv[i], "ZevenbergenThorne"))
+ bZevenbergenThorne = TRUE;
+ else if (!EQUAL(argv[i], "Horn"))
+ {
+ Usage(CPLSPrintf("Wrong value for alg : %s.", argv[i]));
+ }
+ }
+ else if ( eUtilityMode == ASPECT && EQUAL(argv[i], "-trigonometric"))
+ {
+ bAngleAsAzimuth = FALSE;
+ }
+ else if ( eUtilityMode == ASPECT && EQUAL(argv[i], "-zero_for_flat"))
+ {
+ bZeroForFlat = TRUE;
+ }
+ else if ( eUtilityMode == COLOR_RELIEF && EQUAL(argv[i], "-exact_color_entry"))
+ {
+ eColorSelectionMode = COLOR_SELECTION_EXACT_ENTRY;
+ }
+ else if ( eUtilityMode == COLOR_RELIEF && EQUAL(argv[i], "-nearest_color_entry"))
+ {
+ eColorSelectionMode = COLOR_SELECTION_NEAREST_ENTRY;
+ }
+ else if(
+ (EQUAL(argv[i], "--s") ||
+ EQUAL(argv[i], "-s") ||
+ EQUAL(argv[i], "--scale") ||
+ EQUAL(argv[i], "-scale"))
+ )
+ {
+ CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
+ scale = atof(argv[++i]);
+ }
+ else if( eUtilityMode == HILL_SHADE &&
+ (EQUAL(argv[i], "--az") ||
+ EQUAL(argv[i], "-az") ||
+ EQUAL(argv[i], "--azimuth") ||
+ EQUAL(argv[i], "-azimuth"))
+ )
+ {
+ CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
+ az = atof(argv[++i]);
+ }
+ else if( eUtilityMode == HILL_SHADE &&
+ (EQUAL(argv[i], "--alt") ||
+ EQUAL(argv[i], "-alt") ||
+ EQUAL(argv[i], "--alt") ||
+ EQUAL(argv[i], "-alt"))
+ )
+ {
+ CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
+ alt = atof(argv[++i]);
+ }
+ else if( eUtilityMode == HILL_SHADE &&
+ (EQUAL(argv[i], "-combined") ||
+ EQUAL(argv[i], "--combined"))
+ )
+ {
+ bCombined = TRUE;
+ }
+ else if( eUtilityMode == COLOR_RELIEF &&
+ EQUAL(argv[i], "-alpha"))
+ {
+ bAddAlpha = TRUE;
+ }
+ else if( eUtilityMode != COLOR_RELIEF &&
+ EQUAL(argv[i], "-compute_edges"))
+ {
+ bComputeAtEdges = TRUE;
+ }
+ else if( i + 1 < argc &&
+ (EQUAL(argv[i], "--b") ||
+ EQUAL(argv[i], "-b"))
+ )
+ {
+ nBand = atoi(argv[++i]);
+ }
+ else if ( EQUAL(argv[i], "-q") || EQUAL(argv[i], "-quiet") )
+ {
+ pfnProgress = GDALDummyProgress;
+ bQuiet = TRUE;
+ }
+ else if( EQUAL(argv[i],"-co") )
+ {
+ CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
+ papszCreateOptions = CSLAddString( papszCreateOptions, argv[++i] );
+ }
+ else if( EQUAL(argv[i],"-of") )
+ {
+ CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
+ pszFormat = argv[++i];
+ bFormatExplicitelySet = TRUE;
+ }
+ else if( argv[i][0] == '-' )
+ {
+ Usage(CPLSPrintf("Unkown option name '%s'", argv[i]));
+ }
+ else if( pszSrcFilename == NULL )
+ {
+ pszSrcFilename = argv[i];
+ }
+ else if( eUtilityMode == COLOR_RELIEF && pszColorFilename == NULL )
+ {
+ pszColorFilename = argv[i];
+ }
+ else if( pszDstFilename == NULL )
+ {
+ pszDstFilename = argv[i];
+ }
+ else
+ Usage("Too many command options.");
+ }
+
+ if( pszSrcFilename == NULL )
+ {
+ Usage("Missing source.");
+ }
+ if ( eUtilityMode == COLOR_RELIEF && pszColorFilename == NULL )
+ {
+ Usage("Missing color file.");
+ }
+ if( pszDstFilename == NULL )
+ {
+ Usage("Missing destination.");
+ }
+#else
+ map* tmpMap=NULL;
+
+ pfnProgress = GDALDummyProgress;
+ bQuiet = TRUE;
+
+ tmpMap=NULL;
+ tmpMap=getMapFromMaps(inputs,"InputDSN","value");
+ if(tmpMap!=NULL){
+ pszSrcFilename=(char*)CPLMalloc(sizeof(char)*(strlen(tmpMap->value)+1));
+ sprintf((char*)pszSrcFilename,"%s",tmpMap->value);
+ }
+
+ tmpMap=NULL;
+ tmpMap=getMapFromMaps(inputs,"OutputDSN","value");
+ if(tmpMap!=NULL){
+ pszDstFilename=(char*)CPLMalloc(sizeof(char)*(strlen(tmpMap->value)+1));
+ sprintf((char*)pszDstFilename,"%s",tmpMap->value);
+ setMapInMaps(outputs,"Result","value",tmpMap->value);
+ }
+
+ tmpMap=NULL;
+ tmpMap=getMapFromMaps(inputs,"utility","value");
+ if(tmpMap!=NULL){
+ if ( EQUAL(tmpMap->value, "shade") || EQUAL(tmpMap->value, "hillshade") )
+ {
+ eUtilityMode = HILL_SHADE;
+ map* tmpMap1=getMapFromMaps(inputs,"z","value");
+ if(tmpMap1==NULL)
+ return Usage(conf,"Unable to parse your z value");
+ z = atof(tmpMap1->value);
+
+ tmpMap1=getMapFromMaps(inputs,"az","value");
+ if(tmpMap1!=NULL && strcasecmp(tmpMap1->value,"NULL")!=0 )
+ az = atof(tmpMap1->value);
+
+ tmpMap1=getMapFromMaps(inputs,"alt","value");
+ if(tmpMap1!=NULL && strcasecmp(tmpMap1->value,"NULL")!=0 )
+ alt = atof(tmpMap1->value);
+
+ tmpMap1=getMapFromMaps(inputs,"c","value");
+ if(tmpMap1!=NULL && strcasecmp(tmpMap1->value,"true")==0 )
+ bCombined = TRUE;
+ }
+ else if ( EQUAL(tmpMap->value, "slope") )
+ {
+ eUtilityMode = SLOPE;
+ map* tmpMap1=getMapFromMaps(inputs,"p","value");
+ if(tmpMap1!=NULL && strcasecmp(tmpMap1->value,"true")==0)
+ slopeFormat = 0;
+ }
+ else if ( EQUAL(tmpMap->value, "aspect") )
+ {
+ eUtilityMode = ASPECT;
+ map* tmpMap1=getMapFromMaps(inputs,"ece","value");
+ if(tmpMap1!=NULL && strcasecmp(tmpMap1->value,"true")==0)
+ eColorSelectionMode = COLOR_SELECTION_EXACT_ENTRY;
+ tmpMap1=getMapFromMaps(inputs,"nce","value");
+ if(tmpMap1!=NULL && strcasecmp(tmpMap1->value,"true")==0)
+ eColorSelectionMode = COLOR_SELECTION_NEAREST_ENTRY;
+ }
+ else if ( EQUAL(tmpMap->value, "color-relief") )
+ {
+ eUtilityMode = COLOR_RELIEF;
+ map* tmpMap1=getMapFromMaps(inputs,"cfn","value");
+ if(tmpMap!=NULL){
+ pszColorFilename=strdup(tmpMap1->value);
+ }
+ tmpMap1=getMapFromMaps(inputs,"a","value");
+ if(tmpMap1!=NULL && strcasecmp(tmpMap1->value,"true")==0)
+ bAddAlpha = TRUE;
+ tmpMap1=getMapFromMaps(inputs,"ce","value");
+ if(tmpMap1!=NULL && strcasecmp(tmpMap1->value,"true")==0)
+ bComputeAtEdges = TRUE;
+ }
+ else if ( EQUAL(tmpMap->value, "TRI") )
+ {
+ eUtilityMode = TRI;
+ }
+ else if ( EQUAL(tmpMap->value, "TPI") )
+ {
+ eUtilityMode = TPI;
+ }
+ else if ( EQUAL(tmpMap->value, "roughness") )
+ {
+ eUtilityMode = ROUGHNESS;
+ }
+ else
+ {
+ return Usage(conf,"Missing valid utility mention.");
+ }
+ }
+ map* tmpMap1=getMapFromMaps(inputs,"s","value");
+ if(tmpMap1!=NULL && strcasecmp(tmpMap1->value,"NULL")!=0 )
+ scale = atof(tmpMap1->value);
+ map* tmpMap1=getMapFromMaps(inputs,"b","value");
+ if(tmpMap1!=NULL && strcasecmp(tmpMap1->value,"NULL")!=0 )
+ nBand = atoi(tmpMap1->value);
+
+ if( pszSrcFilename == NULL )
+ {
+ return Usage(conf,"Missing source.");
+ }
+ if ( eUtilityMode == COLOR_RELIEF && pszColorFilename == NULL )
+ {
+ return Usage(conf,"Missing color file.");
+ }
+ if( pszDstFilename == NULL )
+ {
+ return Usage(conf,"Missing destination.");
+ }
+#endif
+
+ GDALAllRegister();
+
+ // Open Dataset and get raster band
+ hSrcDataset = GDALOpen( pszSrcFilename, GA_ReadOnly );
+
+ if( hSrcDataset == NULL )
+ {
+ fprintf( stderr,
+ "GDALOpen failed - %d\n%s\n",
+ CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
+ GDALDestroyDriverManager();
+#ifdef ZOO_SERVICE
+ char *tmp=(char*)malloc(1024*sizeof(char));
+ sprintf(tmp,"GDALOpen failed - %d\n%s\n",
+ CPLGetLastErrorNo(), CPLGetLastErrorMsg());
+ return Usage(conf,"GDALOpen failed");
+#else
+ exit( 1 );
+#endif
+ }
+
+ nXSize = GDALGetRasterXSize(hSrcDataset);
+ nYSize = GDALGetRasterYSize(hSrcDataset);
+
+ if( nBand <= 0 || nBand > GDALGetRasterCount(hSrcDataset) )
+ {
+ fprintf( stderr,
+ "Unable to fetch band #%d\n", nBand );
+ GDALDestroyDriverManager();
+#ifdef ZOO_SERVICE
+ char *tmp=(char*)malloc(128*sizeof(char));
+ sprintf(tmp,"Unable to fetch band #%d\n", nBand );
+ return Usage(conf,tmp);
+#else
+ exit( 1 );
+#endif
+ }
+ hSrcBand = GDALGetRasterBand( hSrcDataset, nBand );
+
+ GDALGetGeoTransform(hSrcDataset, adfGeoTransform);
+
+ hDriver = GDALGetDriverByName(pszFormat);
+ if( hDriver == NULL
+ || (GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL &&
+ GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, NULL ) == NULL))
+ {
+ int iDr;
+
+ fprintf( stderr, "Output driver `%s' not recognised to have\n",
+ pszFormat );
+ fprintf( stderr, "output support. The following format drivers are configured\n"
+ "and support output:\n" );
+
+ for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
+ {
+ GDALDriverH hDriver = GDALGetDriver(iDr);
+
+ if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) != NULL ||
+ GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, NULL ) != NULL )
+ {
+ printf( " %s: %s\n",
+ GDALGetDriverShortName( hDriver ),
+ GDALGetDriverLongName( hDriver ) );
+ }
+ }
+ GDALDestroyDriverManager();
+#ifdef ZOO_SERVICE
+ return Usage(conf,"Output driver not recognised to have output support.");
+#else
+ exit( 1 );
+#endif
+ }
+
+#ifndef ZOO_SERVICE
+ if (!bQuiet && !bFormatExplicitelySet)
+ CheckExtensionConsistency(pszDstFilename, pszFormat);
+#endif
+
+ double dfDstNoDataValue = 0;
+ int bDstHasNoData = FALSE;
+ void* pData = NULL;
+ GDALGeneric3x3ProcessingAlg pfnAlg = NULL;
+
+ if (eUtilityMode == HILL_SHADE)
+ {
+ dfDstNoDataValue = 0;
+ bDstHasNoData = TRUE;
+ pData = GDALCreateHillshadeData (adfGeoTransform,
+ z,
+ scale,
+ alt,
+ az,
+ bZevenbergenThorne);
+ if (bZevenbergenThorne)
+ {
+ if(!bCombined)
+ pfnAlg = GDALHillshadeZevenbergenThorneAlg;
+ else
+ pfnAlg = GDALHillshadeZevenbergenThorneCombinedAlg;
+ }
+ else
+ {
+ if(!bCombined)
+ pfnAlg = GDALHillshadeAlg;
+ else
+ pfnAlg = GDALHillshadeCombinedAlg;
+ }
+ }
+ else if (eUtilityMode == SLOPE)
+ {
+ dfDstNoDataValue = -9999;
+ bDstHasNoData = TRUE;
+
+ pData = GDALCreateSlopeData(adfGeoTransform, scale, slopeFormat);
+ if (bZevenbergenThorne)
+ pfnAlg = GDALSlopeZevenbergenThorneAlg;
+ else
+ pfnAlg = GDALSlopeHornAlg;
+ }
+
+ else if (eUtilityMode == ASPECT)
+ {
+ if (!bZeroForFlat)
+ {
+ dfDstNoDataValue = -9999;
+ bDstHasNoData = TRUE;
+ }
+
+ pData = GDALCreateAspectData(bAngleAsAzimuth);
+ if (bZevenbergenThorne)
+ pfnAlg = GDALAspectZevenbergenThorneAlg;
+ else
+ pfnAlg = GDALAspectAlg;
+ }
+ else if (eUtilityMode == TRI)
+ {
+ dfDstNoDataValue = -9999;
+ bDstHasNoData = TRUE;
+ pfnAlg = GDALTRIAlg;
+ }
+ else if (eUtilityMode == TPI)
+ {
+ dfDstNoDataValue = -9999;
+ bDstHasNoData = TRUE;
+ pfnAlg = GDALTPIAlg;
+ }
+ else if (eUtilityMode == ROUGHNESS)
+ {
+ dfDstNoDataValue = -9999;
+ bDstHasNoData = TRUE;
+ pfnAlg = GDALRoughnessAlg;
+ }
+
+ GDALDataType eDstDataType = (eUtilityMode == HILL_SHADE ||
+ eUtilityMode == COLOR_RELIEF) ? GDT_Byte :
+ GDT_Float32;
+
+ if( EQUAL(pszFormat, "VRT") )
+ {
+ if (eUtilityMode == COLOR_RELIEF)
+ {
+ GDALGenerateVRTColorRelief(pszDstFilename,
+ hSrcDataset,
+ hSrcBand,
+ pszColorFilename,
+ eColorSelectionMode,
+ bAddAlpha);
+ GDALClose(hSrcDataset);
+
+ CPLFree(pData);
+
+ GDALDestroyDriverManager();
+#ifndef ZOO_SERVICE
+ CSLDestroy( argv );
+#endif
+ CSLDestroy( papszCreateOptions );
+#ifdef ZOO_SERVICE
+ return SERVICE_SUCCEEDED;
+#else
+ return 0;
+#endif
+ }
+ else
+ {
+ fprintf(stderr, "VRT driver can only be used with color-relief utility\n");
+ GDALDestroyDriverManager();
+#ifdef ZOO_SERVICE
+ return Usage(conf,"VRT driver can only be used with color-relief utility");
+#else
+ exit(1);
+#endif
+ }
+ }
+
+ if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL &&
+ GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, NULL ) != NULL)
+ {
+ GDALDatasetH hIntermediateDataset;
+
+ if (eUtilityMode == COLOR_RELIEF)
+ hIntermediateDataset = (GDALDatasetH)
+ new GDALColorReliefDataset (hSrcDataset,
+ hSrcBand,
+ pszColorFilename,
+ eColorSelectionMode,
+ bAddAlpha);
+ else
+ hIntermediateDataset = (GDALDatasetH)
+ new GDALGeneric3x3Dataset(hSrcDataset, hSrcBand,
+ eDstDataType,
+ bDstHasNoData,
+ dfDstNoDataValue,
+ pfnAlg,
+ pData,
+ bComputeAtEdges);
+
+ GDALDatasetH hOutDS = GDALCreateCopy(
+ hDriver, pszDstFilename, hIntermediateDataset,
+ TRUE, papszCreateOptions,
+ pfnProgress, NULL );
+
+ if( hOutDS != NULL )
+ GDALClose( hOutDS );
+ GDALClose(hIntermediateDataset);
+ GDALClose(hSrcDataset);
+
+ CPLFree(pData);
+
+ GDALDestroyDriverManager();
+#ifndef ZOO_SERVICE
+ CSLDestroy( argv );
+#endif
+ CSLDestroy( papszCreateOptions );
+
+#ifdef ZOO_SERVICE
+ return SERVICE_SUCCEEDED;
+#else
+ return 0;
+#endif
+ }
+
+ int nDstBands;
+ if (eUtilityMode == COLOR_RELIEF)
+ nDstBands = (bAddAlpha) ? 4 : 3;
+ else
+ nDstBands = 1;
+
+ hDstDataset = GDALCreate( hDriver,
+ pszDstFilename,
+ nXSize,
+ nYSize,
+ nDstBands,
+ eDstDataType,
+ papszCreateOptions);
+
+ if( hDstDataset == NULL )
+ {
+ fprintf( stderr,
+ "Unable to create dataset %s %d\n%s\n", pszDstFilename,
+ CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
+ GDALDestroyDriverManager();
+#ifdef ZOO_SERVICE
+ char *tmp=(char*)malloc(1024*sizeof(char));
+ sprintf( tmp,
+ "Unable to create dataset %s %d\n%s\n", pszDstFilename,
+ CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
+ return Usage(conf,tmp);
+#else
+ exit( 1 );
+#endif
+ }
+
+ hDstBand = GDALGetRasterBand( hDstDataset, 1 );
+
+ GDALSetGeoTransform(hDstDataset, adfGeoTransform);
+ GDALSetProjection(hDstDataset, GDALGetProjectionRef(hSrcDataset));
+
+ if (eUtilityMode == COLOR_RELIEF)
+ {
+ GDALColorRelief (hSrcBand,
+ GDALGetRasterBand(hDstDataset, 1),
+ GDALGetRasterBand(hDstDataset, 2),
+ GDALGetRasterBand(hDstDataset, 3),
+ (bAddAlpha) ? GDALGetRasterBand(hDstDataset, 4) : NULL,
+ pszColorFilename,
+ eColorSelectionMode,
+ pfnProgress, NULL);
+ }
+ else
+ {
+ if (bDstHasNoData)
+ GDALSetRasterNoDataValue(hDstBand, dfDstNoDataValue);
+
+ GDALGeneric3x3Processing(hSrcBand, hDstBand,
+ pfnAlg, pData,
+ bComputeAtEdges,
+ pfnProgress, NULL);
+
+ }
+
+ GDALClose(hSrcDataset);
+ GDALClose(hDstDataset);
+ CPLFree(pData);
+
+ GDALDestroyDriverManager();
+#ifdef ZOO_SERVICE
+ CSLDestroy( papszCreateOptions );
+ return SERVICE_SUCCEEDED;
+#else
+ CSLDestroy( argv );
+ CSLDestroy( papszCreateOptions );
+ return 0;
+#endif
+
+}
+
+#ifdef ZOO_SERVICE
+}
+#endif