service.c @ 699

Last change on this file since 699 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

Line
1	/******************************************************************************
2	* $Id$
3	*
4	* Project: GDAL DEM Utilities
5	* Purpose:
6	* Authors: Matthew Perry, perrygeo at gmail.com
7	* Even Rouault, even dot rouault at mines dash paris dot org
8	* Howard Butler, hobu.inc at gmail.com
9	* Chris Yesson, chris dot yesson at ioz dot ac dot uk
10	*
11	******************************************************************************
12	* Copyright (c) 2006, 2009 Matthew Perry
13	* Copyright (c) 2009 Even Rouault
14	* Portions derived from GRASS 4.1 (public domain) See
15	* http://trac.osgeo.org/gdal/ticket/2975 for more information regarding
16	* history of this code
17	*
18	* Permission is hereby granted, free of charge, to any person obtaining a
19	* copy of this software and associated documentation files (the "Software"),
20	* to deal in the Software without restriction, including without limitation
21	* the rights to use, copy, modify, merge, publish, distribute, sublicense,
22	* and/or sell copies of the Software, and to permit persons to whom the
23	* Software is furnished to do so, subject to the following conditions:
24	*
25	* The above copyright notice and this permission notice shall be included
26	* in all copies or substantial portions of the Software.
27	*
28	* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
29	* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
30	* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
31	* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
32	* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
33	* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
34	* DEALINGS IN THE SOFTWARE.
35	****************************************************************************
36	*
37	* Slope and aspect calculations based on original method for GRASS GIS 4.1
38	* by Michael Shapiro, U.S.Army Construction Engineering Research Laboratory
39	* Olga Waupotitsch, U.S.Army Construction Engineering Research Laboratory
40	* Marjorie Larson, U.S.Army Construction Engineering Research Laboratory
41	* as found in GRASS's r.slope.aspect module.
42	*
43	* Horn's formula is used to find the first order derivatives in x and y directions
44	* for slope and aspect calculations: Horn, B. K. P. (1981).
45	* "Hill Shading and the Reflectance Map", Proceedings of the IEEE, 69(1):14-47.
46	*
47	* Other reference :
48	* Burrough, P.A. and McDonell, R.A., 1998. Principles of Geographical Information
49	* Systems. p. 190.
50	*
51	* Shaded relief based on original method for GRASS GIS 4.1 by Jim Westervelt,
52	* U.S. Army Construction Engineering Research Laboratory
53	* as found in GRASS's r.shaded.relief (formerly shade.rel.sh) module.
54	* ref: "r.mapcalc: An Algebra for GIS and Image Processing",
55	* by Michael Shapiro and Jim Westervelt, U.S. Army Construction Engineering
56	* Research Laboratory (March/1991)
57	*
58	* Color table of named colors and lookup code derived from src/libes/gis/named_colr.c
59	* of GRASS 4.1
60	*
61	* TRI - Terrain Ruggedness Index is as descibed in Wilson et al (2007)
62	* this is based on the method of Valentine et al. (2004)
63	*
64	* TPI - Topographic Position Index follows the description in Wilson et al (2007), following Weiss (2001)
65	* The radius is fixed at 1 cell width/height
66	*
67	* Roughness - follows the definition in Wilson et al. (2007), which follows Dartnell (2000)
68	*
69	* References for TRI/TPI/Roughness:
70	* Dartnell, P. 2000. Applying Remote Sensing Techniques to map Seafloor
71	* Geology/Habitat Relationships. Masters Thesis, San Francisco State
72	* University, pp. 108.
73	* Valentine, P. C., S. J. Fuller, L. A. Scully. 2004. Terrain Ruggedness
74	* Analysis and Distribution of Boulder Ridges in the Stellwagen Bank National
75	* Marine Sanctuary Region (poster). Galway, Ireland: 5th International
76	* Symposium on Marine Geological and Biological Habitat Mapping (GeoHAB),
77	* May 2004.
78	* Weiss, A. D. 2001. Topographic Positions and Landforms Analysis (poster),
79	* ESRI International User Conference, July 2001. San Diego, CA: ESRI.
80	* Wilson, M. F. J.; O'Connell, B.; Brown, C.; Guinan, J. C. & Grehan, A. J.
81	* Multiscale terrain analysis of multibeam bathymetry data for habitat mapping
82	* on the continental slope Marine Geodesy, 2007, 30, 3-35
83	****************************************************************************/
84
85	#include <stdlib.h>
86	#include <math.h>
87
88	#include "cpl_conv.h"
89	#include "cpl_string.h"
90	#include "gdal.h"
91	#include "gdal_priv.h"
92	#ifndef ZOO_SERVICE
93	#include "commonutils.h"
94	#else
95	#include "service.h"
96	#endif
97
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);
2563	tmpMap1=getMapFromMaps(inputs,"b","value");
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 @ 699

Download in other formats:

Download now !

Support ZOO-Project

Search

Context Navigation

Main Navigation

ZOO Sponsors

Knowledge partners

Related links