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