1 : /* ****************************************************************************
2 : * $Id: gdal_grid.cpp 18482 2010-01-08 23:22:58Z mloskot $
3 : *
4 : * Project: GDAL Utilities
5 : * Purpose: GDAL scattered data gridding (interpolation) tool
6 : * Author: Andrey Kiselev, dron@ak4719.spb.edu
7 : *
8 : * ****************************************************************************
9 : * Copyright (c) 2007, Andrey Kiselev <dron@ak4719.spb.edu>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include <cstdlib>
31 : #include <vector>
32 : #include <algorithm>
33 :
34 : #include "cpl_string.h"
35 : #include "gdal.h"
36 : #include "gdal_alg.h"
37 : #include "ogr_spatialref.h"
38 : #include "ogr_api.h"
39 : #include "ogrsf_frmts.h"
40 :
41 : CPL_CVSID("$Id: gdal_grid.cpp 18482 2010-01-08 23:22:58Z mloskot $");
42 :
43 : static const char szAlgNameInvDist[] = "invdist";
44 : static const char szAlgNameAverage[] = "average";
45 : static const char szAlgNameNearest[] = "nearest";
46 : static const char szAlgNameMinimum[] = "minimum";
47 : static const char szAlgNameMaximum[] = "maximum";
48 : static const char szAlgNameRange[] = "range";
49 : static const char szAlgNameCount[] = "count";
50 : static const char szAlgNameAverageDistance[] = "average_distance";
51 : static const char szAlgNameAverageDistancePts[] = "average_distance_pts";
52 :
53 : /************************************************************************/
54 : /* Usage() */
55 : /************************************************************************/
56 :
57 0 : static void Usage()
58 :
59 : {
60 : printf(
61 : "Usage: gdal_grid [--help-general] [--formats]\n"
62 : " [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/\n"
63 : " CInt16/CInt32/CFloat32/CFloat64}]\n"
64 : " [-of format] [-co \"NAME=VALUE\"]\n"
65 : " [-zfield field_name]\n"
66 : " [-a_srs srs_def] [-spat xmin ymin xmax ymax]\n"
67 : " [-clipsrc <xmin ymin xmax ymax>|WKT|datasource|spat_extent]\n"
68 : " [-clipsrcsql sql_statement] [-clipsrclayer layer]\n"
69 : " [-clipsrcwhere expression]\n"
70 : " [-l layername]* [-where expression] [-sql select_statement]\n"
71 : " [-txe xmin xmax] [-tye ymin ymax] [-outsize xsize ysize]\n"
72 : " [-a algorithm[:parameter1=value1]*]"
73 : " [-q]\n"
74 : " <src_datasource> <dst_filename>\n"
75 : "\n"
76 : "Available algorithms and parameters with their's defaults:\n"
77 : " Inverse distance to a power (default)\n"
78 : " invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0\n"
79 : " Moving average\n"
80 : " average:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0\n"
81 : " Nearest neighbor\n"
82 : " nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0\n"
83 : " Various data metrics\n"
84 : " <metric name>:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0\n"
85 : " possible metrics are:\n"
86 : " minimum\n"
87 : " maximum\n"
88 : " range\n"
89 : " count\n"
90 : " average_distance\n"
91 : " average_distance_pts\n"
92 0 : "\n");
93 :
94 0 : GDALDestroyDriverManager();
95 0 : exit( 1 );
96 : }
97 :
98 : /************************************************************************/
99 : /* GetAlgorithmName() */
100 : /* */
101 : /* Translates algortihm code into mnemonic name. */
102 : /************************************************************************/
103 :
104 1 : static void PrintAlgorithmAndOptions( GDALGridAlgorithm eAlgorithm,
105 : void *pOptions )
106 : {
107 1 : switch ( eAlgorithm )
108 : {
109 : case GGA_InverseDistanceToAPower:
110 0 : printf( "Algorithm name: \"%s\".\n", szAlgNameInvDist );
111 : printf( "Options are "
112 : "\"power=%f:smoothing=%f:radius1=%f:radius2=%f:angle=%f"
113 : ":max_points=%lu:min_points=%lu:nodata=%f\"\n",
114 : ((GDALGridInverseDistanceToAPowerOptions *)pOptions)->dfPower,
115 : ((GDALGridInverseDistanceToAPowerOptions *)pOptions)->dfSmoothing,
116 : ((GDALGridInverseDistanceToAPowerOptions *)pOptions)->dfRadius1,
117 : ((GDALGridInverseDistanceToAPowerOptions *)pOptions)->dfRadius2,
118 : ((GDALGridInverseDistanceToAPowerOptions *)pOptions)->dfAngle,
119 : (unsigned long)((GDALGridInverseDistanceToAPowerOptions *)pOptions)->nMaxPoints,
120 : (unsigned long)((GDALGridInverseDistanceToAPowerOptions *)pOptions)->nMinPoints,
121 0 : ((GDALGridInverseDistanceToAPowerOptions *)pOptions)->dfNoDataValue);
122 0 : break;
123 : case GGA_MovingAverage:
124 0 : printf( "Algorithm name: \"%s\".\n", szAlgNameAverage );
125 : printf( "Options are "
126 : "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
127 : ":nodata=%f\"\n",
128 : ((GDALGridMovingAverageOptions *)pOptions)->dfRadius1,
129 : ((GDALGridMovingAverageOptions *)pOptions)->dfRadius2,
130 : ((GDALGridMovingAverageOptions *)pOptions)->dfAngle,
131 : (unsigned long)((GDALGridMovingAverageOptions *)pOptions)->nMinPoints,
132 0 : ((GDALGridMovingAverageOptions *)pOptions)->dfNoDataValue);
133 0 : break;
134 : case GGA_NearestNeighbor:
135 1 : printf( "Algorithm name: \"%s\".\n", szAlgNameNearest );
136 : printf( "Options are "
137 : "\"radius1=%f:radius2=%f:angle=%f:nodata=%f\"\n",
138 : ((GDALGridNearestNeighborOptions *)pOptions)->dfRadius1,
139 : ((GDALGridNearestNeighborOptions *)pOptions)->dfRadius2,
140 : ((GDALGridNearestNeighborOptions *)pOptions)->dfAngle,
141 1 : ((GDALGridNearestNeighborOptions *)pOptions)->dfNoDataValue);
142 1 : break;
143 : case GGA_MetricMinimum:
144 0 : printf( "Algorithm name: \"%s\".\n", szAlgNameMinimum );
145 : printf( "Options are "
146 : "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
147 : ":nodata=%f\"\n",
148 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius1,
149 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius2,
150 : ((GDALGridDataMetricsOptions *)pOptions)->dfAngle,
151 : (unsigned long)((GDALGridDataMetricsOptions *)pOptions)->nMinPoints,
152 0 : ((GDALGridDataMetricsOptions *)pOptions)->dfNoDataValue);
153 0 : break;
154 : case GGA_MetricMaximum:
155 0 : printf( "Algorithm name: \"%s\".\n", szAlgNameMaximum );
156 : printf( "Options are "
157 : "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
158 : ":nodata=%f\"\n",
159 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius1,
160 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius2,
161 : ((GDALGridDataMetricsOptions *)pOptions)->dfAngle,
162 : (unsigned long)((GDALGridDataMetricsOptions *)pOptions)->nMinPoints,
163 0 : ((GDALGridDataMetricsOptions *)pOptions)->dfNoDataValue);
164 0 : break;
165 : case GGA_MetricRange:
166 0 : printf( "Algorithm name: \"%s\".\n", szAlgNameRange );
167 : printf( "Options are "
168 : "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
169 : ":nodata=%f\"\n",
170 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius1,
171 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius2,
172 : ((GDALGridDataMetricsOptions *)pOptions)->dfAngle,
173 : (unsigned long)((GDALGridDataMetricsOptions *)pOptions)->nMinPoints,
174 0 : ((GDALGridDataMetricsOptions *)pOptions)->dfNoDataValue);
175 0 : break;
176 : case GGA_MetricCount:
177 0 : printf( "Algorithm name: \"%s\".\n", szAlgNameCount );
178 : printf( "Options are "
179 : "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
180 : ":nodata=%f\"\n",
181 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius1,
182 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius2,
183 : ((GDALGridDataMetricsOptions *)pOptions)->dfAngle,
184 : (unsigned long)((GDALGridDataMetricsOptions *)pOptions)->nMinPoints,
185 0 : ((GDALGridDataMetricsOptions *)pOptions)->dfNoDataValue);
186 0 : break;
187 : case GGA_MetricAverageDistance:
188 0 : printf( "Algorithm name: \"%s\".\n", szAlgNameAverageDistance );
189 : printf( "Options are "
190 : "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
191 : ":nodata=%f\"\n",
192 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius1,
193 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius2,
194 : ((GDALGridDataMetricsOptions *)pOptions)->dfAngle,
195 : (unsigned long)((GDALGridDataMetricsOptions *)pOptions)->nMinPoints,
196 0 : ((GDALGridDataMetricsOptions *)pOptions)->dfNoDataValue);
197 0 : break;
198 : case GGA_MetricAverageDistancePts:
199 0 : printf( "Algorithm name: \"%s\".\n", szAlgNameAverageDistancePts );
200 : printf( "Options are "
201 : "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
202 : ":nodata=%f\"\n",
203 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius1,
204 : ((GDALGridDataMetricsOptions *)pOptions)->dfRadius2,
205 : ((GDALGridDataMetricsOptions *)pOptions)->dfAngle,
206 : (unsigned long)((GDALGridDataMetricsOptions *)pOptions)->nMinPoints,
207 0 : ((GDALGridDataMetricsOptions *)pOptions)->dfNoDataValue);
208 0 : break;
209 : default:
210 0 : printf( "Algorithm is unknown.\n" );
211 : break;
212 : }
213 1 : }
214 :
215 : /************************************************************************/
216 : /* ParseAlgorithmAndOptions() */
217 : /* */
218 : /* Translates mnemonic gridding algorithm names into */
219 : /* GDALGridAlgorithm code, parse control parameters and assign */
220 : /* defaults. */
221 : /************************************************************************/
222 :
223 1 : static CPLErr ParseAlgorithmAndOptions( const char *pszAlgoritm,
224 : GDALGridAlgorithm *peAlgorithm,
225 : void **ppOptions )
226 : {
227 1 : char **papszParms = CSLTokenizeString2( pszAlgoritm, ":", FALSE );
228 :
229 1 : if ( CSLCount(papszParms) < 1 )
230 0 : return CE_Failure;
231 :
232 1 : if ( EQUAL(papszParms[0], szAlgNameInvDist) )
233 0 : *peAlgorithm = GGA_InverseDistanceToAPower;
234 1 : else if ( EQUAL(papszParms[0], szAlgNameAverage) )
235 0 : *peAlgorithm = GGA_MovingAverage;
236 1 : else if ( EQUAL(papszParms[0], szAlgNameNearest) )
237 1 : *peAlgorithm = GGA_NearestNeighbor;
238 0 : else if ( EQUAL(papszParms[0], szAlgNameMinimum) )
239 0 : *peAlgorithm = GGA_MetricMinimum;
240 0 : else if ( EQUAL(papszParms[0], szAlgNameMaximum) )
241 0 : *peAlgorithm = GGA_MetricMaximum;
242 0 : else if ( EQUAL(papszParms[0], szAlgNameRange) )
243 0 : *peAlgorithm = GGA_MetricRange;
244 0 : else if ( EQUAL(papszParms[0], szAlgNameCount) )
245 0 : *peAlgorithm = GGA_MetricCount;
246 0 : else if ( EQUAL(papszParms[0], szAlgNameAverageDistance) )
247 0 : *peAlgorithm = GGA_MetricAverageDistance;
248 0 : else if ( EQUAL(papszParms[0], szAlgNameAverageDistancePts) )
249 0 : *peAlgorithm = GGA_MetricAverageDistancePts;
250 : else
251 : {
252 : fprintf( stderr, "Unsupported gridding method \"%s\".\n",
253 0 : papszParms[0] );
254 0 : CSLDestroy( papszParms );
255 0 : return CE_Failure;
256 : }
257 :
258 : /* -------------------------------------------------------------------- */
259 : /* Parse algorithm parameters and assign defaults. */
260 : /* -------------------------------------------------------------------- */
261 : const char *pszValue;
262 :
263 1 : switch ( *peAlgorithm )
264 : {
265 : case GGA_InverseDistanceToAPower:
266 : default:
267 : *ppOptions =
268 0 : CPLMalloc( sizeof(GDALGridInverseDistanceToAPowerOptions) );
269 :
270 0 : pszValue = CSLFetchNameValue( papszParms, "power" );
271 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
272 0 : dfPower = (pszValue) ? atof(pszValue) : 2.0;
273 :
274 0 : pszValue = CSLFetchNameValue( papszParms, "smoothing" );
275 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
276 0 : dfSmoothing = (pszValue) ? atof(pszValue) : 0.0;
277 :
278 0 : pszValue = CSLFetchNameValue( papszParms, "radius1" );
279 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
280 0 : dfRadius1 = (pszValue) ? atof(pszValue) : 0.0;
281 :
282 0 : pszValue = CSLFetchNameValue( papszParms, "radius2" );
283 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
284 0 : dfRadius2 = (pszValue) ? atof(pszValue) : 0.0;
285 :
286 0 : pszValue = CSLFetchNameValue( papszParms, "angle" );
287 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
288 0 : dfAngle = (pszValue) ? atof(pszValue) : 0.0;
289 :
290 0 : pszValue = CSLFetchNameValue( papszParms, "max_points" );
291 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
292 0 : nMaxPoints = (pszValue) ? atol(pszValue) : 0;
293 :
294 0 : pszValue = CSLFetchNameValue( papszParms, "min_points" );
295 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
296 0 : nMinPoints = (pszValue) ? atol(pszValue) : 0;
297 :
298 0 : pszValue = CSLFetchNameValue( papszParms, "nodata" );
299 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
300 0 : dfNoDataValue = (pszValue) ? atof(pszValue) : 0.0;
301 0 : break;
302 :
303 : case GGA_MovingAverage:
304 : *ppOptions =
305 0 : CPLMalloc( sizeof(GDALGridMovingAverageOptions) );
306 :
307 0 : pszValue = CSLFetchNameValue( papszParms, "radius1" );
308 : ((GDALGridMovingAverageOptions *)*ppOptions)->
309 0 : dfRadius1 = (pszValue) ? atof(pszValue) : 0.0;
310 :
311 0 : pszValue = CSLFetchNameValue( papszParms, "radius2" );
312 : ((GDALGridMovingAverageOptions *)*ppOptions)->
313 0 : dfRadius2 = (pszValue) ? atof(pszValue) : 0.0;
314 :
315 0 : pszValue = CSLFetchNameValue( papszParms, "angle" );
316 : ((GDALGridMovingAverageOptions *)*ppOptions)->
317 0 : dfAngle = (pszValue) ? atof(pszValue) : 0.0;
318 :
319 0 : pszValue = CSLFetchNameValue( papszParms, "min_points" );
320 : ((GDALGridMovingAverageOptions *)*ppOptions)->
321 0 : nMinPoints = (pszValue) ? atol(pszValue) : 0;
322 :
323 0 : pszValue = CSLFetchNameValue( papszParms, "nodata" );
324 : ((GDALGridMovingAverageOptions *)*ppOptions)->
325 0 : dfNoDataValue = (pszValue) ? atof(pszValue) : 0.0;
326 0 : break;
327 :
328 : case GGA_NearestNeighbor:
329 : *ppOptions =
330 1 : CPLMalloc( sizeof(GDALGridNearestNeighborOptions) );
331 :
332 1 : pszValue = CSLFetchNameValue( papszParms, "radius1" );
333 : ((GDALGridNearestNeighborOptions *)*ppOptions)->
334 1 : dfRadius1 = (pszValue) ? atof(pszValue) : 0.0;
335 :
336 1 : pszValue = CSLFetchNameValue( papszParms, "radius2" );
337 : ((GDALGridNearestNeighborOptions *)*ppOptions)->
338 1 : dfRadius2 = (pszValue) ? atof(pszValue) : 0.0;
339 :
340 1 : pszValue = CSLFetchNameValue( papszParms, "angle" );
341 : ((GDALGridNearestNeighborOptions *)*ppOptions)->
342 1 : dfAngle = (pszValue) ? atof(pszValue) : 0.0;
343 :
344 1 : pszValue = CSLFetchNameValue( papszParms, "nodata" );
345 : ((GDALGridNearestNeighborOptions *)*ppOptions)->
346 1 : dfNoDataValue = (pszValue) ? atof(pszValue) : 0.0;
347 1 : break;
348 :
349 : case GGA_MetricMinimum:
350 : case GGA_MetricMaximum:
351 : case GGA_MetricRange:
352 : case GGA_MetricCount:
353 : case GGA_MetricAverageDistance:
354 : case GGA_MetricAverageDistancePts:
355 : *ppOptions =
356 0 : CPLMalloc( sizeof(GDALGridDataMetricsOptions) );
357 :
358 0 : pszValue = CSLFetchNameValue( papszParms, "radius1" );
359 : ((GDALGridDataMetricsOptions *)*ppOptions)->
360 0 : dfRadius1 = (pszValue) ? atof(pszValue) : 0.0;
361 :
362 0 : pszValue = CSLFetchNameValue( papszParms, "radius2" );
363 : ((GDALGridDataMetricsOptions *)*ppOptions)->
364 0 : dfRadius2 = (pszValue) ? atof(pszValue) : 0.0;
365 :
366 0 : pszValue = CSLFetchNameValue( papszParms, "angle" );
367 : ((GDALGridDataMetricsOptions *)*ppOptions)->
368 0 : dfAngle = (pszValue) ? atof(pszValue) : 0.0;
369 :
370 0 : pszValue = CSLFetchNameValue( papszParms, "min_points" );
371 : ((GDALGridDataMetricsOptions *)*ppOptions)->
372 0 : nMinPoints = (pszValue) ? atol(pszValue) : 0;
373 :
374 0 : pszValue = CSLFetchNameValue( papszParms, "nodata" );
375 : ((GDALGridDataMetricsOptions *)*ppOptions)->
376 0 : dfNoDataValue = (pszValue) ? atof(pszValue) : 0.0;
377 : break;
378 :
379 : }
380 :
381 1 : CSLDestroy( papszParms );
382 1 : return CE_None;
383 : }
384 :
385 : /************************************************************************/
386 : /* ProcessGeometry() */
387 : /* */
388 : /* Extract point coordinates from the geometry reference and set the */
389 : /* Z value as requested. Test whther we are in the clipped region */
390 : /* before processing. */
391 : /************************************************************************/
392 :
393 14641 : static void ProcessGeometry( OGRPoint *poGeom, OGRGeometry *poClipSrc,
394 : int iBurnField, double dfBurnValue,
395 : std::vector<double> &adfX,
396 : std::vector<double> &adfY,
397 : std::vector<double> &adfZ)
398 :
399 : {
400 14641 : if ( poClipSrc && !poGeom->Within(poClipSrc) )
401 0 : return;
402 :
403 14641 : adfX.push_back( poGeom->getX() );
404 14641 : adfY.push_back( poGeom->getY() );
405 14641 : if ( iBurnField < 0 )
406 14641 : adfZ.push_back( poGeom->getZ() );
407 : else
408 0 : adfZ.push_back( dfBurnValue );
409 : }
410 :
411 : /************************************************************************/
412 : /* ProcessLayer() */
413 : /* */
414 : /* Process all the features in a layer selection, collecting */
415 : /* geometries and burn values. */
416 : /************************************************************************/
417 :
418 1 : static void ProcessLayer( OGRLayerH hSrcLayer, GDALDatasetH hDstDS,
419 : OGRGeometry *poClipSrc,
420 : GUInt32 nXSize, GUInt32 nYSize, int nBand,
421 : int& bIsXExtentSet, int& bIsYExtentSet,
422 : double& dfXMin, double& dfXMax,
423 : double& dfYMin, double& dfYMax,
424 : const char *pszBurnAttribute,
425 : GDALDataType eType,
426 : GDALGridAlgorithm eAlgorithm, void *pOptions,
427 : int bQuiet, GDALProgressFunc pfnProgress )
428 :
429 : {
430 : /* -------------------------------------------------------------------- */
431 : /* Get field index, and check. */
432 : /* -------------------------------------------------------------------- */
433 1 : int iBurnField = -1;
434 :
435 1 : if ( pszBurnAttribute )
436 : {
437 : iBurnField = OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( hSrcLayer ),
438 0 : pszBurnAttribute );
439 0 : if( iBurnField == -1 )
440 : {
441 : printf( "Failed to find field %s on layer %s, skipping.\n",
442 : pszBurnAttribute,
443 0 : OGR_FD_GetName( OGR_L_GetLayerDefn( hSrcLayer ) ) );
444 0 : return;
445 : }
446 : }
447 :
448 : /* -------------------------------------------------------------------- */
449 : /* Collect the geometries from this layer, and build list of */
450 : /* values to be interpolated. */
451 : /* -------------------------------------------------------------------- */
452 : OGRFeature *poFeat;
453 1 : std::vector<double> adfX, adfY, adfZ;
454 :
455 1 : OGR_L_ResetReading( hSrcLayer );
456 :
457 14643 : while( (poFeat = (OGRFeature *)OGR_L_GetNextFeature( hSrcLayer )) != NULL )
458 : {
459 14641 : OGRGeometry *poGeom = poFeat->GetGeometryRef();
460 :
461 14641 : if ( poGeom != NULL )
462 : {
463 14641 : OGRwkbGeometryType eType = wkbFlatten( poGeom->getGeometryType() );
464 14641 : double dfBurnValue = 0.0;
465 :
466 14641 : if ( iBurnField >= 0 )
467 0 : dfBurnValue = poFeat->GetFieldAsDouble( iBurnField );
468 :
469 14641 : if ( eType == wkbMultiPoint )
470 : {
471 : int iGeom;
472 0 : int nGeomCount = ((OGRMultiPoint *)poGeom)->getNumGeometries();
473 :
474 0 : for ( iGeom = 0; iGeom < nGeomCount; iGeom++ )
475 : {
476 : ProcessGeometry( (OGRPoint *)((OGRMultiPoint *)poGeom)->getGeometryRef(iGeom),
477 : poClipSrc, iBurnField, dfBurnValue,
478 0 : adfX, adfY, adfZ);
479 : }
480 : }
481 14641 : else if ( eType == wkbPoint )
482 : {
483 : ProcessGeometry( (OGRPoint *)poGeom, poClipSrc,
484 14641 : iBurnField, dfBurnValue, adfX, adfY, adfZ);
485 : }
486 : }
487 :
488 14641 : OGRFeature::DestroyFeature( poFeat );
489 : }
490 :
491 1 : if (adfX.size() == 0)
492 : {
493 : printf( "No point geometry found on layer %s, skipping.\n",
494 0 : OGR_FD_GetName( OGR_L_GetLayerDefn( hSrcLayer ) ) );
495 : return;
496 : }
497 :
498 : /* -------------------------------------------------------------------- */
499 : /* Compute grid geometry. */
500 : /* -------------------------------------------------------------------- */
501 :
502 1 : if ( !bIsXExtentSet )
503 : {
504 0 : dfXMin = *std::min_element(adfX.begin(), adfX.end());
505 0 : dfXMax = *std::max_element(adfX.begin(), adfX.end());
506 0 : bIsXExtentSet = TRUE;
507 : }
508 :
509 1 : if ( !bIsYExtentSet )
510 : {
511 0 : dfYMin = *std::min_element(adfY.begin(), adfY.end());
512 0 : dfYMax = *std::max_element(adfY.begin(), adfY.end());
513 0 : bIsYExtentSet = TRUE;
514 : }
515 :
516 : /* -------------------------------------------------------------------- */
517 : /* Perform gridding. */
518 : /* -------------------------------------------------------------------- */
519 :
520 1 : const double dfDeltaX = ( dfXMax - dfXMin ) / nXSize;
521 1 : const double dfDeltaY = ( dfYMax - dfYMin ) / nYSize;
522 :
523 1 : if ( !bQuiet )
524 : {
525 1 : printf( "Grid data type is \"%s\"\n", GDALGetDataTypeName(eType) );
526 : printf( "Grid size = (%lu %lu).\n",
527 1 : (unsigned long)nXSize, (unsigned long)nYSize );
528 : printf( "Corner coordinates = (%f %f)-(%f %f).\n",
529 : dfXMin - dfDeltaX / 2, dfYMax + dfDeltaY / 2,
530 1 : dfXMax + dfDeltaX / 2, dfYMin - dfDeltaY / 2 );
531 1 : printf( "Grid cell size = (%f %f).\n", dfDeltaX, dfDeltaY );
532 1 : printf( "Source point count = %lu.\n", (unsigned long)adfX.size() );
533 1 : PrintAlgorithmAndOptions( eAlgorithm, pOptions );
534 1 : printf("\n");
535 : }
536 :
537 1 : GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, nBand );
538 :
539 1 : if (adfX.size() == 0)
540 : {
541 : // FIXME: Shoulda' set to nodata value instead
542 0 : GDALFillRaster( hBand, 0.0 , 0.0 );
543 : return;
544 : }
545 :
546 : GUInt32 nXOffset, nYOffset;
547 : int nBlockXSize, nBlockYSize;
548 :
549 1 : GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );
550 : void *pData =
551 1 : CPLMalloc( nBlockXSize * nBlockYSize * GDALGetDataTypeSize(eType) );
552 :
553 1 : GUInt32 nBlock = 0;
554 : GUInt32 nBlockCount = ((nXSize + nBlockXSize - 1) / nBlockXSize)
555 1 : * ((nYSize + nBlockYSize - 1) / nBlockYSize);
556 :
557 2 : for ( nYOffset = 0; nYOffset < nYSize; nYOffset += nBlockYSize )
558 : {
559 2 : for ( nXOffset = 0; nXOffset < nXSize; nXOffset += nBlockXSize )
560 : {
561 : void *pScaledProgress;
562 : pScaledProgress =
563 : GDALCreateScaledProgress( 0.0,
564 : (double)++nBlock / nBlockCount,
565 1 : pfnProgress, NULL );
566 :
567 1 : int nXRequest = nBlockXSize;
568 1 : if (nXOffset + nXRequest > nXSize)
569 1 : nXRequest = nXSize - nXOffset;
570 :
571 1 : int nYRequest = nBlockYSize;
572 1 : if (nYOffset + nYRequest > nYSize)
573 1 : nYRequest = nYSize - nYOffset;
574 :
575 : GDALGridCreate( eAlgorithm, pOptions,
576 : adfX.size(), &(adfX[0]), &(adfY[0]), &(adfZ[0]),
577 : dfXMin + dfDeltaX * nXOffset,
578 : dfXMin + dfDeltaX * (nXOffset + nXRequest),
579 : dfYMin + dfDeltaY * nYOffset,
580 : dfYMin + dfDeltaY * (nYOffset + nYRequest),
581 : nXRequest, nYRequest, eType, pData,
582 1 : GDALScaledProgress, pScaledProgress );
583 :
584 : GDALRasterIO( hBand, GF_Write, nXOffset, nYOffset,
585 : nXRequest, nYRequest, pData,
586 1 : nXRequest, nYRequest, eType, 0, 0 );
587 :
588 1 : GDALDestroyScaledProgress( pScaledProgress );
589 : }
590 : }
591 :
592 1 : CPLFree( pData );
593 : }
594 :
595 : /************************************************************************/
596 : /* LoadGeometry() */
597 : /* */
598 : /* Read geometries from the given dataset using specified filters and */
599 : /* returns a collection of read geometries. */
600 : /************************************************************************/
601 :
602 0 : static OGRGeometryCollection* LoadGeometry( const char* pszDS,
603 : const char* pszSQL,
604 : const char* pszLyr,
605 : const char* pszWhere )
606 : {
607 : OGRDataSource *poDS;
608 : OGRLayer *poLyr;
609 : OGRFeature *poFeat;
610 0 : OGRGeometryCollection *poGeom = NULL;
611 :
612 0 : poDS = OGRSFDriverRegistrar::Open( pszDS, FALSE );
613 0 : if ( poDS == NULL )
614 0 : return NULL;
615 :
616 0 : if ( pszSQL != NULL )
617 0 : poLyr = poDS->ExecuteSQL( pszSQL, NULL, NULL );
618 0 : else if ( pszLyr != NULL )
619 0 : poLyr = poDS->GetLayerByName( pszLyr );
620 : else
621 0 : poLyr = poDS->GetLayer(0);
622 :
623 0 : if ( poLyr == NULL )
624 : {
625 : fprintf( stderr,
626 0 : "FAILURE: Failed to identify source layer from datasource.\n" );
627 0 : OGRDataSource::DestroyDataSource( poDS );
628 0 : return NULL;
629 : }
630 :
631 0 : if ( pszWhere )
632 0 : poLyr->SetAttributeFilter( pszWhere );
633 :
634 0 : while ( (poFeat = poLyr->GetNextFeature()) != NULL )
635 : {
636 0 : OGRGeometry* poSrcGeom = poFeat->GetGeometryRef();
637 0 : if ( poSrcGeom )
638 : {
639 : OGRwkbGeometryType eType =
640 0 : wkbFlatten( poSrcGeom->getGeometryType() );
641 :
642 0 : if ( poGeom == NULL )
643 0 : poGeom = new OGRMultiPolygon();
644 :
645 0 : if ( eType == wkbPolygon )
646 0 : poGeom->addGeometry( poSrcGeom );
647 0 : else if ( eType == wkbMultiPolygon )
648 : {
649 : int iGeom;
650 : int nGeomCount =
651 0 : ((OGRMultiPolygon *)poSrcGeom)->getNumGeometries();
652 :
653 0 : for ( iGeom = 0; iGeom < nGeomCount; iGeom++ )
654 : {
655 : poGeom->addGeometry(
656 0 : ((OGRMultiPolygon *)poSrcGeom)->getGeometryRef(iGeom) );
657 : }
658 : }
659 : else
660 : {
661 0 : fprintf( stderr, "FAILURE: Geometry not of polygon type.\n" );
662 0 : OGRGeometryFactory::destroyGeometry( poGeom );
663 0 : OGRFeature::DestroyFeature( poFeat );
664 0 : if ( pszSQL != NULL )
665 0 : poDS->ReleaseResultSet( poLyr );
666 0 : OGRDataSource::DestroyDataSource( poDS );
667 0 : return NULL;
668 : }
669 : }
670 :
671 0 : OGRFeature::DestroyFeature( poFeat );
672 : }
673 :
674 0 : if( pszSQL != NULL )
675 0 : poDS->ReleaseResultSet( poLyr );
676 0 : OGRDataSource::DestroyDataSource( poDS );
677 :
678 0 : return poGeom;
679 : }
680 :
681 : /************************************************************************/
682 : /* main() */
683 : /************************************************************************/
684 :
685 2 : int main( int argc, char ** argv )
686 : {
687 : GDALDriverH hDriver;
688 2 : const char *pszSource=NULL, *pszDest=NULL, *pszFormat = "GTiff";
689 2 : char **papszLayers = NULL;
690 2 : const char *pszBurnAttribute = NULL;
691 2 : const char *pszWHERE = NULL, *pszSQL = NULL;
692 2 : GDALDataType eOutputType = GDT_Float64;
693 2 : char **papszCreateOptions = NULL;
694 2 : GUInt32 nXSize = 0, nYSize = 0;
695 2 : double dfXMin = 0.0, dfXMax = 0.0, dfYMin = 0.0, dfYMax = 0.0;
696 2 : int bIsXExtentSet = FALSE, bIsYExtentSet = FALSE;
697 2 : GDALGridAlgorithm eAlgorithm = GGA_InverseDistanceToAPower;
698 2 : void *pOptions = NULL;
699 2 : char *pszOutputSRS = NULL;
700 2 : int bQuiet = FALSE;
701 2 : GDALProgressFunc pfnProgress = GDALTermProgress;
702 : int i;
703 2 : OGRGeometry *poSpatialFilter = NULL;
704 2 : int bClipSrc = FALSE;
705 2 : OGRGeometry *poClipSrc = NULL;
706 2 : const char *pszClipSrcDS = NULL;
707 2 : const char *pszClipSrcSQL = NULL;
708 2 : const char *pszClipSrcLayer = NULL;
709 2 : const char *pszClipSrcWhere = NULL;
710 :
711 : /* Check strict compilation and runtime library version as we use C++ API */
712 2 : if (! GDAL_CHECK_VERSION(argv[0]))
713 0 : exit(1);
714 :
715 2 : GDALAllRegister();
716 2 : OGRRegisterAll();
717 :
718 2 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
719 2 : if( argc < 1 )
720 0 : exit( -argc );
721 :
722 : /* -------------------------------------------------------------------- */
723 : /* Parse arguments. */
724 : /* -------------------------------------------------------------------- */
725 13 : for( i = 1; i < argc; i++ )
726 : {
727 12 : if( EQUAL(argv[i], "--utility_version") )
728 : {
729 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
730 1 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
731 1 : return 0;
732 : }
733 11 : else if( EQUAL(argv[i],"-of") && i < argc-1 )
734 : {
735 0 : pszFormat = argv[++i];
736 : }
737 :
738 11 : else if( EQUAL(argv[i],"-q") || EQUAL(argv[i],"-quiet") )
739 : {
740 0 : bQuiet = TRUE;
741 0 : pfnProgress = GDALDummyProgress;
742 : }
743 :
744 12 : else if( EQUAL(argv[i],"-ot") && i < argc-1 )
745 : {
746 : int iType;
747 :
748 12 : for( iType = 1; iType < GDT_TypeCount; iType++ )
749 : {
750 22 : if( GDALGetDataTypeName((GDALDataType)iType) != NULL
751 11 : && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
752 : argv[i+1]) )
753 : {
754 1 : eOutputType = (GDALDataType) iType;
755 : }
756 : }
757 :
758 1 : if( eOutputType == GDT_Unknown )
759 : {
760 : fprintf( stderr, "FAILURE: Unknown output pixel type: %s\n",
761 0 : argv[i + 1] );
762 0 : Usage();
763 : }
764 1 : i++;
765 : }
766 :
767 11 : else if( EQUAL(argv[i],"-txe") && i < argc-2 )
768 : {
769 1 : dfXMin = atof(argv[++i]);
770 1 : dfXMax = atof(argv[++i]);
771 1 : bIsXExtentSet = TRUE;
772 : }
773 :
774 10 : else if( EQUAL(argv[i],"-tye") && i < argc-2 )
775 : {
776 1 : dfYMin = atof(argv[++i]);
777 1 : dfYMax = atof(argv[++i]);
778 1 : bIsYExtentSet = TRUE;
779 : }
780 :
781 9 : else if( EQUAL(argv[i],"-outsize") && i < argc-2 )
782 : {
783 1 : nXSize = atoi(argv[++i]);
784 1 : nYSize = atoi(argv[++i]);
785 : }
786 :
787 10 : else if( EQUAL(argv[i],"-co") && i < argc-1 )
788 : {
789 3 : papszCreateOptions = CSLAddString( papszCreateOptions, argv[++i] );
790 : }
791 :
792 4 : else if( EQUAL(argv[i],"-zfield") && i < argc-1 )
793 : {
794 0 : pszBurnAttribute = argv[++i];
795 : }
796 :
797 4 : else if( EQUAL(argv[i],"-where") && i < argc-1 )
798 : {
799 0 : pszWHERE = argv[++i];
800 : }
801 :
802 5 : else if( EQUAL(argv[i],"-l") && i < argc-1 )
803 : {
804 1 : papszLayers = CSLAddString( papszLayers, argv[++i] );
805 : }
806 :
807 3 : else if( EQUAL(argv[i],"-sql") && i < argc-1 )
808 : {
809 0 : pszSQL = argv[++i];
810 : }
811 :
812 3 : else if( EQUAL(argv[i],"-spat")
813 0 : && argv[i+1] != NULL
814 0 : && argv[i+2] != NULL
815 0 : && argv[i+3] != NULL
816 0 : && argv[i+4] != NULL )
817 : {
818 0 : OGRLinearRing oRing;
819 :
820 0 : oRing.addPoint( atof(argv[i+1]), atof(argv[i+2]) );
821 0 : oRing.addPoint( atof(argv[i+1]), atof(argv[i+4]) );
822 0 : oRing.addPoint( atof(argv[i+3]), atof(argv[i+4]) );
823 0 : oRing.addPoint( atof(argv[i+3]), atof(argv[i+2]) );
824 0 : oRing.addPoint( atof(argv[i+1]), atof(argv[i+2]) );
825 :
826 0 : poSpatialFilter = new OGRPolygon();
827 0 : ((OGRPolygon *) poSpatialFilter)->addRing( &oRing );
828 0 : i += 4;
829 : }
830 :
831 3 : else if ( EQUAL(argv[i],"-clipsrc") && i < argc - 1 )
832 : {
833 0 : bClipSrc = TRUE;
834 0 : errno = 0;
835 0 : const double unused = strtod( argv[i + 1], NULL ); // XXX: is it a number or not?
836 0 : if ( errno != 0
837 0 : && argv[i + 2] != NULL
838 0 : && argv[i + 3] != NULL
839 0 : && argv[i + 4] != NULL)
840 : {
841 0 : OGRLinearRing oRing;
842 :
843 0 : oRing.addPoint( atof(argv[i + 1]), atof(argv[i + 2]) );
844 0 : oRing.addPoint( atof(argv[i + 1]), atof(argv[i + 4]) );
845 0 : oRing.addPoint( atof(argv[i + 3]), atof(argv[i + 4]) );
846 0 : oRing.addPoint( atof(argv[i + 3]), atof(argv[i + 2]) );
847 0 : oRing.addPoint( atof(argv[i + 1]), atof(argv[i + 2]) );
848 :
849 0 : poClipSrc = new OGRPolygon();
850 0 : ((OGRPolygon *) poClipSrc)->addRing( &oRing );
851 0 : i += 4;
852 :
853 0 : (void)unused;
854 : }
855 0 : else if (EQUALN(argv[i + 1], "POLYGON", 7)
856 0 : || EQUALN(argv[i + 1], "MULTIPOLYGON", 12))
857 : {
858 0 : OGRGeometryFactory::createFromWkt(&argv[i + 1], NULL, &poClipSrc);
859 0 : if ( poClipSrc == NULL )
860 : {
861 : fprintf( stderr, "FAILURE: Invalid geometry. "
862 0 : "Must be a valid POLYGON or MULTIPOLYGON WKT\n\n");
863 0 : Usage();
864 : }
865 0 : i++;
866 : }
867 0 : else if (EQUAL(argv[i + 1], "spat_extent") )
868 : {
869 0 : i++;
870 : }
871 : else
872 : {
873 0 : pszClipSrcDS = argv[i + 1];
874 0 : i++;
875 : }
876 : }
877 :
878 3 : else if ( EQUAL(argv[i], "-clipsrcsql") && i < argc - 1 )
879 : {
880 0 : pszClipSrcSQL = argv[i + 1];
881 0 : i++;
882 : }
883 :
884 3 : else if ( EQUAL(argv[i], "-clipsrclayer") && i < argc - 1 )
885 : {
886 0 : pszClipSrcLayer = argv[i + 1];
887 0 : i++;
888 : }
889 :
890 3 : else if ( EQUAL(argv[i], "-clipsrcwhere") && i < argc - 1 )
891 : {
892 0 : pszClipSrcWhere = argv[i + 1];
893 0 : i++;
894 : }
895 :
896 3 : else if( EQUAL(argv[i],"-a_srs") && i < argc-1 )
897 : {
898 0 : OGRSpatialReference oOutputSRS;
899 :
900 0 : if( oOutputSRS.SetFromUserInput( argv[i+1] ) != OGRERR_NONE )
901 : {
902 : fprintf( stderr, "Failed to process SRS definition: %s\n",
903 0 : argv[i+1] );
904 0 : GDALDestroyDriverManager();
905 0 : exit( 1 );
906 : }
907 :
908 0 : oOutputSRS.exportToWkt( &pszOutputSRS );
909 0 : i++;
910 : }
911 :
912 4 : else if( EQUAL(argv[i],"-a") && i < argc-1 )
913 : {
914 1 : if ( ParseAlgorithmAndOptions( argv[++i], &eAlgorithm, &pOptions )
915 : != CE_None )
916 : {
917 : fprintf( stderr,
918 0 : "Failed to process algoritm name and parameters.\n" );
919 0 : exit( 1 );
920 : }
921 : }
922 :
923 2 : else if( argv[i][0] == '-' )
924 : {
925 : fprintf( stderr,
926 : "FAILURE: Option %s incomplete, or not recognised.\n\n",
927 0 : argv[i] );
928 0 : Usage();
929 : }
930 :
931 2 : else if( pszSource == NULL )
932 : {
933 1 : pszSource = argv[i];
934 : }
935 :
936 1 : else if( pszDest == NULL )
937 : {
938 1 : pszDest = argv[i];
939 : }
940 :
941 : else
942 : {
943 0 : fprintf( stderr, "FAILURE: Too many command options.\n\n" );
944 0 : Usage();
945 : }
946 : }
947 :
948 1 : if( pszSource == NULL || pszDest == NULL
949 : || (pszSQL == NULL && papszLayers == NULL) )
950 : {
951 0 : Usage();
952 : }
953 :
954 1 : if ( bClipSrc && pszClipSrcDS != NULL )
955 : {
956 : poClipSrc = LoadGeometry( pszClipSrcDS, pszClipSrcSQL,
957 0 : pszClipSrcLayer, pszClipSrcWhere );
958 0 : if ( poClipSrc == NULL )
959 : {
960 0 : fprintf( stderr, "FAILURE: cannot load source clip geometry\n\n" );
961 0 : Usage();
962 : }
963 : }
964 1 : else if ( bClipSrc && poClipSrc == NULL )
965 : {
966 0 : if ( poSpatialFilter )
967 0 : poClipSrc = poSpatialFilter->clone();
968 0 : if ( poClipSrc == NULL )
969 : {
970 : fprintf( stderr,
971 : "FAILURE: -clipsrc must be used with -spat option or \n"
972 : "a bounding box, WKT string or datasource must be "
973 0 : "specified\n\n" );
974 0 : Usage();
975 : }
976 : }
977 :
978 : /* -------------------------------------------------------------------- */
979 : /* Find the output driver. */
980 : /* -------------------------------------------------------------------- */
981 1 : hDriver = GDALGetDriverByName( pszFormat );
982 1 : if( hDriver == NULL )
983 : {
984 : int iDr;
985 :
986 : fprintf( stderr,
987 0 : "FAILURE: Output driver `%s' not recognised.\n", pszFormat );
988 : fprintf( stderr,
989 0 : "The following format drivers are configured and support output:\n" );
990 0 : for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
991 : {
992 0 : GDALDriverH hDriver = GDALGetDriver(iDr);
993 :
994 0 : if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) != NULL
995 : || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY,
996 : NULL ) != NULL )
997 : {
998 : fprintf( stderr, " %s: %s\n",
999 : GDALGetDriverShortName( hDriver ),
1000 0 : GDALGetDriverLongName( hDriver ) );
1001 : }
1002 : }
1003 0 : printf( "\n" );
1004 0 : Usage();
1005 : }
1006 :
1007 : /* -------------------------------------------------------------------- */
1008 : /* Open input datasource. */
1009 : /* -------------------------------------------------------------------- */
1010 : OGRDataSourceH hSrcDS;
1011 :
1012 1 : hSrcDS = OGROpen( pszSource, FALSE, NULL );
1013 1 : if( hSrcDS == NULL )
1014 : {
1015 : fprintf( stderr, "Unable to open input datasource \"%s\".\n",
1016 0 : pszSource );
1017 0 : fprintf( stderr, "%s\n", CPLGetLastErrorMsg() );
1018 0 : exit( 3 );
1019 : }
1020 :
1021 : /* -------------------------------------------------------------------- */
1022 : /* Create target raster file. */
1023 : /* -------------------------------------------------------------------- */
1024 : GDALDatasetH hDstDS;
1025 1 : int nLayerCount = CSLCount(papszLayers);
1026 1 : int nBands = nLayerCount;
1027 :
1028 1 : if ( pszSQL )
1029 0 : nBands++;
1030 :
1031 : // FIXME
1032 1 : if ( nXSize == 0 )
1033 0 : nXSize = 256;
1034 1 : if ( nYSize == 0 )
1035 0 : nYSize = 256;
1036 :
1037 : hDstDS = GDALCreate( hDriver, pszDest, nXSize, nYSize, nBands,
1038 1 : eOutputType, papszCreateOptions );
1039 1 : if ( hDstDS == NULL )
1040 : {
1041 : fprintf( stderr, "Unable to create target dataset \"%s\".\n",
1042 0 : pszDest );
1043 0 : fprintf( stderr, "%s\n", CPLGetLastErrorMsg() );
1044 0 : exit( 3 );
1045 : }
1046 :
1047 : /* -------------------------------------------------------------------- */
1048 : /* If algorithm was not specified assigh default one. */
1049 : /* -------------------------------------------------------------------- */
1050 1 : if ( !pOptions )
1051 0 : ParseAlgorithmAndOptions( szAlgNameInvDist, &eAlgorithm, &pOptions );
1052 :
1053 : /* -------------------------------------------------------------------- */
1054 : /* Process SQL request. */
1055 : /* -------------------------------------------------------------------- */
1056 1 : if( pszSQL != NULL )
1057 : {
1058 : OGRLayerH hLayer;
1059 :
1060 : hLayer = OGR_DS_ExecuteSQL( hSrcDS, pszSQL,
1061 0 : (OGRGeometryH)poSpatialFilter, NULL );
1062 0 : if( hLayer != NULL )
1063 : {
1064 : // Custom layer will be rasterized in the first band.
1065 : ProcessLayer( hLayer, hDstDS, poClipSrc, nXSize, nYSize, 1,
1066 : bIsXExtentSet, bIsYExtentSet,
1067 : dfXMin, dfXMax, dfYMin, dfYMax, pszBurnAttribute,
1068 : eOutputType, eAlgorithm, pOptions,
1069 0 : bQuiet, pfnProgress );
1070 : }
1071 : }
1072 :
1073 : /* -------------------------------------------------------------------- */
1074 : /* Process each layer. */
1075 : /* -------------------------------------------------------------------- */
1076 2 : for( i = 0; i < nLayerCount; i++ )
1077 : {
1078 1 : OGRLayerH hLayer = OGR_DS_GetLayerByName( hSrcDS, papszLayers[i] );
1079 1 : if( hLayer == NULL )
1080 : {
1081 : fprintf( stderr, "Unable to find layer \"%s\", skipping.\n",
1082 0 : papszLayers[i] );
1083 0 : continue;
1084 : }
1085 :
1086 1 : if( pszWHERE )
1087 : {
1088 0 : if( OGR_L_SetAttributeFilter( hLayer, pszWHERE ) != OGRERR_NONE )
1089 0 : break;
1090 : }
1091 :
1092 1 : if ( poSpatialFilter != NULL )
1093 0 : OGR_L_SetSpatialFilter( hLayer, (OGRGeometryH)poSpatialFilter );
1094 :
1095 : // Fetch the first meaningful SRS definition
1096 1 : if ( !pszOutputSRS )
1097 : {
1098 1 : OGRSpatialReferenceH hSRS = OGR_L_GetSpatialRef( hLayer );
1099 1 : if ( hSRS )
1100 0 : OSRExportToWkt( hSRS, &pszOutputSRS );
1101 : }
1102 :
1103 : ProcessLayer( hLayer, hDstDS, poClipSrc, nXSize, nYSize,
1104 : i + 1 + nBands - nLayerCount,
1105 : bIsXExtentSet, bIsYExtentSet,
1106 : dfXMin, dfXMax, dfYMin, dfYMax, pszBurnAttribute,
1107 : eOutputType, eAlgorithm, pOptions,
1108 1 : bQuiet, pfnProgress );
1109 : }
1110 :
1111 : /* -------------------------------------------------------------------- */
1112 : /* Apply geotransformation matrix. */
1113 : /* -------------------------------------------------------------------- */
1114 : double adfGeoTransform[6];
1115 1 : adfGeoTransform[0] = dfXMin;
1116 1 : adfGeoTransform[1] = (dfXMax - dfXMin) / nXSize;
1117 1 : adfGeoTransform[2] = 0.0;
1118 1 : adfGeoTransform[3] = dfYMin;
1119 1 : adfGeoTransform[4] = 0.0;
1120 1 : adfGeoTransform[5] = (dfYMax - dfYMin) / nYSize;
1121 1 : GDALSetGeoTransform( hDstDS, adfGeoTransform );
1122 :
1123 : /* -------------------------------------------------------------------- */
1124 : /* Apply SRS definition if set. */
1125 : /* -------------------------------------------------------------------- */
1126 1 : if ( pszOutputSRS )
1127 : {
1128 0 : GDALSetProjection( hDstDS, pszOutputSRS );
1129 0 : CPLFree( pszOutputSRS );
1130 : }
1131 :
1132 : /* -------------------------------------------------------------------- */
1133 : /* Cleanup */
1134 : /* -------------------------------------------------------------------- */
1135 1 : OGR_DS_Destroy( hSrcDS );
1136 1 : GDALClose( hDstDS );
1137 1 : OGRGeometryFactory::destroyGeometry( poSpatialFilter );
1138 1 : OGRGeometryFactory::destroyGeometry( poClipSrc );
1139 :
1140 1 : CPLFree( pOptions );
1141 1 : CSLDestroy( papszCreateOptions );
1142 1 : CSLDestroy( argv );
1143 1 : CSLDestroy( papszLayers );
1144 :
1145 1 : OGRCleanupAll();
1146 :
1147 1 : GDALDestroyDriverManager();
1148 :
1149 1 : return 0;
1150 : }
1151 :
|
