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