1 : /* ****************************************************************************
2 : * $Id: gdal_grid.cpp 25174 2012-10-21 21:29:42Z 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 25174 2012-10-21 21:29:42Z 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 27 : static void PrintAlgorithmAndOptions( GDALGridAlgorithm eAlgorithm,
97 : void *pOptions )
98 : {
99 27 : switch ( eAlgorithm )
100 : {
101 : case GGA_InverseDistanceToAPower:
102 5 : 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 5 : ((GDALGridInverseDistanceToAPowerOptions *)pOptions)->dfNoDataValue);
114 5 : 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 27 : }
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 25041 : 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 25041 : if ( poClipSrc && !poGeom->Within(poClipSrc) )
223 0 : return;
224 :
225 25041 : adfX.push_back( poGeom->getX() );
226 25041 : adfY.push_back( poGeom->getY() );
227 25041 : if ( iBurnField < 0 )
228 25041 : 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 27 : static CPLErr 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 27 : int iBurnField = -1;
256 :
257 27 : 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 CE_Failure;
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 27 : std::vector<double> adfX, adfY, adfZ;
276 :
277 27 : OGR_L_ResetReading( hSrcLayer );
278 :
279 25095 : while( (poFeat = (OGRFeature *)OGR_L_GetNextFeature( hSrcLayer )) != NULL )
280 : {
281 25041 : OGRGeometry *poGeom = poFeat->GetGeometryRef();
282 :
283 25041 : if ( poGeom != NULL )
284 : {
285 25041 : OGRwkbGeometryType eType = wkbFlatten( poGeom->getGeometryType() );
286 25041 : double dfBurnValue = 0.0;
287 :
288 25041 : if ( iBurnField >= 0 )
289 0 : dfBurnValue = poFeat->GetFieldAsDouble( iBurnField );
290 :
291 25041 : 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 25041 : else if ( eType == wkbPoint )
304 : {
305 : ProcessGeometry( (OGRPoint *)poGeom, poClipSrc,
306 25041 : iBurnField, dfBurnValue, adfX, adfY, adfZ);
307 : }
308 : }
309 :
310 25041 : OGRFeature::DestroyFeature( poFeat );
311 : }
312 :
313 27 : 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 0 : return CE_None;
318 : }
319 :
320 : /* -------------------------------------------------------------------- */
321 : /* Compute grid geometry. */
322 : /* -------------------------------------------------------------------- */
323 27 : 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 27 : const double dfDeltaX = ( dfXMax - dfXMin ) / nXSize;
348 27 : const double dfDeltaY = ( dfYMax - dfYMin ) / nYSize;
349 :
350 27 : if ( !bQuiet )
351 : {
352 27 : printf( "Grid data type is \"%s\"\n", GDALGetDataTypeName(eType) );
353 : printf( "Grid size = (%lu %lu).\n",
354 27 : (unsigned long)nXSize, (unsigned long)nYSize );
355 : printf( "Corner coordinates = (%f %f)-(%f %f).\n",
356 : dfXMin - dfDeltaX / 2, dfYMax + dfDeltaY / 2,
357 27 : dfXMax + dfDeltaX / 2, dfYMin - dfDeltaY / 2 );
358 27 : printf( "Grid cell size = (%f %f).\n", dfDeltaX, dfDeltaY );
359 27 : printf( "Source point count = %lu.\n", (unsigned long)adfX.size() );
360 27 : PrintAlgorithmAndOptions( eAlgorithm, pOptions );
361 27 : printf("\n");
362 : }
363 :
364 27 : GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, nBand );
365 :
366 27 : if (adfX.size() == 0)
367 : {
368 : // FIXME: Shoulda' set to nodata value instead
369 0 : GDALFillRaster( hBand, 0.0 , 0.0 );
370 0 : return CE_None;
371 : }
372 :
373 : GUInt32 nXOffset, nYOffset;
374 : int nBlockXSize, nBlockYSize;
375 27 : int nDataTypeSize = GDALGetDataTypeSize(eType) / 8;
376 :
377 : // Try to grow the work buffer up to 16 MB if it is smaller
378 27 : GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );
379 27 : const GUInt32 nDesiredBufferSize = 16*1024*1024;
380 27 : if( (GUInt32)nBlockXSize < nXSize && (GUInt32)nBlockYSize < nYSize &&
381 : (GUInt32)nBlockXSize < nDesiredBufferSize / (nBlockYSize * nDataTypeSize) )
382 : {
383 0 : int nNewBlockXSize = nDesiredBufferSize / (nBlockYSize * nDataTypeSize);
384 0 : nBlockXSize = (nNewBlockXSize / nBlockXSize) * nBlockXSize;
385 0 : if( (GUInt32)nBlockXSize > nXSize )
386 0 : nBlockXSize = nXSize;
387 : }
388 27 : else if( (GUInt32)nBlockXSize == nXSize && (GUInt32)nBlockYSize < nYSize &&
389 : (GUInt32)nBlockYSize < nDesiredBufferSize / (nXSize * nDataTypeSize) )
390 : {
391 0 : int nNewBlockYSize = nDesiredBufferSize / (nXSize * nDataTypeSize);
392 0 : nBlockYSize = (nNewBlockYSize / nBlockYSize) * nBlockYSize;
393 0 : if( (GUInt32)nBlockYSize > nYSize )
394 0 : nBlockYSize = nYSize;
395 : }
396 27 : CPLDebug("GDAL_GRID", "Work buffer: %d * %d", nBlockXSize, nBlockYSize);
397 :
398 : void *pData =
399 27 : VSIMalloc3( nBlockXSize, nBlockYSize, nDataTypeSize );
400 27 : if( pData == NULL )
401 : {
402 0 : CPLError(CE_Failure, CPLE_OutOfMemory, "Cannot allocate work buffer");
403 0 : return CE_Failure;
404 : }
405 :
406 27 : GUInt32 nBlock = 0;
407 : GUInt32 nBlockCount = ((nXSize + nBlockXSize - 1) / nBlockXSize)
408 27 : * ((nYSize + nBlockYSize - 1) / nBlockYSize);
409 :
410 27 : CPLErr eErr = CE_None;
411 54 : for ( nYOffset = 0; nYOffset < nYSize && eErr == CE_None; nYOffset += nBlockYSize )
412 : {
413 54 : for ( nXOffset = 0; nXOffset < nXSize && eErr == CE_None; nXOffset += nBlockXSize )
414 : {
415 : void *pScaledProgress;
416 : pScaledProgress =
417 : GDALCreateScaledProgress( (double)nBlock / nBlockCount,
418 : (double)(nBlock + 1) / nBlockCount,
419 27 : pfnProgress, NULL );
420 27 : nBlock ++;
421 :
422 27 : int nXRequest = nBlockXSize;
423 27 : if (nXOffset + nXRequest > nXSize)
424 1 : nXRequest = nXSize - nXOffset;
425 :
426 27 : int nYRequest = nBlockYSize;
427 27 : if (nYOffset + nYRequest > nYSize)
428 1 : nYRequest = nYSize - nYOffset;
429 :
430 : eErr = GDALGridCreate( eAlgorithm, pOptions,
431 : adfX.size(), &(adfX[0]), &(adfY[0]), &(adfZ[0]),
432 : dfXMin + dfDeltaX * nXOffset,
433 : dfXMin + dfDeltaX * (nXOffset + nXRequest),
434 : dfYMin + dfDeltaY * nYOffset,
435 : dfYMin + dfDeltaY * (nYOffset + nYRequest),
436 : nXRequest, nYRequest, eType, pData,
437 27 : GDALScaledProgress, pScaledProgress );
438 :
439 27 : if( eErr == CE_None )
440 : eErr = GDALRasterIO( hBand, GF_Write, nXOffset, nYOffset,
441 : nXRequest, nYRequest, pData,
442 27 : nXRequest, nYRequest, eType, 0, 0 );
443 :
444 27 : GDALDestroyScaledProgress( pScaledProgress );
445 : }
446 : }
447 :
448 27 : CPLFree( pData );
449 27 : return eErr;
450 : }
451 :
452 : /************************************************************************/
453 : /* LoadGeometry() */
454 : /* */
455 : /* Read geometries from the given dataset using specified filters and */
456 : /* returns a collection of read geometries. */
457 : /************************************************************************/
458 :
459 0 : static OGRGeometryCollection* LoadGeometry( const char* pszDS,
460 : const char* pszSQL,
461 : const char* pszLyr,
462 : const char* pszWhere )
463 : {
464 : OGRDataSource *poDS;
465 : OGRLayer *poLyr;
466 : OGRFeature *poFeat;
467 0 : OGRGeometryCollection *poGeom = NULL;
468 :
469 0 : poDS = OGRSFDriverRegistrar::Open( pszDS, FALSE );
470 0 : if ( poDS == NULL )
471 0 : return NULL;
472 :
473 0 : if ( pszSQL != NULL )
474 0 : poLyr = poDS->ExecuteSQL( pszSQL, NULL, NULL );
475 0 : else if ( pszLyr != NULL )
476 0 : poLyr = poDS->GetLayerByName( pszLyr );
477 : else
478 0 : poLyr = poDS->GetLayer(0);
479 :
480 0 : if ( poLyr == NULL )
481 : {
482 : fprintf( stderr,
483 0 : "FAILURE: Failed to identify source layer from datasource.\n" );
484 0 : OGRDataSource::DestroyDataSource( poDS );
485 0 : return NULL;
486 : }
487 :
488 0 : if ( pszWhere )
489 0 : poLyr->SetAttributeFilter( pszWhere );
490 :
491 0 : while ( (poFeat = poLyr->GetNextFeature()) != NULL )
492 : {
493 0 : OGRGeometry* poSrcGeom = poFeat->GetGeometryRef();
494 0 : if ( poSrcGeom )
495 : {
496 : OGRwkbGeometryType eType =
497 0 : wkbFlatten( poSrcGeom->getGeometryType() );
498 :
499 0 : if ( poGeom == NULL )
500 0 : poGeom = new OGRMultiPolygon();
501 :
502 0 : if ( eType == wkbPolygon )
503 0 : poGeom->addGeometry( poSrcGeom );
504 0 : else if ( eType == wkbMultiPolygon )
505 : {
506 : int iGeom;
507 : int nGeomCount =
508 0 : ((OGRMultiPolygon *)poSrcGeom)->getNumGeometries();
509 :
510 0 : for ( iGeom = 0; iGeom < nGeomCount; iGeom++ )
511 : {
512 : poGeom->addGeometry(
513 0 : ((OGRMultiPolygon *)poSrcGeom)->getGeometryRef(iGeom) );
514 : }
515 : }
516 : else
517 : {
518 0 : fprintf( stderr, "FAILURE: Geometry not of polygon type.\n" );
519 0 : OGRGeometryFactory::destroyGeometry( poGeom );
520 0 : OGRFeature::DestroyFeature( poFeat );
521 0 : if ( pszSQL != NULL )
522 0 : poDS->ReleaseResultSet( poLyr );
523 0 : OGRDataSource::DestroyDataSource( poDS );
524 0 : return NULL;
525 : }
526 : }
527 :
528 0 : OGRFeature::DestroyFeature( poFeat );
529 : }
530 :
531 0 : if( pszSQL != NULL )
532 0 : poDS->ReleaseResultSet( poLyr );
533 0 : OGRDataSource::DestroyDataSource( poDS );
534 :
535 0 : return poGeom;
536 : }
537 :
538 : /************************************************************************/
539 : /* main() */
540 : /************************************************************************/
541 :
542 28 : int main( int argc, char ** argv )
543 : {
544 : GDALDriverH hDriver;
545 28 : const char *pszSource=NULL, *pszDest=NULL, *pszFormat = "GTiff";
546 28 : int bFormatExplicitelySet = FALSE;
547 28 : char **papszLayers = NULL;
548 28 : const char *pszBurnAttribute = NULL;
549 28 : const char *pszWHERE = NULL, *pszSQL = NULL;
550 28 : GDALDataType eOutputType = GDT_Float64;
551 28 : char **papszCreateOptions = NULL;
552 28 : GUInt32 nXSize = 0, nYSize = 0;
553 28 : double dfXMin = 0.0, dfXMax = 0.0, dfYMin = 0.0, dfYMax = 0.0;
554 28 : int bIsXExtentSet = FALSE, bIsYExtentSet = FALSE;
555 28 : GDALGridAlgorithm eAlgorithm = GGA_InverseDistanceToAPower;
556 28 : void *pOptions = NULL;
557 28 : char *pszOutputSRS = NULL;
558 28 : int bQuiet = FALSE;
559 28 : GDALProgressFunc pfnProgress = GDALTermProgress;
560 : int i;
561 28 : OGRGeometry *poSpatialFilter = NULL;
562 28 : int bClipSrc = FALSE;
563 28 : OGRGeometry *poClipSrc = NULL;
564 28 : const char *pszClipSrcDS = NULL;
565 28 : const char *pszClipSrcSQL = NULL;
566 28 : const char *pszClipSrcLayer = NULL;
567 28 : const char *pszClipSrcWhere = NULL;
568 :
569 : /* Check strict compilation and runtime library version as we use C++ API */
570 28 : if (! GDAL_CHECK_VERSION(argv[0]))
571 0 : exit(1);
572 :
573 28 : GDALAllRegister();
574 28 : OGRRegisterAll();
575 :
576 28 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
577 28 : if( argc < 1 )
578 0 : exit( -argc );
579 :
580 : /* -------------------------------------------------------------------- */
581 : /* Parse arguments. */
582 : /* -------------------------------------------------------------------- */
583 247 : for( i = 1; i < argc; i++ )
584 : {
585 220 : if( EQUAL(argv[i], "--utility_version") )
586 : {
587 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
588 1 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
589 1 : return 0;
590 : }
591 219 : else if( EQUAL(argv[i],"-of") && i < argc-1 )
592 : {
593 0 : pszFormat = argv[++i];
594 0 : bFormatExplicitelySet = TRUE;
595 : }
596 :
597 219 : else if( EQUAL(argv[i],"-q") || EQUAL(argv[i],"-quiet") )
598 : {
599 0 : bQuiet = TRUE;
600 0 : pfnProgress = GDALDummyProgress;
601 : }
602 :
603 246 : else if( EQUAL(argv[i],"-ot") && i < argc-1 )
604 : {
605 : int iType;
606 :
607 324 : for( iType = 1; iType < GDT_TypeCount; iType++ )
608 : {
609 594 : if( GDALGetDataTypeName((GDALDataType)iType) != NULL
610 297 : && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
611 : argv[i+1]) )
612 : {
613 27 : eOutputType = (GDALDataType) iType;
614 : }
615 : }
616 :
617 27 : if( eOutputType == GDT_Unknown )
618 : {
619 : fprintf( stderr, "FAILURE: Unknown output pixel type: %s\n",
620 0 : argv[i + 1] );
621 0 : Usage();
622 : }
623 27 : i++;
624 : }
625 :
626 219 : else if( EQUAL(argv[i],"-txe") && i < argc-2 )
627 : {
628 27 : dfXMin = atof(argv[++i]);
629 27 : dfXMax = atof(argv[++i]);
630 27 : bIsXExtentSet = TRUE;
631 : }
632 :
633 192 : else if( EQUAL(argv[i],"-tye") && i < argc-2 )
634 : {
635 27 : dfYMin = atof(argv[++i]);
636 27 : dfYMax = atof(argv[++i]);
637 27 : bIsYExtentSet = TRUE;
638 : }
639 :
640 165 : else if( EQUAL(argv[i],"-outsize") && i < argc-2 )
641 : {
642 27 : nXSize = atoi(argv[++i]);
643 27 : nYSize = atoi(argv[++i]);
644 : }
645 :
646 114 : else if( EQUAL(argv[i],"-co") && i < argc-1 )
647 : {
648 3 : papszCreateOptions = CSLAddString( papszCreateOptions, argv[++i] );
649 : }
650 :
651 108 : else if( EQUAL(argv[i],"-zfield") && i < argc-1 )
652 : {
653 0 : pszBurnAttribute = argv[++i];
654 : }
655 :
656 108 : else if( EQUAL(argv[i],"-where") && i < argc-1 )
657 : {
658 0 : pszWHERE = argv[++i];
659 : }
660 :
661 135 : else if( EQUAL(argv[i],"-l") && i < argc-1 )
662 : {
663 27 : papszLayers = CSLAddString( papszLayers, argv[++i] );
664 : }
665 :
666 81 : else if( EQUAL(argv[i],"-sql") && i < argc-1 )
667 : {
668 0 : pszSQL = argv[++i];
669 : }
670 :
671 81 : else if( EQUAL(argv[i],"-spat")
672 0 : && argv[i+1] != NULL
673 0 : && argv[i+2] != NULL
674 0 : && argv[i+3] != NULL
675 0 : && argv[i+4] != NULL )
676 : {
677 0 : OGRLinearRing oRing;
678 :
679 0 : oRing.addPoint( atof(argv[i+1]), atof(argv[i+2]) );
680 0 : oRing.addPoint( atof(argv[i+1]), atof(argv[i+4]) );
681 0 : oRing.addPoint( atof(argv[i+3]), atof(argv[i+4]) );
682 0 : oRing.addPoint( atof(argv[i+3]), atof(argv[i+2]) );
683 0 : oRing.addPoint( atof(argv[i+1]), atof(argv[i+2]) );
684 :
685 0 : poSpatialFilter = new OGRPolygon();
686 0 : ((OGRPolygon *) poSpatialFilter)->addRing( &oRing );
687 0 : i += 4;
688 : }
689 :
690 81 : else if ( EQUAL(argv[i],"-clipsrc") && i < argc - 1 )
691 : {
692 0 : bClipSrc = TRUE;
693 0 : errno = 0;
694 0 : const double unused = strtod( argv[i + 1], NULL ); // XXX: is it a number or not?
695 0 : if ( errno != 0
696 0 : && argv[i + 2] != NULL
697 0 : && argv[i + 3] != NULL
698 0 : && argv[i + 4] != NULL)
699 : {
700 0 : OGRLinearRing oRing;
701 :
702 0 : oRing.addPoint( atof(argv[i + 1]), atof(argv[i + 2]) );
703 0 : oRing.addPoint( atof(argv[i + 1]), atof(argv[i + 4]) );
704 0 : oRing.addPoint( atof(argv[i + 3]), atof(argv[i + 4]) );
705 0 : oRing.addPoint( atof(argv[i + 3]), atof(argv[i + 2]) );
706 0 : oRing.addPoint( atof(argv[i + 1]), atof(argv[i + 2]) );
707 :
708 0 : poClipSrc = new OGRPolygon();
709 0 : ((OGRPolygon *) poClipSrc)->addRing( &oRing );
710 0 : i += 4;
711 :
712 0 : (void)unused;
713 : }
714 0 : else if (EQUALN(argv[i + 1], "POLYGON", 7)
715 0 : || EQUALN(argv[i + 1], "MULTIPOLYGON", 12))
716 : {
717 0 : OGRGeometryFactory::createFromWkt(&argv[i + 1], NULL, &poClipSrc);
718 0 : if ( poClipSrc == NULL )
719 : {
720 : fprintf( stderr, "FAILURE: Invalid geometry. "
721 0 : "Must be a valid POLYGON or MULTIPOLYGON WKT\n\n");
722 0 : Usage();
723 : }
724 0 : i++;
725 : }
726 0 : else if (EQUAL(argv[i + 1], "spat_extent") )
727 : {
728 0 : i++;
729 : }
730 : else
731 : {
732 0 : pszClipSrcDS = argv[i + 1];
733 0 : i++;
734 : }
735 : }
736 :
737 81 : else if ( EQUAL(argv[i], "-clipsrcsql") && i < argc - 1 )
738 : {
739 0 : pszClipSrcSQL = argv[i + 1];
740 0 : i++;
741 : }
742 :
743 81 : else if ( EQUAL(argv[i], "-clipsrclayer") && i < argc - 1 )
744 : {
745 0 : pszClipSrcLayer = argv[i + 1];
746 0 : i++;
747 : }
748 :
749 81 : else if ( EQUAL(argv[i], "-clipsrcwhere") && i < argc - 1 )
750 : {
751 0 : pszClipSrcWhere = argv[i + 1];
752 0 : i++;
753 : }
754 :
755 81 : else if( EQUAL(argv[i],"-a_srs") && i < argc-1 )
756 : {
757 0 : OGRSpatialReference oOutputSRS;
758 :
759 0 : if( oOutputSRS.SetFromUserInput( argv[i+1] ) != OGRERR_NONE )
760 : {
761 : fprintf( stderr, "Failed to process SRS definition: %s\n",
762 0 : argv[i+1] );
763 0 : GDALDestroyDriverManager();
764 0 : exit( 1 );
765 : }
766 :
767 0 : oOutputSRS.exportToWkt( &pszOutputSRS );
768 0 : i++;
769 : }
770 :
771 108 : else if( EQUAL(argv[i],"-a") && i < argc-1 )
772 : {
773 27 : if ( ParseAlgorithmAndOptions( argv[++i], &eAlgorithm, &pOptions )
774 : != CE_None )
775 : {
776 : fprintf( stderr,
777 0 : "Failed to process algoritm name and parameters.\n" );
778 0 : exit( 1 );
779 : }
780 : }
781 :
782 54 : else if( argv[i][0] == '-' )
783 : {
784 : fprintf( stderr,
785 : "FAILURE: Option %s incomplete, or not recognised.\n\n",
786 0 : argv[i] );
787 0 : Usage();
788 : }
789 :
790 54 : else if( pszSource == NULL )
791 : {
792 27 : pszSource = argv[i];
793 : }
794 :
795 27 : else if( pszDest == NULL )
796 : {
797 27 : pszDest = argv[i];
798 : }
799 :
800 : else
801 : {
802 0 : fprintf( stderr, "FAILURE: Too many command options.\n\n" );
803 0 : Usage();
804 : }
805 : }
806 :
807 27 : if( pszSource == NULL || pszDest == NULL
808 : || (pszSQL == NULL && papszLayers == NULL) )
809 : {
810 0 : Usage();
811 : }
812 :
813 27 : if ( bClipSrc && pszClipSrcDS != NULL )
814 : {
815 : poClipSrc = LoadGeometry( pszClipSrcDS, pszClipSrcSQL,
816 0 : pszClipSrcLayer, pszClipSrcWhere );
817 0 : if ( poClipSrc == NULL )
818 : {
819 0 : fprintf( stderr, "FAILURE: cannot load source clip geometry\n\n" );
820 0 : Usage();
821 : }
822 : }
823 27 : else if ( bClipSrc && poClipSrc == NULL && !poSpatialFilter )
824 : {
825 : fprintf( stderr,
826 : "FAILURE: -clipsrc must be used with -spat option or \n"
827 : "a bounding box, WKT string or datasource must be "
828 0 : "specified\n\n" );
829 0 : Usage();
830 : }
831 :
832 27 : if ( poSpatialFilter )
833 : {
834 0 : if ( poClipSrc )
835 : {
836 0 : OGRGeometry *poTemp = poSpatialFilter->Intersection( poClipSrc );
837 :
838 0 : if ( poTemp )
839 : {
840 0 : OGRGeometryFactory::destroyGeometry( poSpatialFilter );
841 0 : poSpatialFilter = poTemp;
842 : }
843 :
844 0 : OGRGeometryFactory::destroyGeometry( poClipSrc );
845 0 : poClipSrc = NULL;
846 : }
847 : }
848 : else
849 : {
850 27 : if ( poClipSrc )
851 : {
852 0 : poSpatialFilter = poClipSrc;
853 0 : poClipSrc = NULL;
854 : }
855 : }
856 :
857 : /* -------------------------------------------------------------------- */
858 : /* Find the output driver. */
859 : /* -------------------------------------------------------------------- */
860 27 : hDriver = GDALGetDriverByName( pszFormat );
861 27 : if( hDriver == NULL )
862 : {
863 : int iDr;
864 :
865 : fprintf( stderr,
866 0 : "FAILURE: Output driver `%s' not recognised.\n", pszFormat );
867 : fprintf( stderr,
868 0 : "The following format drivers are configured and support output:\n" );
869 0 : for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
870 : {
871 0 : GDALDriverH hDriver = GDALGetDriver(iDr);
872 :
873 0 : if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) != NULL
874 : || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY,
875 : NULL ) != NULL )
876 : {
877 : fprintf( stderr, " %s: %s\n",
878 : GDALGetDriverShortName( hDriver ),
879 0 : GDALGetDriverLongName( hDriver ) );
880 : }
881 : }
882 0 : printf( "\n" );
883 0 : Usage();
884 : }
885 :
886 : /* -------------------------------------------------------------------- */
887 : /* Open input datasource. */
888 : /* -------------------------------------------------------------------- */
889 : OGRDataSourceH hSrcDS;
890 :
891 27 : hSrcDS = OGROpen( pszSource, FALSE, NULL );
892 27 : if( hSrcDS == NULL )
893 : {
894 : fprintf( stderr, "Unable to open input datasource \"%s\".\n",
895 0 : pszSource );
896 0 : fprintf( stderr, "%s\n", CPLGetLastErrorMsg() );
897 0 : exit( 3 );
898 : }
899 :
900 : /* -------------------------------------------------------------------- */
901 : /* Create target raster file. */
902 : /* -------------------------------------------------------------------- */
903 : GDALDatasetH hDstDS;
904 27 : int nLayerCount = CSLCount(papszLayers);
905 27 : int nBands = nLayerCount;
906 :
907 27 : if ( pszSQL )
908 0 : nBands++;
909 :
910 : // FIXME
911 27 : if ( nXSize == 0 )
912 0 : nXSize = 256;
913 27 : if ( nYSize == 0 )
914 0 : nYSize = 256;
915 :
916 27 : if (!bQuiet && !bFormatExplicitelySet)
917 27 : CheckExtensionConsistency(pszDest, pszFormat);
918 :
919 : hDstDS = GDALCreate( hDriver, pszDest, nXSize, nYSize, nBands,
920 27 : eOutputType, papszCreateOptions );
921 27 : if ( hDstDS == NULL )
922 : {
923 : fprintf( stderr, "Unable to create target dataset \"%s\".\n",
924 0 : pszDest );
925 0 : fprintf( stderr, "%s\n", CPLGetLastErrorMsg() );
926 0 : exit( 3 );
927 : }
928 :
929 : /* -------------------------------------------------------------------- */
930 : /* If algorithm was not specified assigh default one. */
931 : /* -------------------------------------------------------------------- */
932 27 : if ( !pOptions )
933 0 : ParseAlgorithmAndOptions( szAlgNameInvDist, &eAlgorithm, &pOptions );
934 :
935 : /* -------------------------------------------------------------------- */
936 : /* Process SQL request. */
937 : /* -------------------------------------------------------------------- */
938 27 : if( pszSQL != NULL )
939 : {
940 : OGRLayerH hLayer;
941 :
942 : hLayer = OGR_DS_ExecuteSQL( hSrcDS, pszSQL,
943 0 : (OGRGeometryH)poSpatialFilter, NULL );
944 0 : if( hLayer != NULL )
945 : {
946 : // Custom layer will be rasterized in the first band.
947 : ProcessLayer( hLayer, hDstDS, poSpatialFilter, nXSize, nYSize, 1,
948 : bIsXExtentSet, bIsYExtentSet,
949 : dfXMin, dfXMax, dfYMin, dfYMax, pszBurnAttribute,
950 : eOutputType, eAlgorithm, pOptions,
951 0 : bQuiet, pfnProgress );
952 : }
953 : }
954 :
955 : /* -------------------------------------------------------------------- */
956 : /* Process each layer. */
957 : /* -------------------------------------------------------------------- */
958 54 : for( i = 0; i < nLayerCount; i++ )
959 : {
960 27 : OGRLayerH hLayer = OGR_DS_GetLayerByName( hSrcDS, papszLayers[i] );
961 27 : if( hLayer == NULL )
962 : {
963 : fprintf( stderr, "Unable to find layer \"%s\", skipping.\n",
964 0 : papszLayers[i] );
965 0 : continue;
966 : }
967 :
968 27 : if( pszWHERE )
969 : {
970 0 : if( OGR_L_SetAttributeFilter( hLayer, pszWHERE ) != OGRERR_NONE )
971 0 : break;
972 : }
973 :
974 27 : if ( poSpatialFilter != NULL )
975 0 : OGR_L_SetSpatialFilter( hLayer, (OGRGeometryH)poSpatialFilter );
976 :
977 : // Fetch the first meaningful SRS definition
978 27 : if ( !pszOutputSRS )
979 : {
980 27 : OGRSpatialReferenceH hSRS = OGR_L_GetSpatialRef( hLayer );
981 27 : if ( hSRS )
982 26 : OSRExportToWkt( hSRS, &pszOutputSRS );
983 : }
984 :
985 : ProcessLayer( hLayer, hDstDS, poSpatialFilter, nXSize, nYSize,
986 : i + 1 + nBands - nLayerCount,
987 : bIsXExtentSet, bIsYExtentSet,
988 : dfXMin, dfXMax, dfYMin, dfYMax, pszBurnAttribute,
989 : eOutputType, eAlgorithm, pOptions,
990 27 : bQuiet, pfnProgress );
991 : }
992 :
993 : /* -------------------------------------------------------------------- */
994 : /* Apply geotransformation matrix. */
995 : /* -------------------------------------------------------------------- */
996 : double adfGeoTransform[6];
997 27 : adfGeoTransform[0] = dfXMin;
998 27 : adfGeoTransform[1] = (dfXMax - dfXMin) / nXSize;
999 27 : adfGeoTransform[2] = 0.0;
1000 27 : adfGeoTransform[3] = dfYMin;
1001 27 : adfGeoTransform[4] = 0.0;
1002 27 : adfGeoTransform[5] = (dfYMax - dfYMin) / nYSize;
1003 27 : GDALSetGeoTransform( hDstDS, adfGeoTransform );
1004 :
1005 : /* -------------------------------------------------------------------- */
1006 : /* Apply SRS definition if set. */
1007 : /* -------------------------------------------------------------------- */
1008 27 : if ( pszOutputSRS )
1009 : {
1010 26 : GDALSetProjection( hDstDS, pszOutputSRS );
1011 26 : CPLFree( pszOutputSRS );
1012 : }
1013 :
1014 : /* -------------------------------------------------------------------- */
1015 : /* Cleanup */
1016 : /* -------------------------------------------------------------------- */
1017 27 : OGR_DS_Destroy( hSrcDS );
1018 27 : GDALClose( hDstDS );
1019 27 : OGRGeometryFactory::destroyGeometry( poSpatialFilter );
1020 :
1021 27 : CPLFree( pOptions );
1022 27 : CSLDestroy( papszCreateOptions );
1023 27 : CSLDestroy( argv );
1024 27 : CSLDestroy( papszLayers );
1025 :
1026 27 : OGRCleanupAll();
1027 :
1028 27 : GDALDestroyDriverManager();
1029 :
1030 27 : return 0;
1031 : }
1032 :
|
