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