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