PageRenderTime 41ms CodeModel.GetById 17ms RepoModel.GetById 0ms app.codeStats 0ms

/gdal-1.6.1-tcl_patched/apps/gdal_rasterize.cpp

http://tcl-map.googlecode.com/
C++ | 448 lines | 290 code | 64 blank | 94 comment | 96 complexity | 4676540eb91e66ea6a63440ad6c03de9 MD5 | raw file
Possible License(s): LGPL-2.0, GPL-3.0, BSD-3-Clause, LGPL-3.0
  1. /******************************************************************************
  2. * $Id: gdal_rasterize.cpp 15712 2008-11-12 13:59:56Z warmerdam $
  3. *
  4. * Project: GDAL Utilities
  5. * Purpose: Rasterize OGR shapes into a GDAL raster.
  6. * Author: Frank Warmerdam <warmerdam@pobox.com>
  7. *
  8. ******************************************************************************
  9. * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
  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. #include "gdal.h"
  30. #include "gdal_alg.h"
  31. #include "cpl_conv.h"
  32. #include "ogr_api.h"
  33. #include "ogr_srs_api.h"
  34. #include "cpl_string.h"
  35. #include <vector>
  36. CPL_CVSID("$Id: gdal_rasterize.cpp 15712 2008-11-12 13:59:56Z warmerdam $");
  37. /************************************************************************/
  38. /* Usage() */
  39. /************************************************************************/
  40. static void Usage()
  41. {
  42. printf(
  43. "Usage: gdal_rasterize [-b band] [-i]\n"
  44. " [-burn value] | [-a attribute_name] | [-3d]\n"
  45. // " [-of format_driver] [-co key=value]\n"
  46. // " [-te xmin ymin xmax ymax] [-tr xres yres] [-ts width height]\n"
  47. " [-l layername]* [-where expression] [-sql select_statement]\n"
  48. " <src_datasource> <dst_filename>\n" );
  49. exit( 1 );
  50. }
  51. /************************************************************************/
  52. /* InvertGeometries() */
  53. /************************************************************************/
  54. static void InvertGeometries( GDALDatasetH hDstDS,
  55. std::vector<OGRGeometryH> &ahGeometries )
  56. {
  57. OGRGeometryH hCollection =
  58. OGR_G_CreateGeometry( wkbGeometryCollection );
  59. /* -------------------------------------------------------------------- */
  60. /* Create a ring that is a bit outside the raster dataset. */
  61. /* -------------------------------------------------------------------- */
  62. OGRGeometryH hUniversePoly, hUniverseRing;
  63. double adfGeoTransform[6];
  64. int brx = GDALGetRasterXSize( hDstDS ) + 2;
  65. int bry = GDALGetRasterYSize( hDstDS ) + 2;
  66. GDALGetGeoTransform( hDstDS, adfGeoTransform );
  67. hUniverseRing = OGR_G_CreateGeometry( wkbLinearRing );
  68. OGR_G_AddPoint_2D(
  69. hUniverseRing,
  70. adfGeoTransform[0] + -2*adfGeoTransform[1] + -2*adfGeoTransform[2],
  71. adfGeoTransform[3] + -2*adfGeoTransform[4] + -2*adfGeoTransform[5] );
  72. OGR_G_AddPoint_2D(
  73. hUniverseRing,
  74. adfGeoTransform[0] + brx*adfGeoTransform[1] + -2*adfGeoTransform[2],
  75. adfGeoTransform[3] + brx*adfGeoTransform[4] + -2*adfGeoTransform[5] );
  76. OGR_G_AddPoint_2D(
  77. hUniverseRing,
  78. adfGeoTransform[0] + brx*adfGeoTransform[1] + bry*adfGeoTransform[2],
  79. adfGeoTransform[3] + brx*adfGeoTransform[4] + bry*adfGeoTransform[5] );
  80. OGR_G_AddPoint_2D(
  81. hUniverseRing,
  82. adfGeoTransform[0] + -2*adfGeoTransform[1] + bry*adfGeoTransform[2],
  83. adfGeoTransform[3] + -2*adfGeoTransform[4] + bry*adfGeoTransform[5] );
  84. OGR_G_AddPoint_2D(
  85. hUniverseRing,
  86. adfGeoTransform[0] + -2*adfGeoTransform[1] + -2*adfGeoTransform[2],
  87. adfGeoTransform[3] + -2*adfGeoTransform[4] + -2*adfGeoTransform[5] );
  88. hUniversePoly = OGR_G_CreateGeometry( wkbPolygon );
  89. OGR_G_AddGeometryDirectly( hUniversePoly, hUniverseRing );
  90. OGR_G_AddGeometryDirectly( hCollection, hUniversePoly );
  91. /* -------------------------------------------------------------------- */
  92. /* Add the rest of the geometries into our collection. */
  93. /* -------------------------------------------------------------------- */
  94. unsigned int iGeom;
  95. for( iGeom = 0; iGeom < ahGeometries.size(); iGeom++ )
  96. OGR_G_AddGeometryDirectly( hCollection, ahGeometries[iGeom] );
  97. ahGeometries.resize(1);
  98. ahGeometries[0] = hCollection;
  99. }
  100. /************************************************************************/
  101. /* ProcessLayer() */
  102. /* */
  103. /* Process all the features in a layer selection, collecting */
  104. /* geometries and burn values. */
  105. /************************************************************************/
  106. static void ProcessLayer(
  107. OGRLayerH hSrcLayer,
  108. GDALDatasetH hDstDS, std::vector<int> anBandList,
  109. std::vector<double> &adfBurnValues, int b3D, int bInverse,
  110. const char *pszBurnAttribute )
  111. {
  112. /* -------------------------------------------------------------------- */
  113. /* Checkout that SRS are the same. */
  114. /* -------------------------------------------------------------------- */
  115. OGRSpatialReferenceH hDstSRS = NULL;
  116. if( GDALGetProjectionRef( hDstDS ) != NULL )
  117. {
  118. char *pszProjection;
  119. pszProjection = (char *) GDALGetProjectionRef( hDstDS );
  120. hDstSRS = OSRNewSpatialReference(NULL);
  121. if( OSRImportFromWkt( hDstSRS, &pszProjection ) != CE_None )
  122. {
  123. OSRDestroySpatialReference(hDstSRS);
  124. hDstSRS = NULL;
  125. }
  126. }
  127. OGRSpatialReferenceH hSrcSRS = OGR_L_GetSpatialRef(hSrcLayer);
  128. if( hDstSRS != NULL && hSrcSRS != NULL )
  129. {
  130. if( OSRIsSame(hSrcSRS, hDstSRS) == FALSE )
  131. {
  132. fprintf(stderr,
  133. "Warning : the output raster dataset and the input vector layer do not have the same SRS. "
  134. "Results will be probably incorrect.\n");
  135. }
  136. }
  137. else if( hDstSRS != NULL && hSrcSRS == NULL )
  138. {
  139. fprintf(stderr,
  140. "Warning : the output raster dataset has a SRS, but the input vector layer not. "
  141. "Results will be probably incorrect.\n");
  142. }
  143. else if( hDstSRS == NULL && hSrcLayer != NULL )
  144. {
  145. fprintf(stderr,
  146. "Warning : the input vector layer has a SRS, but the output raster dataset not. "
  147. "Results will be probably incorrect.\n");
  148. }
  149. if( hDstSRS != NULL )
  150. {
  151. OSRDestroySpatialReference(hDstSRS);
  152. }
  153. /* -------------------------------------------------------------------- */
  154. /* Get field index, and check. */
  155. /* -------------------------------------------------------------------- */
  156. int iBurnField = -1;
  157. if( pszBurnAttribute )
  158. {
  159. iBurnField = OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( hSrcLayer ),
  160. pszBurnAttribute );
  161. if( iBurnField == -1 )
  162. {
  163. printf( "Failed to find field %s on layer %s, skipping.\n",
  164. pszBurnAttribute,
  165. OGR_FD_GetName( OGR_L_GetLayerDefn( hSrcLayer ) ) );
  166. return;
  167. }
  168. }
  169. /* -------------------------------------------------------------------- */
  170. /* Collect the geometries from this layer, and build list of */
  171. /* burn values. */
  172. /* -------------------------------------------------------------------- */
  173. OGRFeatureH hFeat;
  174. std::vector<OGRGeometryH> ahGeometries;
  175. std::vector<double> adfFullBurnValues;
  176. OGR_L_ResetReading( hSrcLayer );
  177. while( (hFeat = OGR_L_GetNextFeature( hSrcLayer )) != NULL )
  178. {
  179. OGRGeometryH hGeom;
  180. hGeom = OGR_G_Clone( OGR_F_GetGeometryRef( hFeat ) );
  181. ahGeometries.push_back( hGeom );
  182. for( unsigned int iBand = 0; iBand < anBandList.size(); iBand++ )
  183. {
  184. if( adfBurnValues.size() > 0 )
  185. adfFullBurnValues.push_back(
  186. adfBurnValues[MIN(iBand,adfBurnValues.size()-1)] );
  187. else if( b3D )
  188. {
  189. // TODO: get geometry "z" value
  190. adfFullBurnValues.push_back( 0.0 );
  191. }
  192. else if( pszBurnAttribute )
  193. {
  194. adfFullBurnValues.push_back( OGR_F_GetFieldAsDouble( hFeat, iBurnField ) );
  195. }
  196. }
  197. OGR_F_Destroy( hFeat );
  198. }
  199. /* -------------------------------------------------------------------- */
  200. /* If we are in inverse mode, we add one extra ring around the */
  201. /* whole dataset to invert the concept of insideness and then */
  202. /* merge everything into one geomtry collection. */
  203. /* -------------------------------------------------------------------- */
  204. if( bInverse )
  205. {
  206. InvertGeometries( hDstDS, ahGeometries );
  207. }
  208. /* -------------------------------------------------------------------- */
  209. /* Perform the burn. */
  210. /* -------------------------------------------------------------------- */
  211. GDALRasterizeGeometries( hDstDS, anBandList.size(), &(anBandList[0]),
  212. ahGeometries.size(), &(ahGeometries[0]),
  213. NULL, NULL, &(adfFullBurnValues[0]), NULL,
  214. GDALTermProgress, NULL );
  215. /* -------------------------------------------------------------------- */
  216. /* Cleanup geometries. */
  217. /* -------------------------------------------------------------------- */
  218. int iGeom;
  219. for( iGeom = ahGeometries.size()-1; iGeom >= 0; iGeom-- )
  220. OGR_G_DestroyGeometry( ahGeometries[iGeom] );
  221. }
  222. /************************************************************************/
  223. /* main() */
  224. /************************************************************************/
  225. int main( int argc, char ** argv )
  226. {
  227. int i, b3D = FALSE;
  228. int bInverse = FALSE;
  229. const char *pszSrcFilename = NULL;
  230. const char *pszDstFilename = NULL;
  231. char **papszLayers = NULL;
  232. const char *pszSQL = NULL;
  233. const char *pszBurnAttribute = NULL;
  234. const char *pszWHERE = NULL;
  235. std::vector<int> anBandList;
  236. std::vector<double> adfBurnValues;
  237. /* Check that we are running against at least GDAL 1.4 */
  238. /* Note to developers : if we use newer API, please change the requirement */
  239. if (atoi(GDALVersionInfo("VERSION_NUM")) < 1400)
  240. {
  241. fprintf(stderr, "At least, GDAL >= 1.4.0 is required for this version of %s, "
  242. "which was compiled against GDAL %s\n", argv[0], GDAL_RELEASE_NAME);
  243. exit(1);
  244. }
  245. GDALAllRegister();
  246. OGRRegisterAll();
  247. argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
  248. if( argc < 1 )
  249. exit( -argc );
  250. /* -------------------------------------------------------------------- */
  251. /* Parse arguments. */
  252. /* -------------------------------------------------------------------- */
  253. for( i = 1; i < argc; i++ )
  254. {
  255. if( EQUAL(argv[i], "--utility_version") )
  256. {
  257. printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
  258. argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
  259. return 0;
  260. }
  261. else if( EQUAL(argv[i],"-a") && i < argc-1 )
  262. {
  263. pszBurnAttribute = argv[++i];
  264. }
  265. else if( EQUAL(argv[i],"-b") && i < argc-1 )
  266. {
  267. anBandList.push_back( atoi(argv[++i]) );
  268. }
  269. else if( EQUAL(argv[i],"-3d") )
  270. {
  271. b3D = TRUE;
  272. }
  273. else if( EQUAL(argv[i],"-i") )
  274. {
  275. bInverse = TRUE;
  276. }
  277. else if( EQUAL(argv[i],"-burn") && i < argc-1 )
  278. {
  279. adfBurnValues.push_back( atof(argv[++i]) );
  280. }
  281. else if( EQUAL(argv[i],"-where") && i < argc-1 )
  282. {
  283. pszWHERE = argv[++i];
  284. }
  285. else if( EQUAL(argv[i],"-l") && i < argc-1 )
  286. {
  287. papszLayers = CSLAddString( papszLayers, argv[++i] );
  288. }
  289. else if( EQUAL(argv[i],"-sql") && i < argc-1 )
  290. {
  291. pszSQL = argv[++i];
  292. }
  293. else if( pszSrcFilename == NULL )
  294. {
  295. pszSrcFilename = argv[i];
  296. }
  297. else if( pszDstFilename == NULL )
  298. {
  299. pszDstFilename = argv[i];
  300. }
  301. else
  302. Usage();
  303. }
  304. if( pszSrcFilename == NULL || pszDstFilename == NULL )
  305. {
  306. fprintf( stderr, "Missing source or destination.\n\n" );
  307. Usage();
  308. }
  309. if( pszSQL == NULL && papszLayers == NULL )
  310. {
  311. fprintf( stderr, "At least one of -l or -sql required.\n\n" );
  312. Usage();
  313. }
  314. if( adfBurnValues.size() == 0 && pszBurnAttribute == NULL && !b3D )
  315. {
  316. fprintf( stderr, "At least one of -3d, -burn or -a required.\n\n" );
  317. Usage();
  318. }
  319. if( anBandList.size() == 0 )
  320. anBandList.push_back( 1 );
  321. /* -------------------------------------------------------------------- */
  322. /* Open source vector dataset. */
  323. /* -------------------------------------------------------------------- */
  324. OGRDataSourceH hSrcDS;
  325. hSrcDS = OGROpen( pszSrcFilename, FALSE, NULL );
  326. if( hSrcDS == NULL )
  327. {
  328. fprintf( stderr, "Failed to open feature source: %s\n",
  329. pszSrcFilename);
  330. exit( 1 );
  331. }
  332. /* -------------------------------------------------------------------- */
  333. /* Open target raster file. Eventually we will add optional */
  334. /* creation. */
  335. /* -------------------------------------------------------------------- */
  336. GDALDatasetH hDstDS;
  337. hDstDS = GDALOpen( pszDstFilename, GA_Update );
  338. if( hDstDS == NULL )
  339. exit( 2 );
  340. /* -------------------------------------------------------------------- */
  341. /* Process SQL request. */
  342. /* -------------------------------------------------------------------- */
  343. if( pszSQL != NULL )
  344. {
  345. OGRLayerH hLayer;
  346. hLayer = OGR_DS_ExecuteSQL( hSrcDS, pszSQL, NULL, NULL );
  347. if( hLayer != NULL )
  348. {
  349. ProcessLayer( hLayer, hDstDS, anBandList,
  350. adfBurnValues, b3D, bInverse, pszBurnAttribute );
  351. }
  352. }
  353. /* -------------------------------------------------------------------- */
  354. /* Process each layer. */
  355. /* -------------------------------------------------------------------- */
  356. int nLayerCount = CSLCount(papszLayers);
  357. for( i = 0; i < nLayerCount; i++ )
  358. {
  359. OGRLayerH hLayer = OGR_DS_GetLayerByName( hSrcDS, papszLayers[i] );
  360. if( hLayer == NULL )
  361. {
  362. fprintf( stderr, "Unable to find layer %s, skipping.\n",
  363. papszLayers[i] );
  364. continue;
  365. }
  366. if( pszWHERE )
  367. {
  368. if( OGR_L_SetAttributeFilter( hLayer, pszWHERE ) != OGRERR_NONE )
  369. break;
  370. }
  371. ProcessLayer( hLayer, hDstDS, anBandList,
  372. adfBurnValues, b3D, bInverse, pszBurnAttribute );
  373. }
  374. /* -------------------------------------------------------------------- */
  375. /* Cleanup */
  376. /* -------------------------------------------------------------------- */
  377. OGR_DS_Destroy( hSrcDS );
  378. GDALClose( hDstDS );
  379. CSLDestroy( argv );
  380. CSLDestroy( papszLayers );
  381. GDALDestroyDriverManager();
  382. OGRCleanupAll();
  383. return 0;
  384. }