1 /******************************************************************************
2 *
3 * Project: GDAL Utilities
4 * Purpose: Create an OGR datasource from the values of a GDAL dataset
5 * May be useful to test gdal_grid and generate its input OGR file
6 * Author: Even Rouault, <even dot rouault at spatialys.com>
7 *
8 ******************************************************************************
9 * Copyright (c) 2008, Even Rouault <even dot rouault at spatialys.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
30 #include "gdal.h"
31 #include "ogr_api.h"
32 #include "ogr_srs_api.h"
33 #include "cpl_string.h"
34
35 CPL_CVSID("$Id: gdal2ogr.c 355b41831cd2685c85d1aabe5b95665a2c6e99b7 2019-06-19 17:07:04 +0200 Even Rouault $")
36
37 /************************************************************************/
38 /* Usage() */
39 /************************************************************************/
40
Usage()41 void Usage()
42 {
43 int iDriver;
44 int nDriverCount;
45
46 printf( "Usage: gdal2ogr [--help-general] [-f format_name]\n"
47 " [-b band_number] [-l dest_layer_name]\n"
48 " [-t type]\n"
49 " gdal_datasource_src_name ogr_datasource_dst_name\n"
50 "\n"
51 " -f format_name: output file format name, possible values are:\n");
52
53 nDriverCount = OGRGetDriverCount();
54 for( iDriver = 0; iDriver <nDriverCount; iDriver++ )
55 {
56 OGRSFDriverH hDriver = OGRGetDriver(iDriver);
57
58 if( OGR_Dr_TestCapability(hDriver, ODrCCreateDataSource ) )
59 printf( " -f \"%s\"\n", OGR_Dr_GetName(hDriver) );
60 }
61
62 printf( " -b band_number: band number of the GDAL datasource (1 by default)\n"
63 " -l dest_layer_name : name of the layer created in the OGR datasource\n"
64 " (basename of the OGR datasource by default)\n"
65 " -t type: one of POINT, POINT25D (default), POLYGON\n"
66 "\n"
67 "Create an OGR datasource from the values of a GDAL dataset.\n\n");
68
69 exit(1);
70 }
71
72 /************************************************************************/
73 /* main() */
74 /************************************************************************/
75
main(int argc,char * argv[])76 int main(int argc, char* argv[])
77 {
78 const char *pszFormat = "ESRI Shapefile";
79 char *pszLayerName = NULL;
80 const char *pszSrcFilename = NULL, *pszDstFilename = NULL;
81 int iBand = 1;
82 GDALDatasetH hDS;
83 GDALRasterBandH hBand;
84 int nXSize, nYSize;
85 int i, j;
86 FILE *fOut = NULL;
87 double *padfBuffer;
88 double adfGeotransform[6];
89 OGRSFDriverH hOGRDriver;
90 OGRDataSourceH hOGRDS;
91 OGRLayerH hOGRLayer;
92 OGRwkbGeometryType eType = wkbPoint25D;
93 int xStep = 1, yStep = 1;
94
95 OGRRegisterAll();
96 GDALAllRegister();
97
98 argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
99 if( argc < 1 )
100 exit( -argc );
101
102 /* -------------------------------------------------------------------- */
103 /* Parse arguments. */
104 /* -------------------------------------------------------------------- */
105 for( i = 1; i < argc; i++ )
106 {
107 if ( EQUAL(argv[i], "-b") && i < argc - 1)
108 iBand = atoi(argv[++i]);
109 else if ( EQUAL(argv[i], "-f") && i < argc - 1)
110 pszFormat = argv[++i];
111 else if ( EQUAL(argv[i], "-l") && i < argc - 1)
112 pszLayerName = CPLStrdup(argv[++i]);
113 else if ( EQUAL(argv[i], "-t") && i < argc - 1)
114 {
115 i++;
116 if (EQUAL(argv[i], "POLYGON"))
117 eType = wkbPolygon;
118 else if (EQUAL(argv[i], "POINT"))
119 eType = wkbPoint;
120 else if (EQUAL(argv[i], "POINT25D"))
121 eType = wkbPoint25D;
122 else
123 {
124 fprintf(stderr, "unhandled geometry type : %s\n", argv[i]);
125 }
126 }
127 else if ( EQUAL(argv[i], "-step") && i < argc - 1)
128 xStep = yStep = atoi(argv[++i]);
129 else if ( argv[i][0] == '-')
130 Usage();
131 else if( pszSrcFilename == NULL )
132 pszSrcFilename = argv[i];
133 else if( pszDstFilename == NULL )
134 pszDstFilename = argv[i];
135 else
136 Usage();
137 }
138
139 if( pszSrcFilename == NULL || pszDstFilename == NULL)
140 Usage();
141
142 /* -------------------------------------------------------------------- */
143 /* Open GDAL source dataset */
144 /* -------------------------------------------------------------------- */
145 hDS = GDALOpen(pszSrcFilename, GA_ReadOnly);
146 if (hDS == NULL)
147 {
148 fprintf(stderr, "Can't open %s\n", pszSrcFilename);
149 exit(1);
150 }
151
152 hBand = GDALGetRasterBand(hDS, iBand);
153 if (hBand == NULL)
154 {
155 fprintf(stderr, "Can't get band %d\n", iBand);
156 exit(1);
157 }
158
159 if (GDALGetGeoTransform(hDS, adfGeotransform) != CE_None)
160 {
161 fprintf(stderr, "Can't get geotransform\n");
162 exit(1);
163 }
164
165 nXSize = GDALGetRasterXSize(hDS);
166 nYSize = GDALGetRasterYSize(hDS);
167
168 /* -------------------------------------------------------------------- */
169 /* Create OGR destination dataset */
170 /* -------------------------------------------------------------------- */
171 /* Special case for CSV : we generate the appropriate VRT file in the same time */
172 if (EQUAL(pszFormat, "CSV") && EQUAL(CPLGetExtension(pszDstFilename), "CSV"))
173 {
174 FILE* fOutCSVT;
175 FILE* fOutVRT;
176 char* pszDstFilenameCSVT;
177 char* pszDstFilenameVRT;
178
179 fOut = fopen(pszDstFilename, "wt");
180 if (fOut == NULL)
181 {
182 fprintf(stderr, "Can't open %s for writing\n", pszDstFilename);
183 exit(1);
184 }
185 fprintf(fOut, "x,y,z\n");
186
187 pszDstFilenameCSVT = CPLMalloc(strlen(pszDstFilename) + 2);
188 strcpy(pszDstFilenameCSVT, pszDstFilename);
189 strcat(pszDstFilenameCSVT, "t");
190 fOutCSVT = fopen(pszDstFilenameCSVT, "wt");
191 if (fOutCSVT == NULL)
192 {
193 fprintf(stderr, "Can't open %s for writing\n", pszDstFilenameCSVT);
194 exit(1);
195 }
196 CPLFree(pszDstFilenameCSVT);
197 fprintf(fOutCSVT, "Real,Real,Real\n");
198 fclose(fOutCSVT);
199 fOutCSVT = NULL;
200
201 pszDstFilenameVRT = CPLStrdup(pszDstFilename);
202 strcpy(pszDstFilenameVRT + strlen(pszDstFilename) - 3, "vrt");
203 fOutVRT = fopen(pszDstFilenameVRT, "wt");
204 if (fOutVRT == NULL)
205 {
206 fprintf(stderr, "Can't open %s for writing\n", pszDstFilenameVRT);
207 exit(1);
208 }
209 CPLFree(pszDstFilenameVRT);
210 fprintf(fOutVRT, "<OGRVRTDataSource>\n");
211 fprintf(fOutVRT, " <OGRVRTLayer name=\"%s\">\n", CPLGetBasename(pszDstFilename));
212 fprintf(fOutVRT, " <SrcDataSource>%s</SrcDataSource> \n", pszDstFilename);
213 fprintf(fOutVRT, " <GeometryType>wkbPoint</GeometryType>\n");
214 fprintf(fOutVRT, " <GeometryField encoding=\"PointFromColumns\" x=\"x\" y=\"y\" z=\"z\"/>\n");
215 fprintf(fOutVRT, " </OGRVRTLayer>\n");
216 fprintf(fOutVRT, "</OGRVRTDataSource>\n");
217 fclose(fOutVRT);
218 fOutVRT = NULL;
219 }
220 else
221 {
222 OGRSpatialReferenceH hSRS = NULL;
223 const char* pszWkt;
224
225 hOGRDriver = OGRGetDriverByName(pszFormat);
226 if (hOGRDriver == NULL)
227 {
228 fprintf(stderr, "Can't find OGR driver %s\n", pszFormat);
229 exit(1);
230 }
231
232 hOGRDS = OGR_Dr_CreateDataSource(hOGRDriver, pszDstFilename, NULL);
233 if (hOGRDS == NULL)
234 {
235 fprintf(stderr, "Can't create OGR datasource %s\n", pszDstFilename);
236 exit(1);
237 }
238
239 pszWkt = GDALGetProjectionRef(hDS);
240 if (pszWkt && pszWkt[0])
241 {
242 hSRS = OSRNewSpatialReference(pszWkt);
243 }
244
245 if (pszLayerName == NULL)
246 pszLayerName = CPLStrdup(CPLGetBasename(pszDstFilename));
247
248 hOGRLayer = OGR_DS_CreateLayer( hOGRDS, pszLayerName,
249 hSRS, eType, NULL);
250
251 if (hSRS)
252 OSRDestroySpatialReference(hSRS);
253
254 if (hOGRLayer == NULL)
255 {
256 fprintf(stderr, "Can't create layer %s\n", pszLayerName);
257 exit(1);
258 }
259
260 if (eType != wkbPoint25D)
261 {
262 OGRFieldDefnH hFieldDefn = OGR_Fld_Create( "z", OFTReal );
263 OGR_L_CreateField(hOGRLayer, hFieldDefn, 0);
264 OGR_Fld_Destroy( hFieldDefn );
265 }
266 }
267
268
269 padfBuffer = (double*)CPLMalloc(nXSize * sizeof(double));
270
271 #define GET_X(j, i) adfGeotransform[0] + (j) * adfGeotransform[1] + (i) * adfGeotransform[2]
272 #define GET_Y(j, i) adfGeotransform[3] + (j) * adfGeotransform[4] + (i) * adfGeotransform[5]
273 #define GET_XY(j, i) GET_X(j, i), GET_Y(j, i)
274
275 /* -------------------------------------------------------------------- */
276 /* "Translate" the source dataset */
277 /* -------------------------------------------------------------------- */
278 for(i=0;i<nYSize;i+=yStep)
279 {
280 GDALRasterIO( hBand, GF_Read, 0, i, nXSize, 1,
281 padfBuffer, nXSize, 1, GDT_Float64, 0, 0);
282 for(j=0;j<nXSize;j+=xStep)
283 {
284 if (fOut)
285 {
286 fprintf(fOut, "%f,%f,%f\n",
287 GET_XY(j + .5, i + .5), padfBuffer[j]);
288 }
289 else
290 {
291 OGRFeatureH hFeature = OGR_F_Create(OGR_L_GetLayerDefn(hOGRLayer));
292 OGRGeometryH hGeometry = OGR_G_CreateGeometry(eType);
293 if (eType == wkbPoint25D)
294 {
295 OGR_G_SetPoint(hGeometry, 0, GET_XY(j + .5, i + .5),
296 padfBuffer[j]);
297 }
298 else if (eType == wkbPoint)
299 {
300 OGR_G_SetPoint_2D(hGeometry, 0, GET_XY(j + .5, i + .5));
301 OGR_F_SetFieldDouble(hFeature, 0, padfBuffer[j]);
302 }
303 else
304 {
305 OGRGeometryH hLinearRing = OGR_G_CreateGeometry(wkbLinearRing);
306 OGR_G_SetPoint_2D(hLinearRing, 0, GET_XY(j + 0, i + 0));
307 OGR_G_SetPoint_2D(hLinearRing, 1, GET_XY(j + 1, i + 0));
308 OGR_G_SetPoint_2D(hLinearRing, 2, GET_XY(j + 1, i + 1));
309 OGR_G_SetPoint_2D(hLinearRing, 3, GET_XY(j + 0, i + 1));
310 OGR_G_SetPoint_2D(hLinearRing, 4, GET_XY(j + 0, i + 0));
311 OGR_G_AddGeometryDirectly(hGeometry, hLinearRing);
312 OGR_F_SetFieldDouble(hFeature, 0, padfBuffer[j]);
313 }
314 OGR_F_SetGeometryDirectly(hFeature, hGeometry);
315 OGR_L_CreateFeature(hOGRLayer, hFeature);
316 OGR_F_Destroy(hFeature);
317 }
318 }
319 }
320
321 /* -------------------------------------------------------------------- */
322 /* Cleanup */
323 /* -------------------------------------------------------------------- */
324 if (fOut)
325 fclose(fOut);
326 else
327 OGR_DS_Destroy(hOGRDS);
328 GDALClose(hDS);
329
330 CPLFree(padfBuffer);
331 CPLFree(pszLayerName);
332
333 GDALDumpOpenDatasets( stderr );
334 GDALDestroyDriverManager();
335 OGRCleanupAll();
336 CSLDestroy( argv );
337
338 return 0;
339 }
340