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