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