|
1 |
/******************************************************************************
|
|
2 |
* $Id: gdalwarp.cpp,v 1.27 2006/09/27 13:06:13 dron Exp $
|
|
3 |
*
|
|
4 |
* Project: High Performance Image Reprojector
|
|
5 |
* Purpose: Test program for high performance warper API.
|
|
6 |
* Author: Frank Warmerdam <warmerdam@pobox.com>
|
|
7 |
*
|
|
8 |
******************************************************************************
|
|
9 |
* Copyright (c) 2002, i3 - information integration and imaging
|
|
10 |
* Fort Collin, CO
|
|
11 |
*
|
|
12 |
* Permission is hereby granted, free of charge, to any person obtaining a
|
|
13 |
* copy of this software and associated documentation files (the "Software"),
|
|
14 |
* to deal in the Software without restriction, including without limitation
|
|
15 |
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
|
|
16 |
* and/or sell copies of the Software, and to permit persons to whom the
|
|
17 |
* Software is furnished to do so, subject to the following conditions:
|
|
18 |
*
|
|
19 |
* The above copyright notice and this permission notice shall be included
|
|
20 |
* in all copies or substantial portions of the Software.
|
|
21 |
*
|
|
22 |
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
|
|
23 |
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
24 |
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
|
|
25 |
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
26 |
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|
27 |
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
|
|
28 |
* DEALINGS IN THE SOFTWARE.
|
|
29 |
******************************************************************************
|
|
30 |
*
|
|
31 |
* $Log: gdalwarp.cpp,v $
|
|
32 |
* Revision 1.27.1 2007/04/02 11:35:08 maquerol
|
|
33 |
* main() function converted into an accesible warpFunction.
|
|
34 |
*
|
|
35 |
* Revision 1.27 2006/09/27 13:06:13 dron
|
|
36 |
* Memory leak fixed.
|
|
37 |
*
|
|
38 |
* Revision 1.26 2006/07/06 20:30:11 fwarmerdam
|
|
39 |
* use GDALSuggestedWarpOutput2() to avoid approximation implicit in using gt
|
|
40 |
*
|
|
41 |
* Revision 1.25 2006/06/29 21:10:29 fwarmerdam
|
|
42 |
* Avoid a few memory leaks.
|
|
43 |
*
|
|
44 |
* Revision 1.24 2006/06/02 17:31:49 fwarmerdam
|
|
45 |
* Modified -ts to allow width or height to be zero meaning, it should
|
|
46 |
* be computed to retain square pixels.
|
|
47 |
*
|
|
48 |
* Revision 1.23 2006/05/29 17:32:15 fwarmerdam
|
|
49 |
* added preliminary support for controlling longitude wrapping on input dataset
|
|
50 |
*
|
|
51 |
* Revision 1.22 2006/04/25 14:28:33 fwarmerdam
|
|
52 |
* Check for no usable sources, avoid warning.
|
|
53 |
*
|
|
54 |
* Revision 1.21 2006/03/21 21:34:43 fwarmerdam
|
|
55 |
* cleanup headers
|
|
56 |
*
|
|
57 |
* Revision 1.20 2006/01/05 19:51:10 fwarmerdam
|
|
58 |
* fixed checking of targetsrs
|
|
59 |
*
|
|
60 |
* Revision 1.19 2005/12/19 20:20:00 fwarmerdam
|
|
61 |
* preliminary support for multiple input files
|
|
62 |
*
|
|
63 |
* Revision 1.18 2005/09/13 01:57:31 fwarmerdam
|
|
64 |
* Fixed usage message.
|
|
65 |
*
|
|
66 |
* Revision 1.17 2005/09/13 01:24:45 fwarmerdam
|
|
67 |
* Set UNIFIED_SRC_NODATA by default for -srcnodata switch.
|
|
68 |
*
|
|
69 |
* Revision 1.16 2005/09/12 18:06:25 fwarmerdam
|
|
70 |
* Added quotes around -srcnodata and -dstnodata in usage message.
|
|
71 |
*
|
|
72 |
* Revision 1.15 2004/12/26 16:14:53 fwarmerdam
|
|
73 |
* added -tps flag
|
|
74 |
*
|
|
75 |
* Revision 1.14 2004/11/29 15:01:55 fwarmerdam
|
|
76 |
* Fixed typo in printf as per Bug 691.
|
|
77 |
*
|
|
78 |
* Revision 1.13 2004/11/14 04:57:04 fwarmerdam
|
|
79 |
* added -srcalpha switch, and automatic alpha detection
|
|
80 |
*
|
|
81 |
* Revision 1.12 2004/11/05 06:15:08 fwarmerdam
|
|
82 |
* Don't double free the warpoptions array.
|
|
83 |
*
|
|
84 |
* Revision 1.11 2004/11/05 05:53:43 fwarmerdam
|
|
85 |
* Avoid various memory leaks.
|
|
86 |
*
|
|
87 |
* Revision 1.10 2004/10/07 15:53:42 fwarmerdam
|
|
88 |
* added preliminary alpha band support
|
|
89 |
*/
|
|
90 |
#include <jni.h>
|
|
91 |
#include "gdalwarper.h"
|
|
92 |
#include "cpl_string.h"
|
|
93 |
#include "ogr_spatialref.h"
|
|
94 |
#include "gdal.h"
|
|
95 |
//#include "../include/gdalwarp_interfaz.h"
|
|
96 |
|
|
97 |
CPL_CVSID("$Id: gdalwarp.cpp,v 1.27 2006/09/27 13:06:13 dron Exp $");
|
|
98 |
|
|
99 |
static CPLString InsertCenterLong( GDALDatasetH hDS, CPLString osWKT );
|
|
100 |
|
|
101 |
static GDALDatasetH
|
|
102 |
GDALWarpCreateOutput( char **papszSrcFiles, const char *pszFilename,
|
|
103 |
const char *pszFormat, const char *pszSourceSRS,
|
|
104 |
const char *pszTargetSRS, int nOrder,
|
|
105 |
char **papszCreateOptions, GDALDataType eDT );
|
|
106 |
|
|
107 |
int CPL_STDCALL GDALFuncTermProgress( double dfComplete, const char *pszMessage,
|
|
108 |
void * pProgressArg );
|
|
109 |
|
|
110 |
void CPL_STDCALL statusCallBack(int nPercent);
|
|
111 |
|
|
112 |
static double dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0;
|
|
113 |
static double dfXRes=0.0, dfYRes=0.0;
|
|
114 |
static int nForcePixels=0, nForceLines=0, bQuiet = FALSE;
|
|
115 |
static int bEnableDstAlpha = FALSE, bEnableSrcAlpha = FALSE;
|
|
116 |
|
|
117 |
static int bVRT = FALSE;
|
|
118 |
|
|
119 |
|
|
120 |
/************************************************************************/
|
|
121 |
/* SanitizeSRS */
|
|
122 |
/************************************************************************/
|
|
123 |
|
|
124 |
char *SanitizeSRS( const char *pszUserInput )
|
|
125 |
|
|
126 |
{
|
|
127 |
OGRSpatialReferenceH hSRS;
|
|
128 |
char *pszResult = NULL;
|
|
129 |
|
|
130 |
CPLErrorReset();
|
|
131 |
|
|
132 |
hSRS = OSRNewSpatialReference( NULL );
|
|
133 |
if( OSRSetFromUserInput( hSRS, pszUserInput ) == OGRERR_NONE )
|
|
134 |
OSRExportToWkt( hSRS, &pszResult );
|
|
135 |
else
|
|
136 |
{
|
|
137 |
CPLError( CE_Failure, CPLE_AppDefined,
|
|
138 |
"Translating source or target SRS failed:\n%s",
|
|
139 |
pszUserInput );
|
|
140 |
exit( 1 );
|
|
141 |
}
|
|
142 |
|
|
143 |
OSRDestroySpatialReference( hSRS );
|
|
144 |
|
|
145 |
return pszResult;
|
|
146 |
}
|
|
147 |
|
|
148 |
|
|
149 |
/************************************************************************/
|
|
150 |
/* WarpFunction() */
|
|
151 |
/************************************************************************/
|
|
152 |
|
|
153 |
int launchException(JNIEnv *env, char *msg) {
|
|
154 |
env->ExceptionDescribe();
|
|
155 |
env->ExceptionClear();
|
|
156 |
|
|
157 |
jclass newExcCls = env->FindClass("es/gva/cit/jgdal/GdalException");
|
|
158 |
if (newExcCls == 0)
|
|
159 |
return CE_Failure;
|
|
160 |
|
|
161 |
env->ThrowNew(newExcCls, msg);
|
|
162 |
return CE_Failure;
|
|
163 |
}
|
|
164 |
|
|
165 |
int warpFunction(JNIEnv *env, char * s_srs, char * t_srs, char *source, char * dest, char *format) {
|
|
166 |
GDALDatasetH hDstDS;
|
|
167 |
const char *pszFormat = "GTiff";
|
|
168 |
char *pszTargetSRS = NULL;
|
|
169 |
char *pszSourceSRS = NULL;
|
|
170 |
char **papszSrcFiles = NULL;
|
|
171 |
char *pszDstFilename = NULL;
|
|
172 |
int bCreateOutput = FALSE, i, nOrder = 0;
|
|
173 |
void *hTransformArg, *hGenImgProjArg = NULL, *hApproxArg = NULL;
|
|
174 |
char **papszWarpOptions = NULL;
|
|
175 |
double dfErrorThreshold = 0.125;
|
|
176 |
double dfWarpMemoryLimit = 0.0;
|
|
177 |
GDALTransformerFunc pfnTransformer = NULL;
|
|
178 |
char **papszCreateOptions = NULL;
|
|
179 |
GDALDataType eOutputType = GDT_Unknown, eWorkingType = GDT_Unknown;
|
|
180 |
GDALResampleAlg eResampleAlg = GRA_NearestNeighbour;
|
|
181 |
const char *pszSrcNodata = NULL;
|
|
182 |
const char *pszDstNodata = NULL;
|
|
183 |
int bMulti = FALSE;
|
|
184 |
char error[1024];
|
|
185 |
int result = 4;
|
|
186 |
|
|
187 |
GDALAllRegister();
|
|
188 |
|
|
189 |
/* -------------------------------------------------------------------- */
|
|
190 |
/* Parse arguments. */
|
|
191 |
/* -------------------------------------------------------------------- */
|
|
192 |
pszTargetSRS = SanitizeSRS(t_srs);
|
|
193 |
if (s_srs != NULL)
|
|
194 |
pszSourceSRS = SanitizeSRS(s_srs);
|
|
195 |
|
|
196 |
papszSrcFiles = CSLAddString(papszSrcFiles, source);
|
|
197 |
papszSrcFiles = CSLAddString(papszSrcFiles, dest);
|
|
198 |
if (format != NULL)
|
|
199 |
pszFormat = format;
|
|
200 |
bCreateOutput = FALSE;
|
|
201 |
|
|
202 |
/* -------------------------------------------------------------------- */
|
|
203 |
/* The last filename in the file list is really our destination */
|
|
204 |
/* file. */
|
|
205 |
/* -------------------------------------------------------------------- */
|
|
206 |
if (CSLCount(papszSrcFiles) > 1) {
|
|
207 |
pszDstFilename = papszSrcFiles[CSLCount(papszSrcFiles) - 1];
|
|
208 |
papszSrcFiles[CSLCount(papszSrcFiles) - 1] = NULL;
|
|
209 |
}
|
|
210 |
|
|
211 |
/* -------------------------------------------------------------------- */
|
|
212 |
/* Does the output dataset already exist? */
|
|
213 |
/* -------------------------------------------------------------------- */
|
|
214 |
CPLPushErrorHandler(CPLQuietErrorHandler);
|
|
215 |
hDstDS = GDALOpen(pszDstFilename, GA_Update);
|
|
216 |
CPLPopErrorHandler();
|
|
217 |
|
|
218 |
if (hDstDS != NULL && bCreateOutput) {
|
|
219 |
sprintf(error, "warpfunc.cpp (218): Output dataset %s exists, "
|
|
220 |
"but some commandline options were provided indicating a new dataset "
|
|
221 |
"should be created. Please delete existing dataset and run again.", pszDstFilename);
|
|
222 |
return launchException(env, error);
|
|
223 |
}
|
|
224 |
|
|
225 |
/* -------------------------------------------------------------------- */
|
|
226 |
/* If not, we need to create it. */
|
|
227 |
/* -------------------------------------------------------------------- */
|
|
228 |
int bInitDestSetForFirst = FALSE;
|
|
229 |
|
|
230 |
if (hDstDS == NULL) {
|
|
231 |
hDstDS = GDALWarpCreateOutput(papszSrcFiles, pszDstFilename, pszFormat,
|
|
232 |
pszSourceSRS, pszTargetSRS, nOrder, papszCreateOptions, eOutputType);
|
|
233 |
bCreateOutput = TRUE;
|
|
234 |
|
|
235 |
if (CSLFetchNameValue(papszWarpOptions, "INIT_DEST") == NULL && pszDstNodata == NULL) {
|
|
236 |
papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "0");
|
|
237 |
bInitDestSetForFirst = TRUE;
|
|
238 |
} else {
|
|
239 |
if (CSLFetchNameValue(papszWarpOptions, "INIT_DEST") == NULL) {
|
|
240 |
papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "NO_DATA");
|
|
241 |
bInitDestSetForFirst = TRUE;
|
|
242 |
}
|
|
243 |
}
|
|
244 |
|
|
245 |
CSLDestroy(papszCreateOptions);
|
|
246 |
papszCreateOptions = NULL;
|
|
247 |
}
|
|
248 |
|
|
249 |
if (hDstDS == NULL)
|
|
250 |
return launchException(env, "warpfunc.cpp (249):");
|
|
251 |
|
|
252 |
if (pszTargetSRS == NULL)
|
|
253 |
pszTargetSRS = CPLStrdup(GDALGetProjectionRef(hDstDS));
|
|
254 |
|
|
255 |
/* -------------------------------------------------------------------- */
|
|
256 |
/* Loop over all source files, processing each in turn. */
|
|
257 |
/* -------------------------------------------------------------------- */
|
|
258 |
for (int iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++) {
|
|
259 |
CPLString osThisSourceSRS;
|
|
260 |
GDALDatasetH hSrcDS;
|
|
261 |
|
|
262 |
if (pszSourceSRS != NULL)
|
|
263 |
osThisSourceSRS = pszSourceSRS;
|
|
264 |
|
|
265 |
/* -------------------------------------------------------------------- */
|
|
266 |
/* Open this file. */
|
|
267 |
/* -------------------------------------------------------------------- */
|
|
268 |
hSrcDS = GDALOpen(papszSrcFiles[iSrc], GA_ReadOnly);
|
|
269 |
|
|
270 |
if (hSrcDS == NULL)
|
|
271 |
return launchException(env, "warpfunc.cpp (270):");
|
|
272 |
|
|
273 |
// Processing input file papszSrcFiles[iSrc]
|
|
274 |
if (strlen(osThisSourceSRS) == 0) {
|
|
275 |
if (GDALGetProjectionRef(hSrcDS) != NULL && strlen(GDALGetProjectionRef(hSrcDS)) > 0) {
|
|
276 |
osThisSourceSRS = GDALGetProjectionRef(hSrcDS);
|
|
277 |
} else {
|
|
278 |
if (GDALGetGCPProjection(hSrcDS) != NULL && strlen(GDALGetGCPProjection(hSrcDS)) > 0 && GDALGetGCPCount(hSrcDS) > 1) {
|
|
279 |
osThisSourceSRS = GDALGetGCPProjection(hSrcDS);
|
|
280 |
} else {
|
|
281 |
osThisSourceSRS = "";
|
|
282 |
}
|
|
283 |
}
|
|
284 |
|
|
285 |
if (pszTargetSRS != NULL && strlen(pszTargetSRS) > 0 && strlen(osThisSourceSRS) == 0) {
|
|
286 |
return launchException(env, "warpfunc.cpp (285): A target coordinate system was specified, "
|
|
287 |
"but there is no source coordinate system. Consider using -s_srs option to provide a "
|
|
288 |
"source coordinate system. Operation terminated.");
|
|
289 |
}
|
|
290 |
}
|
|
291 |
|
|
292 |
/* -------------------------------------------------------------------- */
|
|
293 |
/* Do we have a source alpha band? */
|
|
294 |
/* -------------------------------------------------------------------- */
|
|
295 |
if (GDALGetRasterColorInterpretation(GDALGetRasterBand(hSrcDS, GDALGetRasterCount(hSrcDS))) == GCI_AlphaBand && !bEnableSrcAlpha) {
|
|
296 |
// Using band GDALGetRasterCount(hSrcDS) of source image as alpha
|
|
297 |
bEnableSrcAlpha = TRUE;
|
|
298 |
}
|
|
299 |
|
|
300 |
/* -------------------------------------------------------------------- */
|
|
301 |
/* If the source coordinate system is geographic, try and */
|
|
302 |
/* insert a CENTER_LONG extension parameter on the GEOGCS to */
|
|
303 |
/* handle wrapping better. */
|
|
304 |
/* -------------------------------------------------------------------- */
|
|
305 |
if (EQUALN(osThisSourceSRS.c_str(), "GEOGCS[", 7))
|
|
306 |
osThisSourceSRS = InsertCenterLong(hSrcDS, osThisSourceSRS);
|
|
307 |
|
|
308 |
/* -------------------------------------------------------------------- */
|
|
309 |
/* Create a transformation object from the source to */
|
|
310 |
/* destination coordinate system. */
|
|
311 |
/* -------------------------------------------------------------------- */
|
|
312 |
hTransformArg = hGenImgProjArg = GDALCreateGenImgProjTransformer(hSrcDS, osThisSourceSRS, hDstDS,
|
|
313 |
pszTargetSRS, TRUE, 1000.0, nOrder);
|
|
314 |
|
|
315 |
if (hTransformArg == NULL)
|
|
316 |
return launchException(env, "warpfunc.cpp (315):");
|
|
317 |
|
|
318 |
pfnTransformer = GDALGenImgProjTransform;
|
|
319 |
|
|
320 |
/* -------------------------------------------------------------------- */
|
|
321 |
/* Warp the transformer with a linear approximator unless the */
|
|
322 |
/* acceptable error is zero. */
|
|
323 |
/* -------------------------------------------------------------------- */
|
|
324 |
if (dfErrorThreshold != 0.0) {
|
|
325 |
hTransformArg = hApproxArg = GDALCreateApproxTransformer(GDALGenImgProjTransform, hGenImgProjArg, dfErrorThreshold);
|
|
326 |
pfnTransformer = GDALApproxTransform;
|
|
327 |
}
|
|
328 |
|
|
329 |
/* -------------------------------------------------------------------- */
|
|
330 |
/* Clear temporary INIT_DEST settings after the first image. */
|
|
331 |
/* -------------------------------------------------------------------- */
|
|
332 |
if (bInitDestSetForFirst && iSrc == 1)
|
|
333 |
papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", NULL);
|
|
334 |
|
|
335 |
/* -------------------------------------------------------------------- */
|
|
336 |
/* Setup warp options. */
|
|
337 |
/* -------------------------------------------------------------------- */
|
|
338 |
GDALWarpOptions *psWO = GDALCreateWarpOptions();
|
|
339 |
|
|
340 |
psWO->papszWarpOptions = CSLDuplicate(papszWarpOptions);
|
|
341 |
psWO->eWorkingDataType = eWorkingType;
|
|
342 |
psWO->eResampleAlg = eResampleAlg;
|
|
343 |
|
|
344 |
psWO->hSrcDS = hSrcDS;
|
|
345 |
psWO->hDstDS = hDstDS;
|
|
346 |
|
|
347 |
psWO->pfnTransformer = pfnTransformer;
|
|
348 |
psWO->pTransformerArg = hTransformArg;
|
|
349 |
|
|
350 |
if (!bQuiet)
|
|
351 |
psWO->pfnProgress = GDALFuncTermProgress;
|
|
352 |
|
|
353 |
if (dfWarpMemoryLimit != 0.0)
|
|
354 |
psWO->dfWarpMemoryLimit = dfWarpMemoryLimit;
|
|
355 |
|
|
356 |
/* -------------------------------------------------------------------- */
|
|
357 |
/* Setup band mapping. */
|
|
358 |
/* -------------------------------------------------------------------- */
|
|
359 |
if (bEnableSrcAlpha)
|
|
360 |
psWO->nBandCount = GDALGetRasterCount(hSrcDS) - 1;
|
|
361 |
else
|
|
362 |
psWO->nBandCount = GDALGetRasterCount(hSrcDS);
|
|
363 |
|
|
364 |
psWO->panSrcBands = (int *) CPLMalloc(psWO->nBandCount * sizeof(int));
|
|
365 |
psWO->panDstBands = (int *) CPLMalloc(psWO->nBandCount * sizeof(int));
|
|
366 |
|
|
367 |
for (i = 0; i < psWO->nBandCount; i++) {
|
|
368 |
psWO->panSrcBands[i] = i + 1;
|
|
369 |
psWO->panDstBands[i] = i + 1;
|
|
370 |
}
|
|
371 |
|
|
372 |
/* -------------------------------------------------------------------- */
|
|
373 |
/* Setup alpha bands used if any. */
|
|
374 |
/* -------------------------------------------------------------------- */
|
|
375 |
if (bEnableSrcAlpha)
|
|
376 |
psWO->nSrcAlphaBand = GDALGetRasterCount(hSrcDS);
|
|
377 |
|
|
378 |
if (!bEnableDstAlpha && GDALGetRasterCount(hDstDS) == psWO->nBandCount + 1 &&
|
|
379 |
GDALGetRasterColorInterpretation(GDALGetRasterBand(hDstDS, GDALGetRasterCount(hDstDS))) == GCI_AlphaBand) {
|
|
380 |
// Using band GDALGetRasterCount(hDstDS) of destination image as alpha
|
|
381 |
bEnableDstAlpha = TRUE;
|
|
382 |
}
|
|
383 |
|
|
384 |
if (bEnableDstAlpha)
|
|
385 |
psWO->nDstAlphaBand = GDALGetRasterCount(hDstDS);
|
|
386 |
|
|
387 |
/* -------------------------------------------------------------------- */
|
|
388 |
/* Setup NODATA options. */
|
|
389 |
/* -------------------------------------------------------------------- */
|
|
390 |
if (pszSrcNodata != NULL) {
|
|
391 |
char **papszTokens = CSLTokenizeString(pszSrcNodata);
|
|
392 |
int nTokenCount = CSLCount(papszTokens);
|
|
393 |
|
|
394 |
psWO->padfSrcNoDataReal = (double *) CPLMalloc(psWO->nBandCount * sizeof(double));
|
|
395 |
psWO->padfSrcNoDataImag = (double *) CPLMalloc(psWO->nBandCount * sizeof(double));
|
|
396 |
|
|
397 |
for (i = 0; i < psWO->nBandCount; i++) {
|
|
398 |
if (i < nTokenCount) {
|
|
399 |
CPLStringToComplex(papszTokens[i], psWO->padfSrcNoDataReal + i, psWO->padfSrcNoDataImag + i);
|
|
400 |
} else {
|
|
401 |
psWO->padfSrcNoDataReal[i] = psWO->padfSrcNoDataReal[i - 1];
|
|
402 |
psWO->padfSrcNoDataImag[i] = psWO->padfSrcNoDataImag[i - 1];
|
|
403 |
}
|
|
404 |
}
|
|
405 |
|
|
406 |
CSLDestroy(papszTokens);
|
|
407 |
|
|
408 |
psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions, "UNIFIED_SRC_NODATA", "YES");
|
|
409 |
}
|
|
410 |
|
|
411 |
/* -------------------------------------------------------------------- */
|
|
412 |
/* If the output dataset was created, and we have a destination */
|
|
413 |
/* nodata value, go through marking the bands with the information.*/
|
|
414 |
/* -------------------------------------------------------------------- */
|
|
415 |
if (pszDstNodata != NULL && bCreateOutput) {
|
|
416 |
char **papszTokens = CSLTokenizeString(pszDstNodata);
|
|
417 |
int nTokenCount = CSLCount(papszTokens);
|
|
418 |
|
|
419 |
psWO->padfDstNoDataReal = (double *) CPLMalloc(psWO->nBandCount * sizeof(double));
|
|
420 |
psWO->padfDstNoDataImag = (double *) CPLMalloc(psWO->nBandCount * sizeof(double));
|
|
421 |
|
|
422 |
for (i = 0; i < psWO->nBandCount; i++) {
|
|
423 |
if (i < nTokenCount) {
|
|
424 |
CPLStringToComplex(papszTokens[i], psWO->padfDstNoDataReal + i, psWO->padfDstNoDataImag + i);
|
|
425 |
} else {
|
|
426 |
psWO->padfDstNoDataReal[i] = psWO->padfDstNoDataReal[i - 1];
|
|
427 |
psWO->padfDstNoDataImag[i] = psWO->padfDstNoDataImag[i - 1];
|
|
428 |
}
|
|
429 |
|
|
430 |
if (bCreateOutput) {
|
|
431 |
GDALSetRasterNoDataValue(GDALGetRasterBand(hDstDS, psWO->panDstBands[i]), psWO->padfDstNoDataReal[i]);
|
|
432 |
}
|
|
433 |
}
|
|
434 |
|
|
435 |
CSLDestroy(papszTokens);
|
|
436 |
}
|
|
437 |
|
|
438 |
/* -------------------------------------------------------------------- */
|
|
439 |
/* If we are producing VRT output, then just initialize it with */
|
|
440 |
/* the warp options and write out now rather than proceeding */
|
|
441 |
/* with the operations. */
|
|
442 |
/* -------------------------------------------------------------------- */
|
|
443 |
if (bVRT) {
|
|
444 |
if (GDALInitializeWarpedVRT(hDstDS, psWO) != CE_None)
|
|
445 |
return launchException(env, "warpfunc.cpp (444):");
|
|
446 |
|
|
447 |
GDALClose(hDstDS);
|
|
448 |
GDALClose(hSrcDS);
|
|
449 |
|
|
450 |
//GDALDumpOpenDatasets( stderr );
|
|
451 |
|
|
452 |
GDALDestroyDriverManager();
|
|
453 |
|
|
454 |
return CE_None;
|
|
455 |
}
|
|
456 |
|
|
457 |
/* -------------------------------------------------------------------- */
|
|
458 |
/* Initialize and execute the warp. */
|
|
459 |
/* -------------------------------------------------------------------- */
|
|
460 |
GDALWarpOperation oWO;
|
|
461 |
result = CE_None;
|
|
462 |
if (oWO.Initialize(psWO) == CE_None) {
|
|
463 |
if (bMulti) {
|
|
464 |
result = oWO.ChunkAndWarpMulti(0, 0, GDALGetRasterXSize(hDstDS), GDALGetRasterYSize(hDstDS));
|
|
465 |
} else {
|
|
466 |
result = oWO.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(hDstDS), GDALGetRasterYSize(hDstDS));
|
|
467 |
}
|
|
468 |
} else
|
|
469 |
result = CE_Failure;
|
|
470 |
|
|
471 |
/* -------------------------------------------------------------------- */
|
|
472 |
/* Cleanup */
|
|
473 |
/* -------------------------------------------------------------------- */
|
|
474 |
if (hApproxArg != NULL)
|
|
475 |
GDALDestroyApproxTransformer(hApproxArg);
|
|
476 |
|
|
477 |
if (hGenImgProjArg != NULL)
|
|
478 |
GDALDestroyGenImgProjTransformer(hGenImgProjArg);
|
|
479 |
|
|
480 |
GDALDestroyWarpOptions(psWO);
|
|
481 |
|
|
482 |
GDALClose(hSrcDS);
|
|
483 |
}
|
|
484 |
|
|
485 |
/* -------------------------------------------------------------------- */
|
|
486 |
/* Final Cleanup. */
|
|
487 |
/* -------------------------------------------------------------------- */
|
|
488 |
GDALClose(hDstDS);
|
|
489 |
|
|
490 |
CPLFree(pszTargetSRS);
|
|
491 |
CPLFree(pszDstFilename);
|
|
492 |
CSLDestroy(papszSrcFiles);
|
|
493 |
CSLDestroy(papszWarpOptions);
|
|
494 |
|
|
495 |
GDALDestroyDriverManager();
|
|
496 |
|
|
497 |
return result;
|
|
498 |
}
|
|
499 |
|
|
500 |
|
|
501 |
|
|
502 |
|
|
503 |
/************************************************************************/
|
|
504 |
/* GDALWarpCreateOutput() */
|
|
505 |
/* */
|
|
506 |
/* Create the output file based on various commandline options, */
|
|
507 |
/* and the input file. */
|
|
508 |
/************************************************************************/
|
|
509 |
|
|
510 |
static GDALDatasetH
|
|
511 |
GDALWarpCreateOutput( char **papszSrcFiles, const char *pszFilename,
|
|
512 |
const char *pszFormat, const char *pszSourceSRS,
|
|
513 |
const char *pszTargetSRS, int nOrder,
|
|
514 |
char **papszCreateOptions, GDALDataType eDT )
|
|
515 |
|
|
516 |
|
|
517 |
{
|
|
518 |
GDALDriverH hDriver;
|
|
519 |
GDALDatasetH hDstDS;
|
|
520 |
void *hTransformArg;
|
|
521 |
GDALColorTableH hCT = NULL;
|
|
522 |
double dfWrkMinX=0, dfWrkMaxX=0, dfWrkMinY=0, dfWrkMaxY=0;
|
|
523 |
double dfWrkResX=0, dfWrkResY=0;
|
|
524 |
int nDstBandCount = 0;
|
|
525 |
|
|
526 |
/* -------------------------------------------------------------------- */
|
|
527 |
/* Find the output driver. */
|
|
528 |
/* -------------------------------------------------------------------- */
|
|
529 |
hDriver = GDALGetDriverByName( pszFormat );
|
|
530 |
if( hDriver == NULL
|
|
531 |
|| GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL )
|
|
532 |
{
|
|
533 |
int iDr;
|
|
534 |
|
|
535 |
printf( "Output driver `%s' not recognised or does not support\n",
|
|
536 |
pszFormat );
|
|
537 |
printf( "direct output file creation. The following format drivers are configured\n"
|
|
538 |
"and support direct output:\n" );
|
|
539 |
|
|
540 |
for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
|
|
541 |
{
|
|
542 |
GDALDriverH hDriver = GDALGetDriver(iDr);
|
|
543 |
|
|
544 |
if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL )
|
|
545 |
{
|
|
546 |
printf( " %s: %s\n",
|
|
547 |
GDALGetDriverShortName( hDriver ),
|
|
548 |
GDALGetDriverLongName( hDriver ) );
|
|
549 |
}
|
|
550 |
}
|
|
551 |
printf( "\n" );
|
|
552 |
exit( 1 );
|
|
553 |
}
|
|
554 |
|
|
555 |
/* -------------------------------------------------------------------- */
|
|
556 |
/* For virtual output files, we have to set a special subclass */
|
|
557 |
/* of dataset to create. */
|
|
558 |
/* -------------------------------------------------------------------- */
|
|
559 |
if( bVRT )
|
|
560 |
papszCreateOptions =
|
|
561 |
CSLSetNameValue( papszCreateOptions, "SUBCLASS",
|
|
562 |
"VRTWarpedDataset" );
|
|
563 |
|
|
564 |
/* -------------------------------------------------------------------- */
|
|
565 |
/* Loop over all input files to collect extents. */
|
|
566 |
/* -------------------------------------------------------------------- */
|
|
567 |
int iSrc;
|
|
568 |
|
|
569 |
for( iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ )
|
|
570 |
{
|
|
571 |
GDALDatasetH hSrcDS;
|
|
572 |
const char *pszThisSourceSRS = pszSourceSRS;
|
|
573 |
|
|
574 |
hSrcDS = GDALOpen( papszSrcFiles[iSrc], GA_ReadOnly );
|
|
575 |
if( hSrcDS == NULL )
|
|
576 |
exit( 1 );
|
|
577 |
|
|
578 |
if( eDT == GDT_Unknown )
|
|
579 |
eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));
|
|
580 |
|
|
581 |
/* -------------------------------------------------------------------- */
|
|
582 |
/* If we are processing the first file, and it has a color */
|
|
583 |
/* table, then we will copy it to the destination file. */
|
|
584 |
/* -------------------------------------------------------------------- */
|
|
585 |
if( iSrc == 0 )
|
|
586 |
{
|
|
587 |
nDstBandCount = GDALGetRasterCount(hSrcDS);
|
|
588 |
hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
|
|
589 |
if( hCT != NULL )
|
|
590 |
{
|
|
591 |
hCT = GDALCloneColorTable( hCT );
|
|
592 |
printf( "Copying color table from %s to new file.\n",
|
|
593 |
papszSrcFiles[iSrc] );
|
|
594 |
}
|
|
595 |
}
|
|
596 |
|
|
597 |
/* -------------------------------------------------------------------- */
|
|
598 |
/* Get the sourcesrs from the dataset, if not set already. */
|
|
599 |
/* -------------------------------------------------------------------- */
|
|
600 |
if( pszThisSourceSRS == NULL )
|
|
601 |
{
|
|
602 |
if( GDALGetProjectionRef( hSrcDS ) != NULL
|
|
603 |
&& strlen(GDALGetProjectionRef( hSrcDS )) > 0 )
|
|
604 |
pszThisSourceSRS = GDALGetProjectionRef( hSrcDS );
|
|
605 |
|
|
606 |
else if( GDALGetGCPProjection( hSrcDS ) != NULL
|
|
607 |
&& strlen(GDALGetGCPProjection(hSrcDS)) > 0
|
|
608 |
&& GDALGetGCPCount( hSrcDS ) > 1 )
|
|
609 |
pszThisSourceSRS = GDALGetGCPProjection( hSrcDS );
|
|
610 |
else
|
|
611 |
pszThisSourceSRS = "";
|
|
612 |
}
|
|
613 |
|
|
614 |
if( pszTargetSRS == NULL )
|
|
615 |
pszTargetSRS = pszThisSourceSRS;
|
|
616 |
|
|
617 |
/* -------------------------------------------------------------------- */
|
|
618 |
/* Create a transformation object from the source to */
|
|
619 |
/* destination coordinate system. */
|
|
620 |
/* -------------------------------------------------------------------- */
|
|
621 |
hTransformArg =
|
|
622 |
GDALCreateGenImgProjTransformer( hSrcDS, pszThisSourceSRS,
|
|
623 |
NULL, pszTargetSRS,
|
|
624 |
TRUE, 1000.0, nOrder );
|
|
625 |
|
|
626 |
if( hTransformArg == NULL )
|
|
627 |
return NULL;
|
|
628 |
|
|
629 |
/* -------------------------------------------------------------------- */
|
|
630 |
/* Get approximate output definition. */
|
|
631 |
/* -------------------------------------------------------------------- */
|
|
632 |
double adfThisGeoTransform[6];
|
|
633 |
double adfExtent[4];
|
|
634 |
int nThisPixels, nThisLines;
|
|
635 |
|
|
636 |
if( GDALSuggestedWarpOutput2( hSrcDS,
|
|
637 |
GDALGenImgProjTransform, hTransformArg,
|
|
638 |
adfThisGeoTransform,
|
|
639 |
&nThisPixels, &nThisLines,
|
|
640 |
adfExtent, 0 ) != CE_None )
|
|
641 |
return NULL;
|
|
642 |
|
|
643 |
/* -------------------------------------------------------------------- */
|
|
644 |
/* Expand the working bounds to include this region, ensure the */
|
|
645 |
/* working resolution is no more than this resolution. */
|
|
646 |
/* -------------------------------------------------------------------- */
|
|
647 |
if( dfWrkMaxX == 0.0 && dfWrkMinX == 0.0 )
|
|
648 |
{
|
|
649 |
dfWrkMinX = adfExtent[0];
|
|
650 |
dfWrkMaxX = adfExtent[2];
|
|
651 |
dfWrkMaxY = adfExtent[3];
|
|
652 |
dfWrkMinY = adfExtent[1];
|
|
653 |
dfWrkResX = adfThisGeoTransform[1];
|
|
654 |
dfWrkResY = ABS(adfThisGeoTransform[5]);
|
|
655 |
}
|
|
656 |
else
|
|
657 |
{
|
|
658 |
dfWrkMinX = MIN(dfWrkMinX,adfExtent[0]);
|
|
659 |
dfWrkMaxX = MAX(dfWrkMaxX,adfExtent[2]);
|
|
660 |
dfWrkMaxY = MAX(dfWrkMaxY,adfExtent[3]);
|
|
661 |
dfWrkMinY = MIN(dfWrkMinY,adfExtent[1]);
|
|
662 |
dfWrkResX = MIN(dfWrkResX,adfThisGeoTransform[1]);
|
|
663 |
dfWrkResY = MIN(dfWrkResY,ABS(adfThisGeoTransform[5]));
|
|
664 |
}
|
|
665 |
|
|
666 |
GDALDestroyGenImgProjTransformer( hTransformArg );
|
|
667 |
|
|
668 |
GDALClose( hSrcDS );
|
|
669 |
}
|
|
670 |
|
|
671 |
/* -------------------------------------------------------------------- */
|
|
672 |
/* Did we have any usable sources? */
|
|
673 |
/* -------------------------------------------------------------------- */
|
|
674 |
if( nDstBandCount == 0 )
|
|
675 |
{
|
|
676 |
CPLError( CE_Failure, CPLE_AppDefined,
|
|
677 |
"No usable source images." );
|
|
678 |
return NULL;
|
|
679 |
}
|
|
680 |
|
|
681 |
/* -------------------------------------------------------------------- */
|
|
682 |
/* Turn the suggested region into a geotransform and suggested */
|
|
683 |
/* number of pixels and lines. */
|
|
684 |
/* -------------------------------------------------------------------- */
|
|
685 |
double adfDstGeoTransform[6];
|
|
686 |
int nPixels, nLines;
|
|
687 |
|
|
688 |
adfDstGeoTransform[0] = dfWrkMinX;
|
|
689 |
adfDstGeoTransform[1] = dfWrkResX;
|
|
690 |
adfDstGeoTransform[2] = 0.0;
|
|
691 |
adfDstGeoTransform[3] = dfWrkMaxY;
|
|
692 |
adfDstGeoTransform[4] = 0.0;
|
|
693 |
adfDstGeoTransform[5] = -1 * dfWrkResY;
|
|
694 |
|
|
695 |
nPixels = (int) ((dfWrkMaxX - dfWrkMinX) / dfWrkResX + 0.5);
|
|
696 |
nLines = (int) ((dfWrkMaxY - dfWrkMinY) / dfWrkResY + 0.5);
|
|
697 |
|
|
698 |
/* -------------------------------------------------------------------- */
|
|
699 |
/* Did the user override some parameters? */
|
|
700 |
/* -------------------------------------------------------------------- */
|
|
701 |
if( dfXRes != 0.0 && dfYRes != 0.0 )
|
|
702 |
{
|
|
703 |
CPLAssert( nForcePixels == 0 && nForceLines == 0 );
|
|
704 |
if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
|
|
705 |
{
|
|
706 |
dfMinX = adfDstGeoTransform[0];
|
|
707 |
dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
|
|
708 |
dfMaxY = adfDstGeoTransform[3];
|
|
709 |
dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
|
|
710 |
}
|
|
711 |
|
|
712 |
nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
|
|
713 |
nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
|
|
714 |
adfDstGeoTransform[0] = dfMinX;
|
|
715 |
adfDstGeoTransform[3] = dfMaxY;
|
|
716 |
adfDstGeoTransform[1] = dfXRes;
|
|
717 |
adfDstGeoTransform[5] = -dfYRes;
|
|
718 |
}
|
|
719 |
|
|
720 |
else if( nForcePixels != 0 && nForceLines != 0 )
|
|
721 |
{
|
|
722 |
if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
|
|
723 |
{
|
|
724 |
dfMinX = dfWrkMinX;
|
|
725 |
dfMaxX = dfWrkMaxX;
|
|
726 |
dfMaxY = dfWrkMaxY;
|
|
727 |
dfMinY = dfWrkMinY;
|
|
728 |
}
|
|
729 |
|
|
730 |
dfXRes = (dfMaxX - dfMinX) / nForcePixels;
|
|
731 |
dfYRes = (dfMaxY - dfMinY) / nForceLines;
|
|
732 |
|
|
733 |
adfDstGeoTransform[0] = dfMinX;
|
|
734 |
adfDstGeoTransform[3] = dfMaxY;
|
|
735 |
adfDstGeoTransform[1] = dfXRes;
|
|
736 |
adfDstGeoTransform[5] = -dfYRes;
|
|
737 |
|
|
738 |
nPixels = nForcePixels;
|
|
739 |
nLines = nForceLines;
|
|
740 |
}
|
|
741 |
|
|
742 |
else if( nForcePixels != 0 )
|
|
743 |
{
|
|
744 |
if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
|
|
745 |
{
|
|
746 |
dfMinX = dfWrkMinX;
|
|
747 |
dfMaxX = dfWrkMaxX;
|
|
748 |
dfMaxY = dfWrkMaxY;
|
|
749 |
dfMinY = dfWrkMinY;
|
|
750 |
}
|
|
751 |
|
|
752 |
dfXRes = (dfMaxX - dfMinX) / nForcePixels;
|
|
753 |
dfYRes = dfXRes;
|
|
754 |
|
|
755 |
adfDstGeoTransform[0] = dfMinX;
|
|
756 |
adfDstGeoTransform[3] = dfMaxY;
|
|
757 |
adfDstGeoTransform[1] = dfXRes;
|
|
758 |
adfDstGeoTransform[5] = -dfYRes;
|
|
759 |
|
|
760 |
nPixels = nForcePixels;
|
|
761 |
nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
|
|
762 |
}
|
|
763 |
|
|
764 |
else if( nForceLines != 0 )
|
|
765 |
{
|
|
766 |
if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
|
|
767 |
{
|
|
768 |
dfMinX = dfWrkMinX;
|
|
769 |
dfMaxX = dfWrkMaxX;
|
|
770 |
dfMaxY = dfWrkMaxY;
|
|
771 |
dfMinY = dfWrkMinY;
|
|
772 |
}
|
|
773 |
|
|
774 |
dfYRes = (dfMaxY - dfMinY) / nForceLines;
|
|
775 |
dfXRes = dfYRes;
|
|
776 |
|
|
777 |
adfDstGeoTransform[0] = dfMinX;
|
|
778 |
adfDstGeoTransform[3] = dfMaxY;
|
|
779 |
adfDstGeoTransform[1] = dfXRes;
|
|
780 |
adfDstGeoTransform[5] = -dfYRes;
|
|
781 |
|
|
782 |
nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
|
|
783 |
nLines = nForceLines;
|
|
784 |
}
|
|
785 |
|
|
786 |
else if( dfMinX != 0.0 || dfMinY != 0.0 || dfMaxX != 0.0 || dfMaxY != 0.0 )
|
|
787 |
{
|
|
788 |
dfXRes = adfDstGeoTransform[1];
|
|
789 |
dfYRes = fabs(adfDstGeoTransform[5]);
|
|
790 |
|
|
791 |
nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
|
|
792 |
nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
|
|
793 |
|
|
794 |
adfDstGeoTransform[0] = dfMinX;
|
|
795 |
adfDstGeoTransform[3] = dfMaxY;
|
|
796 |
}
|
|
797 |
|
|
798 |
/* -------------------------------------------------------------------- */
|
|
799 |
/* Do we want to generate an alpha band in the output file? */
|
|
800 |
/* -------------------------------------------------------------------- */
|
|
801 |
if( bEnableSrcAlpha )
|
|
802 |
nDstBandCount--;
|
|
803 |
|
|
804 |
if( bEnableDstAlpha )
|
|
805 |
nDstBandCount++;
|
|
806 |
|
|
807 |
/* -------------------------------------------------------------------- */
|
|
808 |
/* Create the output file. */
|
|
809 |
/* -------------------------------------------------------------------- */
|
|
810 |
if( !bQuiet )
|
|
811 |
printf( "Creating output file that is %dP x %dL.\n", nPixels, nLines );
|
|
812 |
|
|
813 |
hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines,
|
|
814 |
nDstBandCount, eDT, papszCreateOptions );
|
|
815 |
|
|
816 |
if( hDstDS == NULL )
|
|
817 |
return NULL;
|
|
818 |
|
|
819 |
/* -------------------------------------------------------------------- */
|
|
820 |
/* Write out the projection definition. */
|
|
821 |
/* -------------------------------------------------------------------- */
|
|
822 |
GDALSetProjection( hDstDS, pszTargetSRS );
|
|
823 |
GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
|
|
824 |
|
|
825 |
/* -------------------------------------------------------------------- */
|
|
826 |
/* Try to set color interpretation of output file alpha band. */
|
|
827 |
/* TODO: We should likely try to copy the other bands too. */
|
|
828 |
/* -------------------------------------------------------------------- */
|
|
829 |
if( bEnableDstAlpha )
|
|
830 |
{
|
|
831 |
GDALSetRasterColorInterpretation(
|
|
832 |
GDALGetRasterBand( hDstDS, nDstBandCount ),
|
|
833 |
GCI_AlphaBand );
|
|
834 |
}
|
|
835 |
|
|
836 |
/* -------------------------------------------------------------------- */
|
|
837 |
/* Copy the color table, if required. */
|
|
838 |
/* -------------------------------------------------------------------- */
|
|
839 |
if( hCT != NULL )
|
|
840 |
{
|
|
841 |
GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
|
|
842 |
GDALDestroyColorTable( hCT );
|
|
843 |
}
|
|
844 |
|
|
845 |
return hDstDS;
|
|
846 |
}
|
|
847 |
|
|
848 |
/************************************************************************/
|
|
849 |
/* InsertCenterLong() */
|
|
850 |
/* */
|
|
851 |
/* Insert a CENTER_LONG Extension entry on a GEOGCS to indicate */
|
|
852 |
/* the center longitude of the dataset for wrapping purposes. */
|
|
853 |
/************************************************************************/
|
|
854 |
|
|
855 |
static CPLString InsertCenterLong( GDALDatasetH hDS, CPLString osWKT )
|
|
856 |
|
|
857 |
{
|
|
858 |
if( !EQUALN(osWKT.c_str(), "GEOGCS[", 7) )
|
|
859 |
return osWKT;
|
|
860 |
|
|
861 |
if( strstr(osWKT,"EXTENSION[\"CENTER_LONG") != NULL )
|
|
862 |
return osWKT;
|
|
863 |
|
|
864 |
/* -------------------------------------------------------------------- */
|
|
865 |
/* For now we only do this if we have a geotransform since */
|
|
866 |
/* other forms require a bunch of extra work. */
|
|
867 |
/* -------------------------------------------------------------------- */
|
|
868 |
double adfGeoTransform[6];
|
|
869 |
|
|
870 |
if( GDALGetGeoTransform( hDS, adfGeoTransform ) != CE_None )
|
|
871 |
return osWKT;
|
|
872 |
|
|
873 |
/* -------------------------------------------------------------------- */
|
|
874 |
/* Compute min/max longitude based on testing the four corners. */
|
|
875 |
/* -------------------------------------------------------------------- */
|
|
876 |
double dfMinLong, dfMaxLong;
|
|
877 |
int nXSize = GDALGetRasterXSize( hDS );
|
|
878 |
int nYSize = GDALGetRasterYSize( hDS );
|
|
879 |
|
|
880 |
dfMinLong =
|
|
881 |
MIN(MIN(adfGeoTransform[0] + 0 * adfGeoTransform[1]
|
|
882 |
+ 0 * adfGeoTransform[2],
|
|
883 |
adfGeoTransform[0] + nXSize * adfGeoTransform[1]
|
|
884 |
+ 0 * adfGeoTransform[2]),
|
|
885 |
MIN(adfGeoTransform[0] + 0 * adfGeoTransform[1]
|
|
886 |
+ nYSize * adfGeoTransform[2],
|
|
887 |
adfGeoTransform[0] + nXSize * adfGeoTransform[1]
|
|
888 |
+ nYSize * adfGeoTransform[2]));
|
|
889 |
dfMaxLong =
|
|
890 |
MAX(MAX(adfGeoTransform[0] + 0 * adfGeoTransform[1]
|
|
891 |
+ 0 * adfGeoTransform[2],
|
|
892 |
adfGeoTransform[0] + nXSize * adfGeoTransform[1]
|
|
893 |
+ 0 * adfGeoTransform[2]),
|
|
894 |
MAX(adfGeoTransform[0] + 0 * adfGeoTransform[1]
|
|
895 |
+ nYSize * adfGeoTransform[2],
|
|
896 |
adfGeoTransform[0] + nXSize * adfGeoTransform[1]
|
|
897 |
+ nYSize * adfGeoTransform[2]));
|
|
898 |
|
|
899 |
if( dfMaxLong - dfMinLong > 360.0 )
|
|
900 |
return osWKT;
|
|
901 |
|
|
902 |
/* -------------------------------------------------------------------- */
|
|
903 |
/* Insert center long. */
|
|
904 |
/* -------------------------------------------------------------------- */
|
|
905 |
OGRSpatialReference oSRS( osWKT );
|
|
906 |
double dfCenterLong = (dfMaxLong + dfMinLong) / 2.0;
|
|
907 |
OGR_SRSNode *poExt;
|
|
908 |
|
|
909 |
poExt = new OGR_SRSNode( "EXTENSION" );
|
|
910 |
poExt->AddChild( new OGR_SRSNode( "CENTER_LONG" ) );
|
|
911 |
poExt->AddChild( new OGR_SRSNode( CPLString().Printf("%g",dfCenterLong) ));
|
|
912 |
|
|
913 |
oSRS.GetRoot()->AddChild( poExt );
|
|
914 |
|
|
915 |
/* -------------------------------------------------------------------- */
|
|
916 |
/* Convert back to wkt. */
|
|
917 |
/* -------------------------------------------------------------------- */
|
|
918 |
char *pszWKT = NULL;
|
|
919 |
oSRS.exportToWkt( &pszWKT );
|
|
920 |
|
|
921 |
osWKT = pszWKT;
|
|
922 |
CPLFree( pszWKT );
|
|
923 |
|
|
924 |
return osWKT;
|
|
925 |
}
|
|
926 |
|
|
927 |
|
|
928 |
/************************************************************************/
|
|
929 |
/* GDALTermProgress() */
|
|
930 |
/************************************************************************/
|
|
931 |
|
|
932 |
/**
|
|
933 |
* Simple progress report to terminal.
|
|
934 |
*
|
|
935 |
* This progress reporter prints simple progress report to the
|
|
936 |
* terminal window. The progress report generally looks something like
|
|
937 |
* this:
|
|
938 |
|
|
939 |
\verbatim
|
|
940 |
0...10...20...30...40...50...60...70...80...90...100 - done.
|
|
941 |
\endverbatim
|
|
942 |
|
|
943 |
* Every 2.5% of progress another number or period is emitted. Note that
|
|
944 |
* GDALTermProgress() uses internal static data to keep track of the last
|
|
945 |
* percentage reported and will get confused if two terminal based progress
|
|
946 |
* reportings are active at the same time.
|
|
947 |
*
|
|
948 |
* The GDALTermProgress() function maintains an internal memory of the
|
|
949 |
* last percentage complete reported in a static variable, and this makes
|
|
950 |
* it unsuitable to have multiple GDALTermProgress()'s active eithin a
|
|
951 |
* single thread or across multiple threads.
|
|
952 |
*
|
|
953 |
* @param dfComplete completion ratio from 0.0 to 1.0.
|
|
954 |
* @param pszMessage optional message.
|
|
955 |
* @param pProgressArg ignored callback data argument.
|
|
956 |
*
|
|
957 |
* @return Always returns TRUE indicating the process should continue.
|
|
958 |
*/
|
|
959 |
|
|
960 |
int CPL_STDCALL GDALFuncTermProgress( double dfComplete, const char *pszMessage,
|
|
961 |
void * pProgressArg )
|
|
962 |
|
|
963 |
{
|
|
964 |
static double dfLastComplete = -1.0;
|
|
965 |
|
|
966 |
(void) pProgressArg;
|
|
967 |
|
|
968 |
int nPercent = (int) floor(dfComplete*100);
|
|
969 |
|
|
970 |
statusCallBack(nPercent);
|
|
971 |
fflush( stdout );
|
|
972 |
|
|
973 |
dfLastComplete = dfComplete;
|
|
974 |
|
|
975 |
return TRUE;
|
|
976 |
}
|