service.c @ 767

Last change on this file since 767 was 444, checked in by djay, 11 years ago
Fixing issue on Linux / Mac OS X about language gesture. Fixing strange issue in some services using MapServer? support on WIN32. Use the MapServer? nmake.opt to make sure to get all the definitions settled correctly.
File size: 98.8 KB

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

Note: See TracBrowser for help on using the repository browser.

Open WPS platform

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

Download in other formats:

Download now !

Support ZOO-Project

Search

Context Navigation

Main Navigation

ZOO Sponsors

Knowledge partners

Related links