1 : /******************************************************************************
2 : * $Id: gdaljp2metadata.cpp 24228 2012-04-13 21:12:40Z rouault $
3 : *
4 : * Project: GDAL
5 : * Purpose: GDALJP2Metadata - Read GeoTIFF and/or GML georef info.
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 :
30 : #include "gdaljp2metadata.h"
31 : #include "cpl_string.h"
32 : #include "cpl_minixml.h"
33 : #include "ogr_spatialref.h"
34 : #include "ogr_geometry.h"
35 : #include "ogr_api.h"
36 : #include "gt_wkt_srs_for_gdal.h"
37 :
38 : CPL_CVSID("$Id: gdaljp2metadata.cpp 24228 2012-04-13 21:12:40Z rouault $");
39 :
40 : static const unsigned char msi_uuid2[16] =
41 : {0xb1,0x4b,0xf8,0xbd,0x08,0x3d,0x4b,0x43,
42 : 0xa5,0xae,0x8c,0xd7,0xd5,0xa6,0xce,0x03};
43 :
44 : static const unsigned char msig_uuid[16] =
45 : { 0x96,0xA9,0xF1,0xF1,0xDC,0x98,0x40,0x2D,
46 : 0xA7,0xAE,0xD6,0x8E,0x34,0x45,0x18,0x09 };
47 :
48 : static const unsigned char xmp_uuid[16] =
49 : { 0xBE,0x7A,0xCF,0xCB,0x97,0xA9,0x42,0xE8,
50 : 0x9C,0x71,0x99,0x94,0x91,0xE3,0xAF,0xAC};
51 :
52 : /************************************************************************/
53 : /* GDALJP2Metadata() */
54 : /************************************************************************/
55 :
56 336 : GDALJP2Metadata::GDALJP2Metadata()
57 :
58 : {
59 336 : pszProjection = NULL;
60 :
61 336 : nGCPCount = 0;
62 336 : pasGCPList = NULL;
63 :
64 336 : papszGMLMetadata = NULL;
65 336 : papszMetadata = NULL;
66 :
67 336 : nGeoTIFFSize = 0;
68 336 : pabyGeoTIFFData = NULL;
69 :
70 336 : nMSIGSize = 0;
71 336 : pabyMSIGData = NULL;
72 :
73 336 : pszXMPMetadata = NULL;
74 :
75 336 : bHaveGeoTransform = FALSE;
76 336 : adfGeoTransform[0] = 0.0;
77 336 : adfGeoTransform[1] = 1.0;
78 336 : adfGeoTransform[2] = 0.0;
79 336 : adfGeoTransform[3] = 0.0;
80 336 : adfGeoTransform[4] = 0.0;
81 336 : adfGeoTransform[5] = 1.0;
82 336 : }
83 :
84 : /************************************************************************/
85 : /* ~GDALJP2Metadata() */
86 : /************************************************************************/
87 :
88 336 : GDALJP2Metadata::~GDALJP2Metadata()
89 :
90 : {
91 336 : CPLFree( pszProjection );
92 336 : if( nGCPCount > 0 )
93 : {
94 14 : GDALDeinitGCPs( nGCPCount, pasGCPList );
95 14 : CPLFree( pasGCPList );
96 : }
97 :
98 336 : CPLFree( pabyGeoTIFFData );
99 336 : CPLFree( pabyMSIGData );
100 336 : CSLDestroy( papszGMLMetadata );
101 336 : CSLDestroy( papszMetadata );
102 336 : CPLFree( pszXMPMetadata );
103 336 : }
104 :
105 : /************************************************************************/
106 : /* ReadAndParse() */
107 : /* */
108 : /* Read a JP2 file and try to collect georeferencing */
109 : /* information from the various available forms. Returns TRUE */
110 : /* if anything useful is found. */
111 : /************************************************************************/
112 :
113 262 : int GDALJP2Metadata::ReadAndParse( const char *pszFilename )
114 :
115 : {
116 : VSILFILE *fpLL;
117 :
118 262 : fpLL = VSIFOpenL( pszFilename, "rb" );
119 :
120 262 : if( fpLL == NULL )
121 : {
122 : CPLDebug( "GDALJP2Metadata", "Could not even open %s.",
123 0 : pszFilename );
124 :
125 0 : return FALSE;
126 : }
127 :
128 262 : ReadBoxes( fpLL );
129 262 : VSIFCloseL( fpLL );
130 :
131 : /* -------------------------------------------------------------------- */
132 : /* Try JP2GeoTIFF, GML and finally MSIG to get something. */
133 : /* -------------------------------------------------------------------- */
134 262 : if( !ParseJP2GeoTIFF() && !ParseGMLCoverageDesc() )
135 116 : ParseMSIG();
136 :
137 : /* -------------------------------------------------------------------- */
138 : /* If we still don't have a geotransform, look for a world */
139 : /* file. */
140 : /* -------------------------------------------------------------------- */
141 262 : if( !bHaveGeoTransform )
142 : {
143 : bHaveGeoTransform =
144 : GDALReadWorldFile( pszFilename, NULL, adfGeoTransform )
145 136 : || GDALReadWorldFile( pszFilename, ".wld", adfGeoTransform );
146 : }
147 :
148 : /* -------------------------------------------------------------------- */
149 : /* Return success either either of projection or geotransform */
150 : /* or gcps. */
151 : /* -------------------------------------------------------------------- */
152 : return bHaveGeoTransform
153 : || nGCPCount > 0
154 262 : || (pszProjection != NULL && strlen(pszProjection) > 0);
155 : }
156 :
157 : /************************************************************************/
158 : /* CollectGMLData() */
159 : /* */
160 : /* Read all the asoc boxes after this node, and store the */
161 : /* contain xml documents along with the name from the label. */
162 : /************************************************************************/
163 :
164 56 : void GDALJP2Metadata::CollectGMLData( GDALJP2Box *poGMLData )
165 :
166 : {
167 56 : GDALJP2Box oChildBox( poGMLData->GetFILE() );
168 :
169 56 : oChildBox.ReadFirstChild( poGMLData );
170 :
171 224 : while( strlen(oChildBox.GetType()) > 0 )
172 : {
173 112 : if( EQUAL(oChildBox.GetType(),"asoc") )
174 : {
175 56 : GDALJP2Box oSubChildBox( oChildBox.GetFILE() );
176 :
177 56 : char *pszLabel = NULL;
178 56 : char *pszXML = NULL;
179 :
180 56 : oSubChildBox.ReadFirstChild( &oChildBox );
181 :
182 224 : while( strlen(oSubChildBox.GetType()) > 0 )
183 : {
184 112 : if( EQUAL(oSubChildBox.GetType(),"lbl ") )
185 56 : pszLabel = (char *)oSubChildBox.ReadBoxData();
186 56 : else if( EQUAL(oSubChildBox.GetType(),"xml ") )
187 56 : pszXML = (char *) oSubChildBox.ReadBoxData();
188 :
189 112 : oSubChildBox.ReadNextChild( &oChildBox );
190 : }
191 :
192 56 : if( pszLabel != NULL && pszXML != NULL )
193 : papszGMLMetadata = CSLSetNameValue( papszGMLMetadata,
194 56 : pszLabel, pszXML );
195 56 : CPLFree( pszLabel );
196 56 : CPLFree( pszXML );
197 : }
198 :
199 112 : oChildBox.ReadNextChild( poGMLData );
200 56 : }
201 56 : }
202 :
203 : /************************************************************************/
204 : /* ReadBoxes() */
205 : /************************************************************************/
206 :
207 262 : int GDALJP2Metadata::ReadBoxes( VSILFILE *fpVSIL )
208 :
209 : {
210 262 : GDALJP2Box oBox( fpVSIL );
211 262 : int iBox = 0;
212 :
213 262 : if (!oBox.ReadFirst())
214 0 : return FALSE;
215 :
216 1226 : while( strlen(oBox.GetType()) > 0 )
217 : {
218 : #ifdef DEBUG
219 866 : if (CSLTestBoolean(CPLGetConfigOption("DUMP_JP2_BOXES", "NO")))
220 0 : oBox.DumpReadable(stderr);
221 : #endif
222 :
223 : /* -------------------------------------------------------------------- */
224 : /* Collect geotiff box. */
225 : /* -------------------------------------------------------------------- */
226 866 : if( EQUAL(oBox.GetType(),"uuid")
227 : && memcmp( oBox.GetUUID(), msi_uuid2, 16 ) == 0 )
228 : {
229 138 : nGeoTIFFSize = (int) oBox.GetDataLength();
230 138 : pabyGeoTIFFData = oBox.ReadBoxData();
231 138 : if (pabyGeoTIFFData == NULL)
232 : {
233 0 : CPLDebug("GDALJP2", "Cannot read data for UUID GeoTIFF box");
234 0 : nGeoTIFFSize = 0;
235 : }
236 : }
237 :
238 : /* -------------------------------------------------------------------- */
239 : /* Collect MSIG box. */
240 : /* -------------------------------------------------------------------- */
241 866 : if( EQUAL(oBox.GetType(),"uuid")
242 : && memcmp( oBox.GetUUID(), msig_uuid, 16 ) == 0 )
243 : {
244 0 : nMSIGSize = (int) oBox.GetDataLength();
245 0 : pabyMSIGData = oBox.ReadBoxData();
246 :
247 0 : if( nMSIGSize < 70
248 : || pabyMSIGData == NULL
249 : || memcmp( pabyMSIGData, "MSIG/", 5 ) != 0 )
250 : {
251 0 : CPLFree( pabyMSIGData );
252 0 : pabyMSIGData = NULL;
253 0 : nMSIGSize = 0;
254 : }
255 : }
256 :
257 : /* -------------------------------------------------------------------- */
258 : /* Collect XMP box. */
259 : /* -------------------------------------------------------------------- */
260 866 : if( EQUAL(oBox.GetType(),"uuid")
261 : && memcmp( oBox.GetUUID(), xmp_uuid, 16 ) == 0 &&
262 : pszXMPMetadata == NULL )
263 : {
264 8 : pszXMPMetadata = (char*) oBox.ReadBoxData();
265 : }
266 :
267 : /* -------------------------------------------------------------------- */
268 : /* Process asoc box looking for Labelled GML data. */
269 : /* -------------------------------------------------------------------- */
270 866 : if( EQUAL(oBox.GetType(),"asoc") )
271 : {
272 56 : GDALJP2Box oSubBox( fpVSIL );
273 :
274 56 : oSubBox.ReadFirstChild( &oBox );
275 56 : if( EQUAL(oSubBox.GetType(),"lbl ") )
276 : {
277 56 : char *pszLabel = (char *) oSubBox.ReadBoxData();
278 56 : if( pszLabel != NULL && EQUAL(pszLabel,"gml.data") )
279 : {
280 56 : CollectGMLData( &oBox );
281 : }
282 56 : CPLFree( pszLabel );
283 56 : }
284 : }
285 :
286 : /* -------------------------------------------------------------------- */
287 : /* Process simple xml boxes. */
288 : /* -------------------------------------------------------------------- */
289 866 : if( EQUAL(oBox.GetType(),"xml ") )
290 : {
291 8 : CPLString osBoxName;
292 8 : char *pszXML = (char *) oBox.ReadBoxData();
293 :
294 8 : osBoxName.Printf( "BOX_%d", iBox++ );
295 :
296 : papszGMLMetadata = CSLSetNameValue( papszGMLMetadata,
297 8 : osBoxName, pszXML );
298 8 : CPLFree( pszXML );
299 : }
300 :
301 : /* -------------------------------------------------------------------- */
302 : /* Check for a resd box in jp2h. */
303 : /* -------------------------------------------------------------------- */
304 866 : if( EQUAL(oBox.GetType(),"jp2h") )
305 : {
306 164 : GDALJP2Box oSubBox( fpVSIL );
307 :
308 574 : for( oSubBox.ReadFirstChild( &oBox );
309 : strlen(oSubBox.GetType()) > 0;
310 : oSubBox.ReadNextChild( &oBox ) )
311 : {
312 410 : if( EQUAL(oSubBox.GetType(),"res ") )
313 : {
314 8 : GDALJP2Box oResBox( fpVSIL );
315 :
316 8 : oResBox.ReadFirstChild( &oSubBox );
317 :
318 : // we will use either the resd or resc box, which ever
319 : // happens to be first. Should we prefer resd?
320 8 : unsigned char *pabyResData = NULL;
321 8 : if( oResBox.GetDataLength() == 10 &&
322 : (pabyResData = oResBox.ReadBoxData()) != NULL )
323 : {
324 : int nVertNum, nVertDen, nVertExp;
325 : int nHorzNum, nHorzDen, nHorzExp;
326 :
327 8 : nVertNum = pabyResData[0] * 256 + pabyResData[1];
328 8 : nVertDen = pabyResData[2] * 256 + pabyResData[3];
329 8 : nHorzNum = pabyResData[4] * 256 + pabyResData[5];
330 8 : nHorzDen = pabyResData[6] * 256 + pabyResData[7];
331 8 : nVertExp = pabyResData[8];
332 8 : nHorzExp = pabyResData[9];
333 :
334 : // compute in pixels/cm
335 : double dfVertRes =
336 8 : (nVertNum/(double)nVertDen) * pow(10.0,nVertExp)/100;
337 : double dfHorzRes =
338 8 : (nHorzNum/(double)nHorzDen) * pow(10.0,nHorzExp)/100;
339 8 : CPLString osFormatter;
340 :
341 : papszMetadata = CSLSetNameValue(
342 : papszMetadata,
343 : "TIFFTAG_XRESOLUTION",
344 8 : osFormatter.Printf("%g",dfHorzRes) );
345 :
346 : papszMetadata = CSLSetNameValue(
347 : papszMetadata,
348 : "TIFFTAG_YRESOLUTION",
349 8 : osFormatter.Printf("%g",dfVertRes) );
350 : papszMetadata = CSLSetNameValue(
351 : papszMetadata,
352 : "TIFFTAG_RESOLUTIONUNIT",
353 8 : "3 (pixels/cm)" );
354 :
355 8 : CPLFree( pabyResData );
356 8 : }
357 : }
358 164 : }
359 : }
360 :
361 866 : if (!oBox.ReadNext())
362 164 : break;
363 : }
364 :
365 262 : return TRUE;
366 : }
367 :
368 : /************************************************************************/
369 : /* ParseJP2GeoTIFF() */
370 : /************************************************************************/
371 :
372 262 : int GDALJP2Metadata::ParseJP2GeoTIFF()
373 :
374 : {
375 262 : if( nGeoTIFFSize < 1 )
376 124 : return FALSE;
377 :
378 : /* -------------------------------------------------------------------- */
379 : /* Convert raw data into projection and geotransform. */
380 : /* -------------------------------------------------------------------- */
381 138 : int bSuccess = TRUE;
382 :
383 138 : if( GTIFWktFromMemBuf( nGeoTIFFSize, pabyGeoTIFFData,
384 : &pszProjection, adfGeoTransform,
385 : &nGCPCount, &pasGCPList ) != CE_None )
386 : {
387 0 : bSuccess = FALSE;
388 : }
389 :
390 138 : if( pszProjection == NULL || strlen(pszProjection) == 0 )
391 0 : bSuccess = FALSE;
392 :
393 138 : if( bSuccess )
394 : CPLDebug( "GDALJP2Metadata",
395 : "Got projection from GeoJP2 (geotiff) box: %s",
396 138 : pszProjection );
397 :
398 238 : if( adfGeoTransform[0] != 0
399 20 : || adfGeoTransform[1] != 1
400 20 : || adfGeoTransform[2] != 0
401 20 : || adfGeoTransform[3] != 0
402 20 : || adfGeoTransform[4] != 0
403 20 : || adfGeoTransform[5] != 1 )
404 118 : bHaveGeoTransform = TRUE;
405 :
406 138 : return bSuccess;;
407 : }
408 :
409 : /************************************************************************/
410 : /* ParseMSIG() */
411 : /************************************************************************/
412 :
413 116 : int GDALJP2Metadata::ParseMSIG()
414 :
415 : {
416 116 : if( nMSIGSize < 70 )
417 116 : return FALSE;
418 :
419 : /* -------------------------------------------------------------------- */
420 : /* Try and extract worldfile parameters and adjust. */
421 : /* -------------------------------------------------------------------- */
422 0 : memcpy( adfGeoTransform + 0, pabyMSIGData + 22 + 8 * 4, 8 );
423 0 : memcpy( adfGeoTransform + 1, pabyMSIGData + 22 + 8 * 0, 8 );
424 0 : memcpy( adfGeoTransform + 2, pabyMSIGData + 22 + 8 * 2, 8 );
425 0 : memcpy( adfGeoTransform + 3, pabyMSIGData + 22 + 8 * 5, 8 );
426 0 : memcpy( adfGeoTransform + 4, pabyMSIGData + 22 + 8 * 1, 8 );
427 0 : memcpy( adfGeoTransform + 5, pabyMSIGData + 22 + 8 * 3, 8 );
428 :
429 : // data is in LSB (little endian) order in file.
430 : CPL_LSBPTR64( adfGeoTransform + 0 );
431 : CPL_LSBPTR64( adfGeoTransform + 1 );
432 : CPL_LSBPTR64( adfGeoTransform + 2 );
433 : CPL_LSBPTR64( adfGeoTransform + 3 );
434 : CPL_LSBPTR64( adfGeoTransform + 4 );
435 : CPL_LSBPTR64( adfGeoTransform + 5 );
436 :
437 : // correct for center of pixel vs. top left of pixel
438 0 : adfGeoTransform[0] -= 0.5 * adfGeoTransform[1];
439 0 : adfGeoTransform[0] -= 0.5 * adfGeoTransform[2];
440 0 : adfGeoTransform[3] -= 0.5 * adfGeoTransform[4];
441 0 : adfGeoTransform[3] -= 0.5 * adfGeoTransform[5];
442 :
443 0 : bHaveGeoTransform = TRUE;
444 :
445 0 : return TRUE;
446 : }
447 :
448 : /************************************************************************/
449 : /* GetDictionaryItem() */
450 : /************************************************************************/
451 :
452 : static CPLXMLNode *
453 0 : GetDictionaryItem( char **papszGMLMetadata, const char *pszURN )
454 :
455 : {
456 : char *pszLabel;
457 0 : const char *pszFragmentId = NULL;
458 : int i;
459 :
460 :
461 0 : if( EQUALN(pszURN,"urn:jp2k:xml:", 13) )
462 0 : pszLabel = CPLStrdup( pszURN + 13 );
463 0 : else if( EQUALN(pszURN,"urn:ogc:tc:gmljp2:xml:", 22) )
464 0 : pszLabel = CPLStrdup( pszURN + 22 );
465 0 : else if( EQUALN(pszURN,"gmljp2://xml/",13) )
466 0 : pszLabel = CPLStrdup( pszURN + 13 );
467 : else
468 0 : pszLabel = CPLStrdup( pszURN );
469 :
470 : /* -------------------------------------------------------------------- */
471 : /* Split out label and fragment id. */
472 : /* -------------------------------------------------------------------- */
473 0 : for( i = 0; pszLabel[i] != '#'; i++ )
474 : {
475 0 : if( pszLabel[i] == '\0' )
476 0 : return NULL;
477 : }
478 :
479 0 : pszFragmentId = pszLabel + i + 1;
480 0 : pszLabel[i] = '\0';
481 :
482 : /* -------------------------------------------------------------------- */
483 : /* Can we find an XML box with the desired label? */
484 : /* -------------------------------------------------------------------- */
485 : const char *pszDictionary =
486 0 : CSLFetchNameValue( papszGMLMetadata, pszLabel );
487 :
488 0 : if( pszDictionary == NULL )
489 0 : return NULL;
490 :
491 : /* -------------------------------------------------------------------- */
492 : /* Try and parse the dictionary. */
493 : /* -------------------------------------------------------------------- */
494 0 : CPLXMLNode *psDictTree = CPLParseXMLString( pszDictionary );
495 :
496 0 : if( psDictTree == NULL )
497 : {
498 0 : CPLDestroyXMLNode( psDictTree );
499 0 : return NULL;
500 : }
501 :
502 0 : CPLStripXMLNamespace( psDictTree, NULL, TRUE );
503 :
504 0 : CPLXMLNode *psDictRoot = CPLSearchXMLNode( psDictTree, "=Dictionary" );
505 :
506 0 : if( psDictRoot == NULL )
507 : {
508 0 : CPLDestroyXMLNode( psDictTree );
509 0 : return NULL;
510 : }
511 :
512 : /* -------------------------------------------------------------------- */
513 : /* Search for matching id. */
514 : /* -------------------------------------------------------------------- */
515 0 : CPLXMLNode *psEntry, *psHit = NULL;
516 0 : for( psEntry = psDictRoot->psChild;
517 : psEntry != NULL && psHit == NULL;
518 : psEntry = psEntry->psNext )
519 : {
520 : const char *pszId;
521 :
522 0 : if( psEntry->eType != CXT_Element )
523 0 : continue;
524 :
525 0 : if( !EQUAL(psEntry->pszValue,"dictionaryEntry") )
526 0 : continue;
527 :
528 0 : if( psEntry->psChild == NULL )
529 0 : continue;
530 :
531 0 : pszId = CPLGetXMLValue( psEntry->psChild, "id", "" );
532 :
533 0 : if( EQUAL(pszId, pszFragmentId) )
534 0 : psHit = CPLCloneXMLTree( psEntry->psChild );
535 : }
536 :
537 : /* -------------------------------------------------------------------- */
538 : /* Cleanup */
539 : /* -------------------------------------------------------------------- */
540 0 : CPLFree( pszLabel );
541 0 : CPLDestroyXMLNode( psDictTree );
542 :
543 0 : return psHit;
544 : }
545 :
546 :
547 : /************************************************************************/
548 : /* GMLSRSLookup() */
549 : /* */
550 : /* Lookup an SRS in a dictionary inside this file. We will get */
551 : /* something like: */
552 : /* urn:jp2k:xml:CRSDictionary.xml#crs1112 */
553 : /* */
554 : /* We need to split the filename from the fragment id, and */
555 : /* lookup the fragment in the file if we can find it our */
556 : /* list of labelled xml boxes. */
557 : /************************************************************************/
558 :
559 0 : int GDALJP2Metadata::GMLSRSLookup( const char *pszURN )
560 :
561 : {
562 0 : CPLXMLNode *psDictEntry = GetDictionaryItem( papszGMLMetadata, pszURN );
563 :
564 0 : if( psDictEntry == NULL )
565 0 : return FALSE;
566 :
567 : /* -------------------------------------------------------------------- */
568 : /* Reserialize this fragment. */
569 : /* -------------------------------------------------------------------- */
570 0 : char *pszDictEntryXML = CPLSerializeXMLTree( psDictEntry );
571 0 : CPLDestroyXMLNode( psDictEntry );
572 :
573 : /* -------------------------------------------------------------------- */
574 : /* Try to convert into an OGRSpatialReference. */
575 : /* -------------------------------------------------------------------- */
576 0 : OGRSpatialReference oSRS;
577 0 : int bSuccess = FALSE;
578 :
579 0 : if( oSRS.importFromXML( pszDictEntryXML ) == OGRERR_NONE )
580 : {
581 0 : CPLFree( pszProjection );
582 0 : pszProjection = NULL;
583 :
584 0 : oSRS.exportToWkt( &pszProjection );
585 0 : bSuccess = TRUE;
586 : }
587 :
588 0 : CPLFree( pszDictEntryXML );
589 :
590 0 : return bSuccess;
591 : }
592 :
593 : /************************************************************************/
594 : /* ParseGMLCoverageDesc() */
595 : /************************************************************************/
596 :
597 124 : int GDALJP2Metadata::ParseGMLCoverageDesc()
598 :
599 : {
600 : /* -------------------------------------------------------------------- */
601 : /* Do we have an XML doc that is apparently a coverage */
602 : /* description? */
603 : /* -------------------------------------------------------------------- */
604 : const char *pszCoverage = CSLFetchNameValue( papszGMLMetadata,
605 124 : "gml.root-instance" );
606 :
607 124 : if( pszCoverage == NULL )
608 116 : return FALSE;
609 :
610 8 : CPLDebug( "GDALJP2Metadata", "Found GML Box:\n%s", pszCoverage );
611 :
612 : /* -------------------------------------------------------------------- */
613 : /* Try parsing the XML. Wipe any namespace prefixes. */
614 : /* -------------------------------------------------------------------- */
615 8 : CPLXMLNode *psXML = CPLParseXMLString( pszCoverage );
616 :
617 8 : if( psXML == NULL )
618 0 : return FALSE;
619 :
620 8 : CPLStripXMLNamespace( psXML, NULL, TRUE );
621 :
622 : /* -------------------------------------------------------------------- */
623 : /* Isolate RectifiedGrid. Eventually we will need to support */
624 : /* other georeferencing objects. */
625 : /* -------------------------------------------------------------------- */
626 8 : CPLXMLNode *psRG = CPLSearchXMLNode( psXML, "=RectifiedGrid" );
627 8 : CPLXMLNode *psOriginPoint = NULL;
628 8 : const char *pszOffset1=NULL, *pszOffset2=NULL;
629 :
630 8 : if( psRG != NULL )
631 : {
632 8 : psOriginPoint = CPLGetXMLNode( psRG, "origin.Point" );
633 :
634 :
635 8 : CPLXMLNode *psOffset1 = CPLGetXMLNode( psRG, "offsetVector" );
636 8 : if( psOffset1 != NULL )
637 : {
638 8 : pszOffset1 = CPLGetXMLValue( psOffset1, "", NULL );
639 : pszOffset2 = CPLGetXMLValue( psOffset1->psNext, "=offsetVector",
640 8 : NULL );
641 : }
642 : }
643 :
644 : /* -------------------------------------------------------------------- */
645 : /* If we are missing any of the origin or 2 offsets then give up. */
646 : /* -------------------------------------------------------------------- */
647 8 : if( psOriginPoint == NULL || pszOffset1 == NULL || pszOffset2 == NULL )
648 : {
649 0 : CPLDestroyXMLNode( psXML );
650 0 : return FALSE;
651 : }
652 :
653 : /* -------------------------------------------------------------------- */
654 : /* Extract origin location. */
655 : /* -------------------------------------------------------------------- */
656 8 : OGRPoint *poOriginGeometry = NULL;
657 8 : const char *pszSRSName = NULL;
658 :
659 8 : if( psOriginPoint != NULL )
660 : {
661 : poOriginGeometry = (OGRPoint *)
662 8 : OGR_G_CreateFromGMLTree( psOriginPoint );
663 :
664 16 : if( poOriginGeometry != NULL
665 8 : && wkbFlatten(poOriginGeometry->getGeometryType()) != wkbPoint )
666 : {
667 0 : delete poOriginGeometry;
668 0 : poOriginGeometry = NULL;
669 : }
670 :
671 : // SRS?
672 8 : pszSRSName = CPLGetXMLValue( psOriginPoint, "srsName", NULL );
673 : }
674 :
675 : /* -------------------------------------------------------------------- */
676 : /* Extract offset(s) */
677 : /* -------------------------------------------------------------------- */
678 8 : char **papszOffset1Tokens = NULL;
679 8 : char **papszOffset2Tokens = NULL;
680 8 : int bSuccess = FALSE;
681 :
682 : papszOffset1Tokens =
683 8 : CSLTokenizeStringComplex( pszOffset1, " ,", FALSE, FALSE );
684 : papszOffset2Tokens =
685 8 : CSLTokenizeStringComplex( pszOffset2, " ,", FALSE, FALSE );
686 :
687 8 : if( CSLCount(papszOffset1Tokens) >= 2
688 : && CSLCount(papszOffset2Tokens) >= 2
689 : && poOriginGeometry != NULL )
690 : {
691 8 : adfGeoTransform[0] = poOriginGeometry->getX();
692 8 : adfGeoTransform[1] = atof(papszOffset1Tokens[0]);
693 8 : adfGeoTransform[2] = atof(papszOffset2Tokens[0]);
694 8 : adfGeoTransform[3] = poOriginGeometry->getY();
695 8 : adfGeoTransform[4] = atof(papszOffset1Tokens[1]);
696 8 : adfGeoTransform[5] = atof(papszOffset2Tokens[1]);
697 :
698 : // offset from center of pixel.
699 8 : adfGeoTransform[0] -= adfGeoTransform[1]*0.5;
700 8 : adfGeoTransform[0] -= adfGeoTransform[2]*0.5;
701 8 : adfGeoTransform[3] -= adfGeoTransform[4]*0.5;
702 8 : adfGeoTransform[3] -= adfGeoTransform[5]*0.5;
703 :
704 8 : bSuccess = TRUE;
705 8 : bHaveGeoTransform = TRUE;
706 : }
707 :
708 8 : CSLDestroy( papszOffset1Tokens );
709 8 : CSLDestroy( papszOffset2Tokens );
710 :
711 8 : if( poOriginGeometry != NULL )
712 8 : delete poOriginGeometry;
713 :
714 : /* -------------------------------------------------------------------- */
715 : /* If we still don't have an srsName, check for it on the */
716 : /* boundedBy Envelope. Some products */
717 : /* (ie. EuropeRasterTile23.jpx) use this as the only srsName */
718 : /* delivery vehicle. */
719 : /* -------------------------------------------------------------------- */
720 8 : if( pszSRSName == NULL )
721 : {
722 : pszSRSName =
723 : CPLGetXMLValue( psXML,
724 : "=FeatureCollection.boundedBy.Envelope.srsName",
725 8 : NULL );
726 : }
727 :
728 : /* -------------------------------------------------------------------- */
729 : /* If we have gotten a geotransform, then try to interprete the */
730 : /* srsName. */
731 : /* -------------------------------------------------------------------- */
732 8 : int bNeedAxisFlip = FALSE;
733 :
734 8 : if( bSuccess && pszSRSName != NULL
735 : && (pszProjection == NULL || strlen(pszProjection) == 0) )
736 : {
737 8 : OGRSpatialReference oSRS;
738 :
739 8 : if( EQUALN(pszSRSName,"epsg:",5) )
740 : {
741 0 : if( oSRS.SetFromUserInput( pszSRSName ) == OGRERR_NONE )
742 0 : oSRS.exportToWkt( &pszProjection );
743 : }
744 8 : else if( EQUALN(pszSRSName,"urn:",4)
745 : && strstr(pszSRSName,":def:") != NULL
746 : && oSRS.importFromURN(pszSRSName) == OGRERR_NONE )
747 : {
748 8 : const char *pszCode = strrchr(pszSRSName,':') + 1;
749 :
750 8 : oSRS.exportToWkt( &pszProjection );
751 :
752 : // Per #2131
753 8 : if( atoi(pszCode) >= 4000 && atoi(pszCode) <= 4999 )
754 : {
755 : CPLDebug( "GMLJP2", "Request axis flip for SRS=%s",
756 8 : pszSRSName );
757 8 : bNeedAxisFlip = TRUE;
758 : }
759 : }
760 0 : else if( !GMLSRSLookup( pszSRSName ) )
761 : {
762 : CPLDebug( "GDALJP2Metadata",
763 : "Unable to evaluate SRSName=%s",
764 0 : pszSRSName );
765 8 : }
766 : }
767 :
768 8 : if( pszProjection )
769 : CPLDebug( "GDALJP2Metadata",
770 : "Got projection from GML box: %s",
771 8 : pszProjection );
772 :
773 8 : CPLDestroyXMLNode( psXML );
774 8 : psXML = NULL;
775 :
776 : /* -------------------------------------------------------------------- */
777 : /* Do we need to flip the axes? */
778 : /* -------------------------------------------------------------------- */
779 8 : if( bNeedAxisFlip
780 : && CSLTestBoolean( CPLGetConfigOption( "GDAL_IGNORE_AXIS_ORIENTATION",
781 : "FALSE" ) ) )
782 : {
783 0 : bNeedAxisFlip = FALSE;
784 0 : CPLDebug( "GMLJP2", "Supressed axis flipping based on GDAL_IGNORE_AXIS_ORIENTATION." );
785 : }
786 :
787 8 : if( bNeedAxisFlip )
788 : {
789 : double dfTemp;
790 :
791 : CPLDebug( "GMLJP2",
792 8 : "Flipping axis orientation in GMLJP2 coverage description." );
793 :
794 8 : dfTemp = adfGeoTransform[0];
795 8 : adfGeoTransform[0] = adfGeoTransform[3];
796 8 : adfGeoTransform[3] = dfTemp;
797 :
798 8 : dfTemp = adfGeoTransform[1];
799 8 : adfGeoTransform[1] = adfGeoTransform[4];
800 8 : adfGeoTransform[4] = dfTemp;
801 :
802 8 : dfTemp = adfGeoTransform[2];
803 8 : adfGeoTransform[2] = adfGeoTransform[5];
804 8 : adfGeoTransform[5] = dfTemp;
805 : }
806 :
807 8 : return pszProjection != NULL && bSuccess;
808 : }
809 :
810 : /************************************************************************/
811 : /* SetProjection() */
812 : /************************************************************************/
813 :
814 72 : void GDALJP2Metadata::SetProjection( const char *pszWKT )
815 :
816 : {
817 72 : CPLFree( pszProjection );
818 72 : pszProjection = CPLStrdup(pszWKT);
819 72 : }
820 :
821 : /************************************************************************/
822 : /* SetGCPs() */
823 : /************************************************************************/
824 :
825 48 : void GDALJP2Metadata::SetGCPs( int nCount, const GDAL_GCP *pasGCPsIn )
826 :
827 : {
828 48 : if( nGCPCount > 0 )
829 : {
830 0 : GDALDeinitGCPs( nGCPCount, pasGCPList );
831 0 : CPLFree( pasGCPList );
832 : }
833 :
834 48 : nGCPCount = nCount;
835 48 : pasGCPList = GDALDuplicateGCPs(nGCPCount, pasGCPsIn);
836 48 : }
837 :
838 : /************************************************************************/
839 : /* SetGeoTransform() */
840 : /************************************************************************/
841 :
842 72 : void GDALJP2Metadata::SetGeoTransform( double *padfGT )
843 :
844 : {
845 72 : memcpy( adfGeoTransform, padfGT, sizeof(double) * 6 );
846 72 : }
847 :
848 : /************************************************************************/
849 : /* CreateJP2GeoTIFF() */
850 : /************************************************************************/
851 :
852 66 : GDALJP2Box *GDALJP2Metadata::CreateJP2GeoTIFF()
853 :
854 : {
855 : /* -------------------------------------------------------------------- */
856 : /* Prepare the memory buffer containing the degenerate GeoTIFF */
857 : /* file. */
858 : /* -------------------------------------------------------------------- */
859 66 : int nGTBufSize = 0;
860 66 : unsigned char *pabyGTBuf = NULL;
861 :
862 66 : if( GTIFMemBufFromWkt( pszProjection, adfGeoTransform,
863 : nGCPCount, pasGCPList,
864 : &nGTBufSize, &pabyGTBuf ) != CE_None )
865 0 : return NULL;
866 :
867 66 : if( nGTBufSize == 0 )
868 0 : return NULL;
869 :
870 : /* -------------------------------------------------------------------- */
871 : /* Write to a box on the JP2 file. */
872 : /* -------------------------------------------------------------------- */
873 : GDALJP2Box *poBox;
874 :
875 66 : poBox = GDALJP2Box::CreateUUIDBox( msi_uuid2, nGTBufSize, pabyGTBuf );
876 :
877 66 : CPLFree( pabyGTBuf );
878 :
879 66 : return poBox;
880 : }
881 :
882 : /************************************************************************/
883 : /* PrepareCoverageBox() */
884 : /************************************************************************/
885 :
886 48 : GDALJP2Box *GDALJP2Metadata::CreateGMLJP2( int nXSize, int nYSize )
887 :
888 : {
889 : /* -------------------------------------------------------------------- */
890 : /* This is a backdoor to let us embed a literal gmljp2 chunk */
891 : /* supplied by the user as an external file. This is mostly */
892 : /* for preparing test files with exotic contents. */
893 : /* -------------------------------------------------------------------- */
894 48 : if( CPLGetConfigOption( "GMLJP2OVERRIDE", NULL ) != NULL )
895 : {
896 0 : VSILFILE *fp = VSIFOpenL( CPLGetConfigOption( "GMLJP2OVERRIDE",""), "r" );
897 0 : char *pszGML = NULL;
898 :
899 0 : if( fp == NULL )
900 : {
901 : CPLError( CE_Failure, CPLE_AppDefined,
902 0 : "Unable to open GMLJP2OVERRIDE file." );
903 0 : return NULL;
904 : }
905 :
906 0 : VSIFSeekL( fp, 0, SEEK_END );
907 0 : int nLength = (int) VSIFTellL( fp );
908 0 : pszGML = (char *) CPLCalloc(1,nLength+1);
909 0 : VSIFSeekL( fp, 0, SEEK_SET );
910 0 : VSIFReadL( pszGML, 1, nLength, fp );
911 0 : VSIFCloseL( fp );
912 :
913 : GDALJP2Box *apoGMLBoxes[2];
914 :
915 0 : apoGMLBoxes[0] = GDALJP2Box::CreateLblBox( "gml.data" );
916 : apoGMLBoxes[1] =
917 : GDALJP2Box::CreateLabelledXMLAssoc( "gml.root-instance",
918 0 : pszGML );
919 :
920 0 : GDALJP2Box *poGMLData = GDALJP2Box::CreateAsocBox( 2, apoGMLBoxes);
921 :
922 0 : delete apoGMLBoxes[0];
923 0 : delete apoGMLBoxes[1];
924 :
925 0 : CPLFree( pszGML );
926 :
927 0 : return poGMLData;
928 : }
929 :
930 : /* -------------------------------------------------------------------- */
931 : /* Try do determine a PCS or GCS code we can use. */
932 : /* -------------------------------------------------------------------- */
933 48 : OGRSpatialReference oSRS;
934 48 : char *pszWKTCopy = (char *) pszProjection;
935 48 : int nEPSGCode = 0;
936 : char szSRSName[100];
937 48 : int bNeedAxisFlip = FALSE;
938 :
939 48 : if( oSRS.importFromWkt( &pszWKTCopy ) != OGRERR_NONE )
940 4 : return NULL;
941 :
942 44 : if( oSRS.IsProjected() )
943 : {
944 12 : const char *pszAuthName = oSRS.GetAuthorityName( "PROJCS" );
945 :
946 12 : if( pszAuthName != NULL && EQUAL(pszAuthName,"epsg") )
947 : {
948 12 : nEPSGCode = atoi(oSRS.GetAuthorityCode( "PROJCS" ));
949 : }
950 : }
951 32 : else if( oSRS.IsGeographic() )
952 : {
953 32 : const char *pszAuthName = oSRS.GetAuthorityName( "GEOGCS" );
954 :
955 32 : if( pszAuthName != NULL && EQUAL(pszAuthName,"epsg") )
956 : {
957 24 : nEPSGCode = atoi(oSRS.GetAuthorityCode( "GEOGCS" ));
958 24 : bNeedAxisFlip = TRUE;
959 : }
960 : }
961 :
962 44 : if( nEPSGCode != 0 )
963 36 : sprintf( szSRSName, "urn:ogc:def:crs:EPSG::%d", nEPSGCode );
964 : else
965 : strcpy( szSRSName,
966 8 : "gmljp2://xml/CRSDictionary.gml#ogrcrs1" );
967 :
968 : /* -------------------------------------------------------------------- */
969 : /* Prepare coverage origin and offset vectors. Take axis */
970 : /* order into account if needed. */
971 : /* -------------------------------------------------------------------- */
972 : double adfOrigin[2];
973 : double adfXVector[2];
974 : double adfYVector[2];
975 :
976 44 : adfOrigin[0] = adfGeoTransform[0] + adfGeoTransform[1] * 0.5
977 44 : + adfGeoTransform[4] * 0.5;
978 44 : adfOrigin[1] = adfGeoTransform[3] + adfGeoTransform[2] * 0.5
979 44 : + adfGeoTransform[5] * 0.5;
980 44 : adfXVector[0] = adfGeoTransform[1];
981 44 : adfXVector[1] = adfGeoTransform[2];
982 :
983 44 : adfYVector[0] = adfGeoTransform[4];
984 44 : adfYVector[1] = adfGeoTransform[5];
985 :
986 44 : if( bNeedAxisFlip
987 : && CSLTestBoolean( CPLGetConfigOption( "GDAL_IGNORE_AXIS_ORIENTATION",
988 : "FALSE" ) ) )
989 : {
990 0 : bNeedAxisFlip = FALSE;
991 0 : CPLDebug( "GMLJP2", "Supressed axis flipping on write based on GDAL_IGNORE_AXIS_ORIENTATION." );
992 : }
993 :
994 44 : if( bNeedAxisFlip )
995 : {
996 : double dfTemp;
997 :
998 24 : CPLDebug( "GMLJP2", "Flipping GML coverage axis order." );
999 :
1000 24 : dfTemp = adfOrigin[0];
1001 24 : adfOrigin[0] = adfOrigin[1];
1002 24 : adfOrigin[1] = dfTemp;
1003 :
1004 24 : dfTemp = adfXVector[0];
1005 24 : adfXVector[0] = adfXVector[1];
1006 24 : adfXVector[1] = dfTemp;
1007 :
1008 24 : dfTemp = adfYVector[0];
1009 24 : adfYVector[0] = adfYVector[1];
1010 24 : adfYVector[1] = dfTemp;
1011 : }
1012 :
1013 : /* -------------------------------------------------------------------- */
1014 : /* For now we hardcode for a minimal instance format. */
1015 : /* -------------------------------------------------------------------- */
1016 44 : CPLString osDoc;
1017 :
1018 : osDoc.Printf(
1019 : "<gml:FeatureCollection\n"
1020 : " xmlns:gml=\"http://www.opengis.net/gml\"\n"
1021 : " xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"\n"
1022 : " xsi:schemaLocation=\"http://www.opengis.net/gml http://schemas.opengis.net/gml/3.1.1/profiles/gmlJP2Profile/1.0.0/gmlJP2Profile.xsd\">\n"
1023 : " <gml:boundedBy>\n"
1024 : " <gml:Null>withheld</gml:Null>\n"
1025 : " </gml:boundedBy>\n"
1026 : " <gml:featureMember>\n"
1027 : " <gml:FeatureCollection>\n"
1028 : " <gml:featureMember>\n"
1029 : " <gml:RectifiedGridCoverage dimension=\"2\" gml:id=\"RGC0001\">\n"
1030 : " <gml:rectifiedGridDomain>\n"
1031 : " <gml:RectifiedGrid dimension=\"2\">\n"
1032 : " <gml:limits>\n"
1033 : " <gml:GridEnvelope>\n"
1034 : " <gml:low>0 0</gml:low>\n"
1035 : " <gml:high>%d %d</gml:high>\n"
1036 : " </gml:GridEnvelope>\n"
1037 : " </gml:limits>\n"
1038 : " <gml:axisName>x</gml:axisName>\n"
1039 : " <gml:axisName>y</gml:axisName>\n"
1040 : " <gml:origin>\n"
1041 : " <gml:Point gml:id=\"P0001\" srsName=\"%s\">\n"
1042 : " <gml:pos>%.15g %.15g</gml:pos>\n"
1043 : " </gml:Point>\n"
1044 : " </gml:origin>\n"
1045 : " <gml:offsetVector srsName=\"%s\">%.15g %.15g</gml:offsetVector>\n"
1046 : " <gml:offsetVector srsName=\"%s\">%.15g %.15g</gml:offsetVector>\n"
1047 : " </gml:RectifiedGrid>\n"
1048 : " </gml:rectifiedGridDomain>\n"
1049 : " <gml:rangeSet>\n"
1050 : " <gml:File>\n"
1051 : " <gml:fileName>gmljp2://codestream/0</gml:fileName>\n"
1052 : " <gml:fileStructure>Record Interleaved</gml:fileStructure>\n"
1053 : " </gml:File>\n"
1054 : " </gml:rangeSet>\n"
1055 : " </gml:RectifiedGridCoverage>\n"
1056 : " </gml:featureMember>\n"
1057 : " </gml:FeatureCollection>\n"
1058 : " </gml:featureMember>\n"
1059 : "</gml:FeatureCollection>\n",
1060 : nXSize-1, nYSize-1, szSRSName, adfOrigin[0], adfOrigin[1],
1061 : szSRSName, adfXVector[0], adfXVector[1],
1062 44 : szSRSName, adfYVector[0], adfYVector[1] );
1063 :
1064 : /* -------------------------------------------------------------------- */
1065 : /* If we need a user defined CRSDictionary entry, prepare it */
1066 : /* here. */
1067 : /* -------------------------------------------------------------------- */
1068 44 : CPLString osDictBox;
1069 :
1070 44 : if( nEPSGCode == 0 )
1071 : {
1072 8 : char *pszGMLDef = NULL;
1073 :
1074 8 : if( oSRS.exportToXML( &pszGMLDef, NULL ) == OGRERR_NONE )
1075 : {
1076 : osDictBox.Printf(
1077 : "<gml:Dictionary gml:id=\"CRSU1\" \n"
1078 : " xmlns:gml=\"http://www.opengis.net/gml\"\n"
1079 : " xmlns:xlink=\"http://www.w3.org/1999/xlink\"\n"
1080 : " xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\">\n"
1081 : " <gml:dictionaryEntry>\n"
1082 : "%s\n"
1083 : " </gml:dictionaryEntry>\n"
1084 : "</gml:Dictionary>\n",
1085 8 : pszGMLDef );
1086 : }
1087 8 : CPLFree( pszGMLDef );
1088 : }
1089 :
1090 : /* -------------------------------------------------------------------- */
1091 : /* Setup the gml.data label. */
1092 : /* -------------------------------------------------------------------- */
1093 : GDALJP2Box *apoGMLBoxes[5];
1094 44 : int nGMLBoxes = 0;
1095 :
1096 44 : apoGMLBoxes[nGMLBoxes++] = GDALJP2Box::CreateLblBox( "gml.data" );
1097 :
1098 : /* -------------------------------------------------------------------- */
1099 : /* Setup gml.root-instance. */
1100 : /* -------------------------------------------------------------------- */
1101 88 : apoGMLBoxes[nGMLBoxes++] =
1102 44 : GDALJP2Box::CreateLabelledXMLAssoc( "gml.root-instance", osDoc );
1103 :
1104 : /* -------------------------------------------------------------------- */
1105 : /* Add optional dictionary. */
1106 : /* -------------------------------------------------------------------- */
1107 44 : if( strlen(osDictBox) > 0 )
1108 16 : apoGMLBoxes[nGMLBoxes++] =
1109 : GDALJP2Box::CreateLabelledXMLAssoc( "CRSDictionary.gml",
1110 8 : osDictBox );
1111 :
1112 : /* -------------------------------------------------------------------- */
1113 : /* Bundle gml.data boxes into an association. */
1114 : /* -------------------------------------------------------------------- */
1115 44 : GDALJP2Box *poGMLData = GDALJP2Box::CreateAsocBox( nGMLBoxes, apoGMLBoxes);
1116 :
1117 : /* -------------------------------------------------------------------- */
1118 : /* Cleanup working boxes. */
1119 : /* -------------------------------------------------------------------- */
1120 184 : while( nGMLBoxes > 0 )
1121 96 : delete apoGMLBoxes[--nGMLBoxes];
1122 :
1123 44 : return poGMLData;
1124 : }
1125 :
1126 :
|
