VeniVidiVici
/
gdal
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
gdal/apps/gdal2ogr.c

333 lines
12 KiB
C
Raw Permalink Blame History

/******************************************************************************
*
* Project: GDAL Utilities
* Purpose: Create an OGR datasource from the values of a GDAL dataset
* May be useful to test gdal_grid and generate its input OGR file
* Author: Even Rouault, <even dot rouault at spatialys.com>
*
******************************************************************************
* Copyright (c) 2008, Even Rouault <even dot rouault at spatialys.com>
*
* SPDX-License-Identifier: MIT
****************************************************************************/
#include "gdal.h"
#include "ogr_api.h"
#include "ogr_srs_api.h"
#include "cpl_string.h"
/************************************************************************/
/* Usage() */
/************************************************************************/
void Usage()
{
int iDriver;
int nDriverCount;
printf("Usage: gdal2ogr [--help] [--help-general] [-f format_name]\n"
" [-b band_number] [-l dest_layer_name]\n"
" [-t type]\n"
" gdal_datasource_src_name ogr_datasource_dst_name\n"
"\n"
" -f format_name: output file format name, possible values are:\n");
nDriverCount = OGRGetDriverCount();
for (iDriver = 0; iDriver < nDriverCount; iDriver++)
{
OGRSFDriverH hDriver = OGRGetDriver(iDriver);
if (OGR_Dr_TestCapability(hDriver, ODrCCreateDataSource))
printf(" -f \"%s\"\n", OGR_Dr_GetName(hDriver));
}
printf(
" -b band_number: band number of the GDAL datasource (1 by default)\n"
" -l dest_layer_name : name of the layer created in the OGR "
"datasource\n"
" (basename of the OGR datasource by default)\n"
" -t type: one of POINT, POINT25D (default), POLYGON\n"
"\n"
"Create an OGR datasource from the values of a GDAL dataset.\n\n");
exit(1);
}
/************************************************************************/
/* main() */
/************************************************************************/
int main(int argc, char *argv[])
{
const char *pszFormat = "ESRI Shapefile";
char *pszLayerName = NULL;
const char *pszSrcFilename = NULL, *pszDstFilename = NULL;
int iBand = 1;
GDALDatasetH hDS;
GDALRasterBandH hBand;
int nXSize, nYSize;
int i, j;
FILE *fOut = NULL;
double *padfBuffer;
double adfGeotransform[6];
OGRSFDriverH hOGRDriver;
OGRDataSourceH hOGRDS;
OGRLayerH hOGRLayer;
OGRwkbGeometryType eType = wkbPoint25D;
int xStep = 1, yStep = 1;
OGRRegisterAll();
GDALAllRegister();
argc = GDALGeneralCmdLineProcessor(argc, &argv, 0);
if (argc < 1)
exit(-argc);
/* -------------------------------------------------------------------- */
/* Parse arguments. */
/* -------------------------------------------------------------------- */
for (i = 1; i < argc; i++)
{
if (EQUAL(argv[i], "--help"))
Usage();
else if (EQUAL(argv[i], "-b") && i < argc - 1)
iBand = atoi(argv[++i]);
else if (EQUAL(argv[i], "-f") && i < argc - 1)
pszFormat = argv[++i];
else if (EQUAL(argv[i], "-l") && i < argc - 1)
pszLayerName = CPLStrdup(argv[++i]);
else if (EQUAL(argv[i], "-t") && i < argc - 1)
{
i++;
if (EQUAL(argv[i], "POLYGON"))
eType = wkbPolygon;
else if (EQUAL(argv[i], "POINT"))
eType = wkbPoint;
else if (EQUAL(argv[i], "POINT25D"))
eType = wkbPoint25D;
else
{
fprintf(stderr, "unhandled geometry type : %s\n", argv[i]);
}
}
else if (EQUAL(argv[i], "-step") && i < argc - 1)
xStep = yStep = atoi(argv[++i]);
else if (argv[i][0] == '-')
Usage();
else if (pszSrcFilename == NULL)
pszSrcFilename = argv[i];
else if (pszDstFilename == NULL)
pszDstFilename = argv[i];
else
Usage();
}
if (pszSrcFilename == NULL || pszDstFilename == NULL)
Usage();
/* -------------------------------------------------------------------- */
/* Open GDAL source dataset */
/* -------------------------------------------------------------------- */
hDS = GDALOpen(pszSrcFilename, GA_ReadOnly);
if (hDS == NULL)
{
fprintf(stderr, "Can't open %s\n", pszSrcFilename);
exit(1);
}
hBand = GDALGetRasterBand(hDS, iBand);
if (hBand == NULL)
{
fprintf(stderr, "Can't get band %d\n", iBand);
exit(1);
}
if (GDALGetGeoTransform(hDS, adfGeotransform) != CE_None)
{
fprintf(stderr, "Can't get geotransform\n");
exit(1);
}
nXSize = GDALGetRasterXSize(hDS);
nYSize = GDALGetRasterYSize(hDS);
/* -------------------------------------------------------------------- */
/* Create OGR destination dataset */
/* -------------------------------------------------------------------- */
/* Special case for CSV : we generate the appropriate VRT file in the same
* time */
if (EQUAL(pszFormat, "CSV") &&
EQUAL(CPLGetExtension(pszDstFilename), "CSV"))
{
FILE *fOutCSVT;
FILE *fOutVRT;
char *pszDstFilenameCSVT;
char *pszDstFilenameVRT;
fOut = fopen(pszDstFilename, "wt");
if (fOut == NULL)
{
fprintf(stderr, "Can't open %s for writing\n", pszDstFilename);
exit(1);
}
fprintf(fOut, "x,y,z\n");
pszDstFilenameCSVT = CPLMalloc(strlen(pszDstFilename) + 2);
strcpy(pszDstFilenameCSVT, pszDstFilename);
strcat(pszDstFilenameCSVT, "t");
fOutCSVT = fopen(pszDstFilenameCSVT, "wt");
if (fOutCSVT == NULL)
{
fprintf(stderr, "Can't open %s for writing\n", pszDstFilenameCSVT);
exit(1);
}
CPLFree(pszDstFilenameCSVT);
fprintf(fOutCSVT, "Real,Real,Real\n");
fclose(fOutCSVT);
fOutCSVT = NULL;
pszDstFilenameVRT = CPLStrdup(pszDstFilename);
strcpy(pszDstFilenameVRT + strlen(pszDstFilename) - 3, "vrt");
fOutVRT = fopen(pszDstFilenameVRT, "wt");
if (fOutVRT == NULL)
{
fprintf(stderr, "Can't open %s for writing\n", pszDstFilenameVRT);
exit(1);
}
CPLFree(pszDstFilenameVRT);
fprintf(fOutVRT, "<OGRVRTDataSource>\n");
fprintf(fOutVRT, " <OGRVRTLayer name=\"%s\">\n",
CPLGetBasename(pszDstFilename));
fprintf(fOutVRT, " <SrcDataSource>%s</SrcDataSource> \n",
pszDstFilename);
fprintf(fOutVRT, " <GeometryType>wkbPoint</GeometryType>\n");
fprintf(fOutVRT, " <GeometryField encoding=\"PointFromColumns\" "
"x=\"x\" y=\"y\" z=\"z\"/>\n");
fprintf(fOutVRT, " </OGRVRTLayer>\n");
fprintf(fOutVRT, "</OGRVRTDataSource>\n");
fclose(fOutVRT);
fOutVRT = NULL;
}
else
{
OGRSpatialReferenceH hSRS = NULL;
const char *pszWkt;
hOGRDriver = OGRGetDriverByName(pszFormat);
if (hOGRDriver == NULL)
{
fprintf(stderr, "Can't find OGR driver %s\n", pszFormat);
exit(1);
}
hOGRDS = OGR_Dr_CreateDataSource(hOGRDriver, pszDstFilename, NULL);
if (hOGRDS == NULL)
{
fprintf(stderr, "Can't create OGR datasource %s\n", pszDstFilename);
exit(1);
}
pszWkt = GDALGetProjectionRef(hDS);
if (pszWkt && pszWkt[0])
{
hSRS = OSRNewSpatialReference(pszWkt);
}
if (pszLayerName == NULL)
pszLayerName = CPLStrdup(CPLGetBasename(pszDstFilename));
hOGRLayer = OGR_DS_CreateLayer(hOGRDS, pszLayerName, hSRS, eType, NULL);
if (hSRS)
OSRDestroySpatialReference(hSRS);
if (hOGRLayer == NULL)
{
fprintf(stderr, "Can't create layer %s\n", pszLayerName);
exit(1);
}
if (eType != wkbPoint25D)
{
OGRFieldDefnH hFieldDefn = OGR_Fld_Create("z", OFTReal);
OGR_L_CreateField(hOGRLayer, hFieldDefn, 0);
OGR_Fld_Destroy(hFieldDefn);
}
}
padfBuffer = (double *)CPLMalloc(nXSize * sizeof(double));
#define GET_X(j, i) \
adfGeotransform[0] + (j)*adfGeotransform[1] + (i)*adfGeotransform[2]
#define GET_Y(j, i) \
adfGeotransform[3] + (j)*adfGeotransform[4] + (i)*adfGeotransform[5]
#define GET_XY(j, i) GET_X(j, i), GET_Y(j, i)
/* -------------------------------------------------------------------- */
/* "Translate" the source dataset */
/* -------------------------------------------------------------------- */
for (i = 0; i < nYSize; i += yStep)
{
GDALRasterIO(hBand, GF_Read, 0, i, nXSize, 1, padfBuffer, nXSize, 1,
GDT_Float64, 0, 0);
for (j = 0; j < nXSize; j += xStep)
{
if (fOut)
{
fprintf(fOut, "%f,%f,%f\n", GET_XY(j + .5, i + .5),
padfBuffer[j]);
}
else
{
OGRFeatureH hFeature =
OGR_F_Create(OGR_L_GetLayerDefn(hOGRLayer));
OGRGeometryH hGeometry = OGR_G_CreateGeometry(eType);
if (eType == wkbPoint25D)
{
OGR_G_SetPoint(hGeometry, 0, GET_XY(j + .5, i + .5),
padfBuffer[j]);
}
else if (eType == wkbPoint)
{
OGR_G_SetPoint_2D(hGeometry, 0, GET_XY(j + .5, i + .5));
OGR_F_SetFieldDouble(hFeature, 0, padfBuffer[j]);
}
else
{
OGRGeometryH hLinearRing =
OGR_G_CreateGeometry(wkbLinearRing);
OGR_G_SetPoint_2D(hLinearRing, 0, GET_XY(j + 0, i + 0));
OGR_G_SetPoint_2D(hLinearRing, 1, GET_XY(j + 1, i + 0));
OGR_G_SetPoint_2D(hLinearRing, 2, GET_XY(j + 1, i + 1));
OGR_G_SetPoint_2D(hLinearRing, 3, GET_XY(j + 0, i + 1));
OGR_G_SetPoint_2D(hLinearRing, 4, GET_XY(j + 0, i + 0));
OGR_G_AddGeometryDirectly(hGeometry, hLinearRing);
OGR_F_SetFieldDouble(hFeature, 0, padfBuffer[j]);
}
OGR_F_SetGeometryDirectly(hFeature, hGeometry);
OGR_L_CreateFeature(hOGRLayer, hFeature);
OGR_F_Destroy(hFeature);
}
}
}
/* -------------------------------------------------------------------- */
/* Cleanup */
/* -------------------------------------------------------------------- */
if (fOut)
fclose(fOut);
else
OGR_DS_Destroy(hOGRDS);
GDALClose(hDS);
CPLFree(padfBuffer);
CPLFree(pszLayerName);
GDALDumpOpenDatasets(stderr);
GDALDestroyDriverManager();
OGRCleanupAll();
CSLDestroy(argv);
return 0;
}
Reference in New Issue View Git Blame Copy Permalink